| 229 | #define P_HIGH 0.97575 |
| 230 | |
| 231 | static double normsinv(double p) { |
| 232 | if ((0 < p ) && (p < P_LOW)) { |
| 233 | const double q = sqrt(-2 * log(p)); |
| 234 | return (((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) / |
| 235 | ((((D1 * q + D2) * q + D3) * q + D4) * q + 1); |
| 236 | } else if ((P_LOW <= p) && (p <= P_HIGH)) { |
| 237 | const double q = p - 0.5; |
| 238 | const double r = q * q; |
| 239 | return (((((A1 * r + A2) * r + A3) * r + A4) * r + A5) * r + A6) * q / |
| 240 | (((((B1 * r + B2) * r + B3) * r + B4) * r + B5) * r + 1); |
| 241 | } else if ((P_HIGH < p) && (p < 1)) { |
| 242 | const double q = sqrt(-2 * log(1 - p)); |
| 243 | return -(((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6) / |
| 244 | ((((D1 * q + D2) * q + D3) * q + D4) * q + 1); |
| 245 | } else { |
| 246 | return INFINITY; |
| 247 | } |
| 248 | |
| 249 | /* Under UNIX OR LINUX, The return value 'x' could be corrected by the following for better accuracy. |
| 250 | if(( 0 < p)&&(p < 1)){ |
| 251 | e = 0.5 * erfc(-x/sqrt(2)) - p; |
| 252 | u = e * sqrt(2*M_PI) * exp(x*x/2); |
| 253 | x = x - u/(1 + x*u/2); |
| 254 | } |
| 255 | */ |
| 256 | |
| 257 | } |
| 258 | |
| 259 | static float normsinvf(float p) { |
| 260 | return static_cast<float>(normsinv(p)); |