Function gamma_p 

pub fn gamma_p(a: f64, x: f64) -> f64

The regularized lower incomplete gamma function $P(a, x) = \frac{1}{\Gamma(a)} \int_{0}^{x} u^{a-1} e^{-u}, du$, $a > 0$, $x \ge 0$: the distribution function of a gamma variable of shape $a$ and unit scale (a $\chi^{2}$ variable with $2a$ degrees of freedom, halved). The series $e^{-x} x^{a} \sum_{k} x^{k} / \Gamma(a + k + 1)$ for $x < a + 1$, and $1 - Q(a, x)$ with the continued fraction of $Q$ (modified Lentz) beyond (Press et al., Numerical Recipes, §6.2); relative error about $10^{-14}$ on the side computed directly.

Examples

use comnum_math::gamma_p;
// P(1, x) = 1 - e^{-x}
assert!((gamma_p(1.0, 0.5) - (1.0 - (-0.5f64).exp())).abs() < 1e-15);
assert!((gamma_p(20.0, 5.0) / 3.452_135_820_914_455e-7 - 1.0).abs() < 1e-12);