special_audit_dump/
special_audit_dump.rs

1//! Dumps `gam-math` special-function values as TSV for an external
2//! arbitrary-precision comparison. Not part of the build; run by hand:
3//!
4//! ```text
5//! cargo run -p gam-math --release --example special_audit_dump > /tmp/dump.tsv
6//! ```
7//!
8//! Each line is `channel<TAB>arg<TAB>value` (or `channel<TAB>arg<TAB>arg2<TAB>value`
9//! for the two-argument channels), with every float printed at full `{:e}`
10//! precision so the reader can round-trip it.
11
12use gam_math::probability as prob;
13use gam_math::special;
14
15fn emit(channel: &str, arg: f64, value: f64) {
16    println!("{channel}\t{arg:e}\t{value:e}");
17}
18
19fn geometric_grid(lo: f64, hi: f64, count: usize) -> Vec<f64> {
20    let log_lo = lo.ln();
21    let log_hi = hi.ln();
22    (0..count)
23        .map(|i| (log_lo + (log_hi - log_lo) * (i as f64) / ((count - 1) as f64)).exp())
24        .collect()
25}
26
27fn linear_grid(lo: f64, hi: f64, count: usize) -> Vec<f64> {
28    (0..count)
29        .map(|i| lo + (hi - lo) * (i as f64) / ((count - 1) as f64))
30        .collect()
31}
32
33fn main() {
34    // ---- Bessel channels -------------------------------------------------
35    for eta in geometric_grid(1e-6, 1e12, 400) {
36        let (centered, ratio, d1) = special::bessel_i0_centered_terms(eta);
37        emit("bessel_centered_log", eta, centered);
38        emit("bessel_ratio", eta, ratio);
39        emit("bessel_d1", eta, d1);
40        emit(
41            "bessel_d2",
42            eta,
43            special::bessel_i0_centered_second_log_derivative_from_log_abs(eta.ln()),
44        );
45    }
46    // Dense sweep across the ascending/asymptotic seam at 20.
47    for eta in linear_grid(15.0, 25.0, 201) {
48        let (centered, ratio, d1) = special::bessel_i0_centered_terms(eta);
49        emit("bessel_centered_log", eta, centered);
50        emit("bessel_ratio", eta, ratio);
51        emit("bessel_d1", eta, d1);
52        emit(
53            "bessel_d2",
54            eta,
55            special::bessel_i0_centered_second_log_derivative_from_log_abs(eta.ln()),
56        );
57    }
58
59    // ---- Polygamma family ------------------------------------------------
60    for x in geometric_grid(1e-8, 1e10, 400) {
61        emit("digamma", x, special::digamma(x));
62        emit("trigamma", x, special::trigamma(x));
63        emit("tetragamma", x, special::tetragamma(x));
64        emit("pentagamma", x, special::pentagamma(x));
65    }
66    for x in linear_grid(0.5, 30.0, 300) {
67        emit("digamma", x, special::digamma(x));
68        emit("trigamma", x, special::trigamma(x));
69        emit("tetragamma", x, special::tetragamma(x));
70        emit("pentagamma", x, special::pentagamma(x));
71    }
72
73    // ---- Normal-distribution channels ------------------------------------
74    let mut normal_args: Vec<f64> = linear_grid(-40.0, 40.0, 801);
75    normal_args.extend(geometric_grid(1e-8, 1e2, 200));
76    normal_args.extend(geometric_grid(1e-8, 1e2, 200).into_iter().map(|v| -v));
77    for x in normal_args {
78        emit("normal_pdf", x, prob::normal_pdf(x));
79        emit("normal_cdf", x, prob::normal_cdf(x));
80        emit("normal_logcdf", x, prob::normal_logcdf(x));
81        emit("normal_logsf", x, prob::normal_logsf(x));
82        let (log_cdf, mills) = prob::signed_probit_logcdf_and_mills_ratio(x);
83        emit("probit_logcdf", x, log_cdf);
84        emit("probit_mills", x, mills);
85        let derivatives = prob::normal_logcdf_derivatives(x);
86        for (order, value) in derivatives.iter().enumerate() {
87            emit(&format!("logcdf_d{order}"), x, *value);
88        }
89        if x >= 0.0 {
90            emit("erfcx", x, prob::erfcx_nonnegative(x));
91        }
92    }
93    for x in geometric_grid(1e-3, 1e6, 400) {
94        emit("erfcx", x, prob::erfcx_nonnegative(x));
95    }
96
97    // ---- log1mexp --------------------------------------------------------
98    for a in geometric_grid(1e-14, 50.0, 300) {
99        emit("log1mexp", a, prob::log1mexp_positive(a));
100    }
101
102    // ---- Normal quantile -------------------------------------------------
103    for p in geometric_grid(1e-300, 0.5, 400) {
104        if let Ok(q) = prob::standard_normal_quantile(p) {
105            emit("normal_quantile", p, q);
106        }
107    }
108    for p in linear_grid(0.001, 0.999, 400) {
109        if let Ok(q) = prob::standard_normal_quantile(p) {
110            emit("normal_quantile", p, q);
111        }
112    }
113    for log_p in linear_grid(-700.0, -1e-6, 400) {
114        if let Ok(q) = prob::standard_normal_quantile_from_log_cdf(log_p) {
115            emit("normal_quantile_from_log", log_p, q);
116        }
117    }
118
119    // ---- Gauss-Legendre --------------------------------------------------
120    for n in [
121        3usize, 4, 5, 7, 8, 12, 15, 16, 20, 24, 31, 32, 40, 48, 63, 64, 80, 96, 100, 127, 128, 160,
122        200, 256,
123    ] {
124        let (nodes, weights) = special::gauss_legendre(n);
125        for (i, (node, weight)) in nodes.iter().zip(weights.iter()).enumerate() {
126            println!("gl_node\t{n}\t{i}\t{node:e}");
127            println!("gl_weight\t{n}\t{i}\t{weight:e}");
128        }
129    }
130
131    // ---- Binomial coefficient --------------------------------------------
132    for n in 0usize..=60 {
133        for k in 0..=n {
134            println!(
135                "binomial\t{n}\t{k}\t{:e}",
136                special::binomial_coefficient_f64(n, k)
137            );
138        }
139    }
140}
special_audit_dump/special_audit_dump.rs

special_audit_dump/
special_audit_dump.rs
