eic-opticks/sysrap/njuffa_erfcinvf.h

0001 /**
0002 njuffa_erfcinvf.h
0003 ====================
0004
0005 https://stackoverflow.com/questions/60472139/computing-the-inverse-of-the-complementary-error-function-erfcinv-in-c
0006
0007 This C CPU implementation by njuffa (from stackoverflow) of the inverse of the complementary
0008 error function gives results very close to those from the CUDA erfcinvf function.
0009
0010 **/
0011
0012
0013 #include <math.h>
0014 #include <stdint.h>
0015
0016 #define PORTABLE  (1)
0017
0018 float njuffa_logf (float a);
0019 #if !PORTABLE
0020 #include "immintrin.h"
0021 float sse_recipf (float a);
0022 float sse_rsqrtf (float a);
0023 #endif // !PORTABLE
0024
0025 /* Compute inverse of the complementary error function. For the central region,
0026    re-use the polynomial approximation for erfinv. For the tail regions, use an
0027    approximation based on the observation that erfcinv(x) is very approximately
0028    sqrt(-log(x)).
0029
0030    PORTABLE=1 max. ulp err. = 3.12017
0031    PORTABLE=0 max. ulp err. = 3.13523
0032 */
0033 float njuffa_erfcinvf (float a)
0034 {
0035     float r;
0036
0037     if ((a >= 2.1875e-3f) && (a <= 1.998125f)) { // max. ulp err. = 2.77667
0038         float p, t;
0039         t = fmaf (-a, a, a + a);
0040         t = njuffa_logf (t);
0041         p =              5.43877832e-9f;  //  0x1.75c000p-28
0042         p = fmaf (p, t,  1.43286059e-7f); //  0x1.33b458p-23
0043         p = fmaf (p, t,  1.22775396e-6f); //  0x1.49929cp-20
0044         p = fmaf (p, t,  1.12962631e-7f); //  0x1.e52bbap-24
0045         p = fmaf (p, t, -5.61531961e-5f); // -0x1.d70c12p-15
0046         p = fmaf (p, t, -1.47697705e-4f); // -0x1.35be9ap-13
0047         p = fmaf (p, t,  2.31468701e-3f); //  0x1.2f6402p-9
0048         p = fmaf (p, t,  1.15392562e-2f); //  0x1.7a1e4cp-7
0049         p = fmaf (p, t, -2.32015476e-1f); // -0x1.db2aeep-3
0050         t = fmaf (p, t,  8.86226892e-1f); //  0x1.c5bf88p-1
0051         r = fmaf (t, -a, t);
0052     } else {
0053         float p, q, s, t;
0054         t = (a >= 1.0f) ? (2.0f - a) : a;
0055         t = 0.0f - njuffa_logf (t);
0056 #if PORTABLE
0057         s = sqrtf (1.0f / t);
0058 #else // PORTABLE
0059         s = sse_rsqrtf (t);
0060 #endif // PORTABLE
0061         p =              2.23100796e+1f;  //  0x1.64f616p+4
0062         p = fmaf (p, s, -5.23008537e+1f); // -0x1.a26826p+5
0063         p = fmaf (p, s,  5.44409714e+1f); //  0x1.b3871cp+5
0064         p = fmaf (p, s, -3.35030403e+1f); // -0x1.0c063ap+5
0065         p = fmaf (p, s,  1.38580027e+1f); //  0x1.bb74c2p+3
0066         p = fmaf (p, s, -4.37277269e+0f); // -0x1.17db82p+2
0067         p = fmaf (p, s,  1.53075826e+0f); //  0x1.87dfc6p+0
0068         p = fmaf (p, s,  2.97993328e-2f); //  0x1.e83b76p-6
0069         p = fmaf (p, s, -3.71997419e-4f); // -0x1.86114cp-12
0070         p = fmaf (p, s, s);
0071 #if PORTABLE
0072         r = 1.0f / p;
0073 #else // PORTABLE
0074         r = sse_recipf (p);
0075         if (t == INFINITY) r = t;
0076 #endif // PORTABLE
0077         if (a >= 1.0f) r = 0.0f - r;
0078     }
0079     return r;
0080 }
0081
0082 /* Compute inverse of the CDF of the standard normal distribution.
0083    max ulp err = 4.08385
0084 */
0085 float njuffa_normcdfinvf (float a)
0086 {
0087     return fmaf (-1.41421356f, njuffa_erfcinvf (a + a), 0.0f);
0088 }
0089
0090 /* natural logarithm. max ulp err = 0.85089 */
0091 float njuffa_logf (float a)
0092 {
0093     float m, r, s, t, i, f;
0094     int32_t e;
0095     const float cutoff_f = 0.666666667f;
0096     if ((a > 0.0f) && (a <=  0x1.fffffep+127f)) { // 3.40282347e+38
0097         m = frexpf (a, &e);
0098         if (m < cutoff_f) {
0099             m = m + m;
0100             e = e - 1;
0101         }
0102         i = (float)e;
0103         f = m - 1.0f;
0104         s = f * f;
0105         /* Compute log1p(f) for f in [-1/3, 1/3] */
0106         r =             -0x1.0ae000p-3f;  // -0.130310059
0107         t =              0x1.208000p-3f;  //  0.140869141
0108         r = fmaf (r, s, -0x1.f1988ap-4f); // -0.121483363
0109         t = fmaf (t, s,  0x1.1e5740p-3f); //  0.139814854
0110         r = fmaf (r, s, -0x1.55b36ep-3f); // -0.166846141
0111         t = fmaf (t, s,  0x1.99d8b2p-3f); //  0.200120345
0112         r = fmaf (r, s, -0x1.fffe02p-3f); // -0.249996200
0113         r = fmaf (t, f, r);
0114         r = fmaf (r, f,  0x1.5554fap-2f); //  0.333331972
0115         r = fmaf (r, f, -0x1.000000p-1f); // -0.500000000
0116         r = fmaf (r, s, f);
0117         r = fmaf (i, 0x1.62e430p-01f, r); //  0.693147182 // log(2)
0118     } else {
0119         r = a + a;  // silence NaNs if necessary
0120         if (a  < 0.0f) r = 0.0f/0.0f; // QNaN INDEFINITE
0121         if (a == 0.0f) r = -INFINITY; // -INF
0122     }
0123     return r;
0124 }
0125
0126 #if !PORTABLE
0127 float sse_recipf (float a)
0128 {
0129     __m128 t;
0130     float e, r;
0131     t = _mm_set_ss (a);
0132     t = _mm_rcp_ss (t);
0133     _mm_store_ss (&r, t);
0134     e = fmaf (0.0f - a, r, 1.0f);
0135     e = fmaf (e, e, e);
0136     r = fmaf (e, r, r);
0137     return r;
0138 }
0139
0140
0141 float sse_rsqrtf (float a)
0142 {
0143     __m128 t;
0144     float e, r;
0145     t = _mm_set_ss (a);
0146     t = _mm_rsqrt_ss (t);
0147     _mm_store_ss (&r, t);
0148     e = fmaf (0.0f - a, r * r, 1.0f);
0149     r = fmaf (fmaf (0.375f, e, 0.5f), e * r, r);
0150     return r;
0151 }
0152 #endif
0153