use mathr::special;
fn main() {
println!("Gamma function:");
for x in [0.5, 1.0, 1.5, 2.0, 3.0, 5.0, 10.0, 0.1] {
println!(" Γ({:5.1}) = {:.10}", x, special::gamma(x));
}
let sqrt_pi = std::f64::consts::PI.sqrt();
println!(" Γ(0.5) = √π = {:.10} (error: {:.2e})", special::gamma(0.5), (special::gamma(0.5) - sqrt_pi).abs());
println!("\nlog-Gamma:");
for x in [1.0, 2.0, 10.0, 100.0] {
println!(" ln Γ({:5.1}) = {:.10}", x, special::log_gamma(x));
}
println!("\nBeta function:");
println!(" B(1, 1) = {:.10} (expected 1.0)", special::beta(1.0, 1.0));
println!(" B(0.5, 0.5) = {:.10} (expected π)", special::beta(0.5, 0.5));
println!(" B(2, 3) = {:.10} (expected 0.0833...)", special::beta(2.0, 3.0));
println!("\nError function:");
for x in [-2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 3.0] {
println!(" erf({:5.1}) = {:12.10} erfc({:5.1}) = {:12.10}",
x, special::erf(x), x, special::erfc(x));
}
println!("\nSinc function:");
for x in [0.0, 0.5, 1.0, 2.0, 3.0] {
println!(" sinc({:5.1}) = {:.10}", x, special::sinc(x));
}
println!("\nIncomplete gamma P(2, x) = χ² CDF (2 dof):");
for x in [1.0, 2.0, 4.0, 6.0, 10.0] {
let p = special::incomplete_gamma_p(1.0, x / 2.0);
println!(" P(1, {:4.1}) = {:.10}", x / 2.0, p);
}
println!("\nZeta, polygamma, harmonic:");
println!(" zeta(2) = {:.12} (= π²/6)", special::zeta(2.0));
println!(" zeta(-1) = {:.12} (= −1/12)", special::zeta(-1.0));
println!(" zeta(3) = {:.12} (Apéry)", special::zeta(3.0));
println!(" harmonic(10) = {:.12}", special::harmonic(10));
println!(" polygamma(2, 1) = {:.12} (= −2·ζ(3))", special::polygamma(2, 1.0));
println!("\nElliptic integrals (modulus k):");
for k in [0.0, 0.5, 0.9] {
println!(" K({}) = {:.12} E({}) = {:.12}", k, special::elliptic_k(k), k, special::elliptic_e(k));
}
}