Incomplete gamma integral The function is defined by x - 1 | | -t a-1 igam(a,x) = ----- | e t dt. - | | | (a) - 0 In this implementation both arguments must be positive. The integral is evaluated by either a power series or continued fr
| 57 | Copyright 1985, 1987, 2000 by Stephen L. Moshier |
| 58 | *************************************************************************/ |
| 59 | double incompletegamma(double a, double x) |
| 60 | { |
| 61 | double result; |
| 62 | double igammaepsilon; |
| 63 | double ans; |
| 64 | double ax; |
| 65 | double c; |
| 66 | double r; |
| 67 | double tmp; |
| 68 | |
| 69 | igammaepsilon = 0.000000000000001; |
| 70 | if( ap::fp_less_eq(x,0)||ap::fp_less_eq(a,0) ) |
| 71 | { |
| 72 | result = 0; |
| 73 | return result; |
| 74 | } |
| 75 | if( ap::fp_greater(x,1)&&ap::fp_greater(x,a) ) |
| 76 | { |
| 77 | result = 1-incompletegammac(a, x); |
| 78 | return result; |
| 79 | } |
| 80 | ax = a*log(x)-x-lngamma(a, tmp); |
| 81 | if( ap::fp_less(ax,-709.78271289338399) ) |
| 82 | { |
| 83 | result = 0; |
| 84 | return result; |
| 85 | } |
| 86 | ax = exp(ax); |
| 87 | r = a; |
| 88 | c = 1; |
| 89 | ans = 1; |
| 90 | do |
| 91 | { |
| 92 | r = r+1; |
| 93 | c = c*x/r; |
| 94 | ans = ans+c; |
| 95 | } |
| 96 | while(ap::fp_greater(c/ans,igammaepsilon)); |
| 97 | result = ans*ax/a; |
| 98 | return result; |
| 99 | } |
| 100 | |
| 101 | |
| 102 | /************************************************************************* |
no test coverage detected