Inverse of imcomplete beta integral Given y, the function finds x such that incbet( a, b, x ) = y . The routine performs interval halving or Newton iterations to find the root of incbet(a,b,x) - y = 0. ACCURACY: Relative error: x a,b arithmetic domain domain # trials peak rms IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13 IE
| 236 | Copyright 1984, 1996, 2000 by Stephen L. Moshier |
| 237 | *************************************************************************/ |
| 238 | double invincompletebeta(double a, double b, double y) |
| 239 | { |
| 240 | double result; |
| 241 | double aaa; |
| 242 | double bbb; |
| 243 | double y0; |
| 244 | double d; |
| 245 | double yyy; |
| 246 | double x; |
| 247 | double x0; |
| 248 | double x1; |
| 249 | double lgm; |
| 250 | double yp; |
| 251 | double di; |
| 252 | double dithresh; |
| 253 | double yl; |
| 254 | double yh; |
| 255 | double xt; |
| 256 | int i; |
| 257 | int rflg; |
| 258 | int dir; |
| 259 | int nflg; |
| 260 | double s; |
| 261 | int mainlooppos; |
| 262 | int ihalve; |
| 263 | int ihalvecycle; |
| 264 | int newt; |
| 265 | int newtcycle; |
| 266 | int breaknewtcycle; |
| 267 | int breakihalvecycle; |
| 268 | |
| 269 | i = 0; |
| 270 | ap::ap_error::make_assertion(ap::fp_greater_eq(y,0)&&ap::fp_less_eq(y,1), "Domain error in InvIncompleteBeta"); |
| 271 | if( ap::fp_eq(y,0) ) |
| 272 | { |
| 273 | result = 0; |
| 274 | return result; |
| 275 | } |
| 276 | if( ap::fp_eq(y,1.0) ) |
| 277 | { |
| 278 | result = 1; |
| 279 | return result; |
| 280 | } |
| 281 | x0 = 0.0; |
| 282 | yl = 0.0; |
| 283 | x1 = 1.0; |
| 284 | yh = 1.0; |
| 285 | nflg = 0; |
| 286 | mainlooppos = 0; |
| 287 | ihalve = 1; |
| 288 | ihalvecycle = 2; |
| 289 | newt = 3; |
| 290 | newtcycle = 4; |
| 291 | breaknewtcycle = 5; |
| 292 | breakihalvecycle = 6; |
| 293 | while(true) |
| 294 | { |
| 295 |
no test coverage detected