Inverse error functions

Functions
real	inv_erf (real x)
real	inv_erf_nsqrt (real x)

Detailed Description

 
real inv_erf(real x) {
  if (is_nan(x) || is_inf(x) || x < 0 || x > 1) 
    reject("inv_erf: x must be finite and between 0 and 1; ", "found x = ", x);
  
  real w = -log1m(square(x));
  real s = sqrt(w);
  real p;
  
  if (w < 5.0) {
    w = w - 2.5;
    p = 2.81022636e-08;
    p = fma(p, w, 3.43273939e-07);
    p = fma(p, w, -3.5233877e-06);
    p = fma(p, w, -4.39150654e-06);
    p = fma(p, w, 0.00021858087);
    p = fma(p, w, -0.00125372503);
    p = fma(p, w, -0.00417768164);
    p = fma(p, w, 0.246640727);
    p = fma(p, w, 1.50140941);
  } else {
    w = s - 3.000000;
    p = -0.000200214257;
    p = fma(p, w, 0.000100950558);
    p = fma(p, w, 0.00134934322);
    p = fma(p, w, -0.00367342844);
    p = fma(p, w, 0.00573950773);
    p = fma(p, w, -0.0076224613);
    p = fma(p, w, 0.00943887047);
    p = fma(p, w, 1.00167406);
    p = fma(p, w, 2.83297682);
  }
  return p * x;
}
 
real inv_erf_nsqrt(real x) {
  if (is_nan(x) || is_inf(x) || x < 0 || x > 1) 
    reject("inv_erf_nsqrt: x must be finite and between 0 and 1; ", "found x = ", x);
  
  real p;
  real t;
  t = fma(x, -x, 1.0);
  t = log(t);
  if (fabs(t) > 6.125) {
    // maximum ulp error = 2.35793
    p = 3.03697567e-10; //  0x1.4deb44p-32 
    p = fma(p, t, 2.93243101e-8); //  0x1.f7c9aep-26 
    p = fma(p, t, 1.22150334e-6); //  0x1.47e512p-20 
    p = fma(p, t, 2.84108955e-5); //  0x1.dca7dep-16 
    p = fma(p, t, 3.93552968e-4); //  0x1.9cab92p-12 
    p = fma(p, t, 3.02698812e-3); //  0x1.8cc0dep-9 
    p = fma(p, t, 4.83185798e-3); //  0x1.3ca920p-8 
    p = fma(p, t, -2.64646143e-1); // -0x1.0eff66p-2 
    p = fma(p, t, 8.40016484e-1); //  0x1.ae16a4p-1 
  } else {
    // maximum ulp error = 2.35002
    p = 5.43877832e-9; //  0x1.75c000p-28 
    p = fma(p, t, 1.43285448e-7); //  0x1.33b402p-23 
    p = fma(p, t, 1.22774793e-6); //  0x1.499232p-20 
    p = fma(p, t, 1.12963626e-7); //  0x1.e52cd2p-24 
    p = fma(p, t, -5.61530760e-5); // -0x1.d70bd0p-15 
    p = fma(p, t, -1.47697632e-4); // -0x1.35be90p-13 
    p = fma(p, t, 2.31468678e-3); //  0x1.2f6400p-9 
    p = fma(p, t, 1.15392581e-2); //  0x1.7a1e50p-7 
    p = fma(p, t, -2.32015476e-1); // -0x1.db2aeep-3 
    p = fma(p, t, 8.86226892e-1); //  0x1.c5bf88p-1 
  }
  return p * x;
}

Function Documentation

inv_erf()

real inv_erf

(

real

x

)

Inverse error function

Copyright Sean Pinkney, Feb. 2021

Giles, Mike. (2012). "Approximating the Erfinv Function." GPU Computing Gems Jade Edition. 10.1016/B978-0-12-385963-1.00010-1. accessed Feb. 15, 2021. at https://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf

Parameters

x	Real number on [0,1]
return	inverse error result

inv_erf_nsqrt()

real inv_erf_nsqrt

(

real

x

)

Inverse error function no sqrt

Copyright njuffa, Sean Pinkney, 2021

Compute the inverse error functions with maximum error of 2.35793 ul Originally implemented by njuffa and accessed June 11, 2021 at https://stackoverflow.com/a/49743348

Parameters

x	Real number on [0,1]
return	inverse error result