The inverse erf and erfc functions share a common implementation, this version is for 80-bit long double's and smaller:
| 72 | // this version is for 80-bit long double's and smaller: |
| 73 | // |
| 74 | inline double erfinv_imp(double p, double q) { |
| 75 | using namespace std; |
| 76 | |
| 77 | double result = 0; |
| 78 | |
| 79 | if (p <= 0.5) { |
| 80 | // |
| 81 | // Evaluate inverse erf using the rational approximation: |
| 82 | // |
| 83 | // x = p(p+10)(Y+R(p)) |
| 84 | // |
| 85 | // Where Y is a constant, and R(p) is optimised for a low |
| 86 | // absolute error compared to |Y|. |
| 87 | // |
| 88 | // double: Max error found: 2.001849e-18 |
| 89 | // long double: Max error found: 1.017064e-20 |
| 90 | // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21 |
| 91 | // |
| 92 | static const float Y = 0.0891314744949340820313f; |
| 93 | static const double P[] = { |
| 94 | -0.000508781949658280665617, -0.00836874819741736770379, |
| 95 | 0.0334806625409744615033, -0.0126926147662974029034, |
| 96 | -0.0365637971411762664006, 0.0219878681111168899165, |
| 97 | 0.00822687874676915743155, -0.00538772965071242932965}; |
| 98 | static const double Q[] = { |
| 99 | 1.0, |
| 100 | -0.970005043303290640362, |
| 101 | -1.56574558234175846809, |
| 102 | 1.56221558398423026363, |
| 103 | 0.662328840472002992063, |
| 104 | -0.71228902341542847553, |
| 105 | -0.0527396382340099713954, |
| 106 | 0.0795283687341571680018, |
| 107 | -0.00233393759374190016776, |
| 108 | 0.000886216390456424707504}; |
| 109 | double g = p * (p + 10); |
| 110 | double r = evaluate_polynomial(P, p) / evaluate_polynomial(Q, p); |
| 111 | result = g * Y + g * r; |
| 112 | } else if (q >= 0.25) { |
| 113 | // |
| 114 | // Rational approximation for 0.5 > q >= 0.25 |
| 115 | // |
| 116 | // x = sqrt(-2*log(q)) / (Y + R(q)) |
| 117 | // |
| 118 | // Where Y is a constant, and R(q) is optimised for a low |
| 119 | // absolute error compared to Y. |
| 120 | // |
| 121 | // double : Max error found: 7.403372e-17 |
| 122 | // long double : Max error found: 6.084616e-20 |
| 123 | // Maximum Deviation Found (error term) 4.811e-20 |
| 124 | // |
| 125 | static const float Y = 2.249481201171875f; |
| 126 | static const double P[] = {-0.202433508355938759655, 0.105264680699391713268, |
| 127 | 8.37050328343119927838, 17.6447298408374015486, |
| 128 | -18.8510648058714251895, -44.6382324441786960818, |
| 129 | 17.445385985570866523, 21.1294655448340526258, |
| 130 | -3.67192254707729348546}; |
| 131 | static const double Q[] = { |
no test coverage detected