Incomplete beta integral Returns incomplete beta integral of the arguments, evaluated from zero to x. The function is defined as x - - | (a+b) | | a-1 b-1 ----------- | t (1-t) dt. - - | | | (a) | (b) - 0 The domain of definition is 0 <= x <= 1. In this implementation a and b are restricted to positive values.
| 82 | Copyright 1984, 1995, 2000 by Stephen L. Moshier |
| 83 | *************************************************************************/ |
| 84 | double incompletebeta(double a, double b, double x) |
| 85 | { |
| 86 | double result; |
| 87 | double t; |
| 88 | double xc; |
| 89 | double w; |
| 90 | double y; |
| 91 | int flag; |
| 92 | double sg; |
| 93 | double big; |
| 94 | double biginv; |
| 95 | double maxgam; |
| 96 | double minlog; |
| 97 | double maxlog; |
| 98 | |
| 99 | big = 4.503599627370496e15; |
| 100 | biginv = 2.22044604925031308085e-16; |
| 101 | maxgam = 171.624376956302725; |
| 102 | minlog = log(ap::minrealnumber); |
| 103 | maxlog = log(ap::maxrealnumber); |
| 104 | ap::ap_error::make_assertion(ap::fp_greater(a,0)&&ap::fp_greater(b,0), "Domain error in IncompleteBeta"); |
| 105 | ap::ap_error::make_assertion(ap::fp_greater_eq(x,0)&&ap::fp_less_eq(x,1), "Domain error in IncompleteBeta"); |
| 106 | if( ap::fp_eq(x,0) ) |
| 107 | { |
| 108 | result = 0; |
| 109 | return result; |
| 110 | } |
| 111 | if( ap::fp_eq(x,1) ) |
| 112 | { |
| 113 | result = 1; |
| 114 | return result; |
| 115 | } |
| 116 | flag = 0; |
| 117 | if( ap::fp_less_eq(b*x,1.0)&&ap::fp_less_eq(x,0.95) ) |
| 118 | { |
| 119 | result = incompletebetaps(a, b, x, maxgam); |
| 120 | return result; |
| 121 | } |
| 122 | w = 1.0-x; |
| 123 | if( ap::fp_greater(x,a/(a+b)) ) |
| 124 | { |
| 125 | flag = 1; |
| 126 | t = a; |
| 127 | a = b; |
| 128 | b = t; |
| 129 | xc = x; |
| 130 | x = w; |
| 131 | } |
| 132 | else |
| 133 | { |
| 134 | xc = w; |
| 135 | } |
| 136 | if( flag==1&&ap::fp_less_eq(b*x,1.0)&&ap::fp_less_eq(x,0.95) ) |
| 137 | { |
| 138 | t = incompletebetaps(a, b, x, maxgam); |
| 139 | if( ap::fp_less_eq(t,ap::machineepsilon) ) |
| 140 | { |
| 141 | result = 1.0-ap::machineepsilon; |
no test coverage detected