| 301 | } |
| 302 | |
| 303 | inline double erfinv(double z) { |
| 304 | // |
| 305 | // Begin by testing for domain errors, and other special cases: |
| 306 | // |
| 307 | if ((z < -1) || (z > 1)) |
| 308 | return NAN; |
| 309 | if (z == 1) |
| 310 | return INFINITY; |
| 311 | if (z == -1) |
| 312 | return -INFINITY; |
| 313 | if (z == 0) |
| 314 | return 0; |
| 315 | // |
| 316 | // Normalise the input, so it's in the range [0,1], we will |
| 317 | // negate the result if z is outside that range. This is a simple |
| 318 | // application of the erf reflection formula: erf(-z) = -erf(z) |
| 319 | // |
| 320 | double p, q, s; |
| 321 | if (z < 0) { |
| 322 | p = -z; |
| 323 | q = 1 - p; |
| 324 | s = -1; |
| 325 | } else { |
| 326 | p = z; |
| 327 | q = 1 - z; |
| 328 | s = 1; |
| 329 | } |
| 330 | |
| 331 | // |
| 332 | // And get the result, negating where required: |
| 333 | // |
| 334 | return s * erfinv_imp(p, q); |
| 335 | } |
| 336 | |
| 337 | inline float erfcinvf(float z) { |
| 338 | return erfcinv(z); |