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);
}
}