Computes the regularized gamma function Q(a,z) = Γ(a,z)/Γ(a). This method is more numerically stable and accurate than computing it via the direct method, and is always in the range [0,1]. Note: The this method returns Double#NaN for a<0, though real values of Q(a,z) do exis
(double a, double z)
| 636 | * @return Q(a,z) |
| 637 | */ |
| 638 | public static double gammaQ(double a, double z) |
| 639 | { |
| 640 | if(z<= 0 || a < 0 ) |
| 641 | return Double.NaN; |
| 642 | |
| 643 | if(z < a+1) |
| 644 | return 1-gammaPSeries(a, z); |
| 645 | |
| 646 | /** |
| 647 | * On the range of x from 0.5 to 50 |
| 648 | * a=0.15, |rel error| is ~ 3e-15 for most values of x, with a bad spot of |rel error| ~ 3e-11 when x ~= 5.75 |
| 649 | * a=0.5, max |rel error| is 3.9e-15 |
| 650 | * a=1, max |rel error| is 4e-15, rel error groining as x -> infinity |
| 651 | * a=5, max |rel error| is 7e-13, but only near x ~= 0, most is in the range 3e-15 |
| 652 | * a=10, max |rel error| is 4e-6, but only near x ~= 0, most is in the range 3e-15 |
| 653 | */ |
| 654 | return exp(a*log(z)-z-lnGamma(a))/gammaQ.lentz(a, z); |
| 655 | } |
| 656 | |
| 657 | public static double gammaPSeries(double a, double z) |
| 658 | { |