Inverse Error Function in C

Question

Is it possible to calculate the inverse error function in C?

I can find erf(x) in <math.h> which calculates the error function, but I can't find anything to do the inverse.

Answer 1

I don't think it's a standard implementation in <math.h>, but there are other C math libraries that have implement the inverse error function erfinv(x), that you can use.

Answer 2

I wrote another method that uses the fast-converging Newton-Rhapson method, which is an iterative method to find the root of a function. It starts with an initial guess and then iteratively improves the guess by using the derivative of the function. The Newton-Raphson method requires the function, its derivative, an initial guess and a stopping criteria.

In this case, the function we are trying to find the root of is erf(x) - x. And the derivative of this function is 2.0 / sqrt(pi) * exp(-x**2). The initial guess is the input value for x. The stopping criteria is a tolerance value, in this case it's 1.0e-16. Here is the code:

/*
============================================
 Compile and execute with:
    $ gcc inverf.c -o inverf -lm
    $ ./inverf
============================================
*/
#include <stdio.h>
#include <math.h>

int main() {
    double x, result, fx, dfx, dx, xold;
    double tolerance = 1.0e-16;
    double pi = 4.0 * atan(1.0);
    int iteration, i;

    // input value for x
    printf("Calculator for inverse error function.\n");
    printf("Enter the value for x: ");
    scanf("%lf", &x);

    // check the input value is between -1 and 1
    if (x < -1.0 || x > 1.0) {
        printf("Invalid input, x must be between -1 and 1.");
        return 0;
    }

    // initial guess
    result = x;
    xold = 0.0;
    iteration = 0;

    // iterate until the solution converges
    do {
        xold = result;
        fx = erf(result) - x;
        dfx = 2.0 / sqrt(pi) * exp(-pow(result, 2.0));
        dx = fx / dfx;

        // update the solution
        result = result - dx;

        iteration = iteration + 1;
    } while (fabs(result - xold) >= tolerance);

    // output the result
    printf("The inverse error function of %lf is %lf\n", x, result);
    printf("Number of iterations: %d\n", iteration);

    return 0;
}

In the terminal it should look something like this:

 Calculator for inverse error function.
 Enter the value for x: 0.5
 The inverse error function of 0.500000 is 0.476936
 Number of iterations: 5

Answer 3

At this time, the ISO C standard math library does not include erfinv(), or its single-precision variant erfinvf(). However, it is not too difficult to create one's own version, which I demonstrate below with an implementation of erfinvf() of reasonable accuracy and performance.

Looking at the graph of the inverse error function we observe that it is highly non-linear and is therefore difficult to approximate with a polynomial. One strategy do deal with this scenario is to "linearize" such a function by compositing it from simpler elementary functions (which can themselves be computed at high performance and excellent accuracy) and a fairly linear function which is more easily amenable to polynomial approximations or rational approximations of low degree.

Here are some approaches to erfinv linearization known from the literature, all of which are based on logarithms. Typically, authors differentiate between a main, fairly linear portion of the inverse error function from zero to a switchover point very roughly around 0.9 and a tail portion from the switchover point to unity. In the following, log() denotes the natural logarithm, R() denotes a rational approximation, and P() denotes a polynomial approximation.

A. J. Strecok, "On the Calculation of the Inverse of the Error Function." Mathematics of Computation, Vol. 22, No. 101 (Jan. 1968), pp. 144-158 (online)

β(x) = (-log(1-x²]))^½; erfinv(x) = x · R(x²) [main]; R(x) · β(x) [tail]

J. M. Blair, C. A. Edwards, J. H. Johnson, "Rational Chebyshev Approximations for the Inverse of the Error Function." Mathematics of Computation, Vol. 30, No. 136 (Oct. 1976), pp. 827-830 (online)

ξ = (-log(1-x))^-½; erfinv(x) = x · R(x²) [main]; ξ^-1 · R(ξ) [tail]

M. Giles, "Approximating the erfinv function." In GPU Computing Gems Jade Edition, pp. 109-116. 2011. (online)

w = -log(1-x²); s = √w; erfinv(x) = x · P(w) [main]; x · P(s) [tail]

The solution below generally follows the approach by Giles, but simplifies it in not requiring the square root for the tail portion, i.e. it uses two approximations of the type x · P(w). The code takes maximum advantage of the fused multiply-add operation FMA, which is exposed via the standard math functions fma() and fmaf() in C. Many common compute platforms, such as IBM Power, Arm64, x86-64, and GPUs offer this operation in hardware. Where no hardware support exists, the use of fma{f}() will likely make the code below unacceptably slow as the operation needs to be emulated by the standard math library. Also, functionally incorrect emulations of FMA are known to exist.

The accuracy of the standard math library's logarithm function logf() will have some impact on the accuracy of my_erfinvf() below. As long as the library provides a faithfully-rounded implementation with error < 1 ulp, the stated error bound should hold and it did for the few libraries I tried. For improved reproducability, I have included my own portable faithfully-rounded implementation, my_logf().

#include <math.h>
float my_logf (float);

/* compute inverse error functions with maximum error of 2.35793 ulp */
float my_erfinvf (float a)
{
    float p, r, t;
    t = fmaf (a, 0.0f - a, 1.0f);
    t = my_logf (t);
    if (fabsf(t) > 6.125f) { // maximum ulp error = 2.35793
        p =              3.03697567e-10f; //  0x1.4deb44p-32 
        p = fmaf (p, t,  2.93243101e-8f); //  0x1.f7c9aep-26 
        p = fmaf (p, t,  1.22150334e-6f); //  0x1.47e512p-20 
        p = fmaf (p, t,  2.84108955e-5f); //  0x1.dca7dep-16 
        p = fmaf (p, t,  3.93552968e-4f); //  0x1.9cab92p-12 
        p = fmaf (p, t,  3.02698812e-3f); //  0x1.8cc0dep-9 
        p = fmaf (p, t,  4.83185798e-3f); //  0x1.3ca920p-8 
        p = fmaf (p, t, -2.64646143e-1f); // -0x1.0eff66p-2 
        p = fmaf (p, t,  8.40016484e-1f); //  0x1.ae16a4p-1 
    } else { // maximum ulp error = 2.35002
        p =              5.43877832e-9f;  //  0x1.75c000p-28 
        p = fmaf (p, t,  1.43285448e-7f); //  0x1.33b402p-23 
        p = fmaf (p, t,  1.22774793e-6f); //  0x1.499232p-20 
        p = fmaf (p, t,  1.12963626e-7f); //  0x1.e52cd2p-24 
        p = fmaf (p, t, -5.61530760e-5f); // -0x1.d70bd0p-15 
        p = fmaf (p, t, -1.47697632e-4f); // -0x1.35be90p-13 
        p = fmaf (p, t,  2.31468678e-3f); //  0x1.2f6400p-9 
        p = fmaf (p, t,  1.15392581e-2f); //  0x1.7a1e50p-7 
        p = fmaf (p, t, -2.32015476e-1f); // -0x1.db2aeep-3 
        p = fmaf (p, t,  8.86226892e-1f); //  0x1.c5bf88p-1 
    }
    r = a * p;
    return r;
}

/* compute natural logarithm with a maximum error of 0.85089 ulp */
float my_logf (float a)
{
    float i, m, r, s, t;
    int e;

    m = frexpf (a, &e);
    if (m < 0.666666667f) { // 0x1.555556p-1
        m = m + m;
        e = e - 1;
    }
    i = (float)e;
    /* m in [2/3, 4/3] */
    m = m - 1.0f;
    s = m * m;
    /* Compute log1p(m) for m in [-1/3, 1/3] */
    r =             -0.130310059f;  // -0x1.0ae000p-3
    t =              0.140869141f;  //  0x1.208000p-3
    r = fmaf (r, s, -0.121484190f); // -0x1.f19968p-4
    t = fmaf (t, s,  0.139814854f); //  0x1.1e5740p-3
    r = fmaf (r, s, -0.166846052f); // -0x1.55b362p-3
    t = fmaf (t, s,  0.200120345f); //  0x1.99d8b2p-3
    r = fmaf (r, s, -0.249996200f); // -0x1.fffe02p-3
    r = fmaf (t, m, r);
    r = fmaf (r, m,  0.333331972f); //  0x1.5554fap-2
    r = fmaf (r, m, -0.500000000f); // -0x1.000000p-1
    r = fmaf (r, s, m);
    r = fmaf (i,  0.693147182f, r); //  0x1.62e430p-1 // log(2)
    if (!((a > 0.0f) && (a <= 3.40282346e+38f))) { // 0x1.fffffep+127
        r = a + a;  // silence NaNs if necessary
        if (a  < 0.0f) r = ( 0.0f / 0.0f); //  NaN
        if (a == 0.0f) r = (-1.0f / 0.0f); // -Inf
    }
    return r;
}

Answer 4

Quick & dirty, tolerance under +-6e-3. Work based on "A handy approximation for the error function and its inverse" by Sergei Winitzki.

C/C++ CODE:

float myErfInv2(float x){
   float tt1, tt2, lnx, sgn;
   sgn = (x < 0) ? -1.0f : 1.0f;

   x = (1 - x)*(1 + x);        // x = 1 - x*x;
   lnx = logf(x);

   tt1 = 2/(PI*0.147) + 0.5f * lnx;
   tt2 = 1/(0.147) * lnx;

   return(sgn*sqrtf(-tt1 + sqrtf(tt1*tt1 - tt2)));
}

MATLAB sanity check:

clear all,  close all, clc 

x = linspace(-1, 1,10000);
% x = 1 - logspace(-8,-15,1000);

a = 0.15449436008930206298828125;
% a = 0.147;

u = log(1-x.^2);  
u1 = 2/(pi*a) + u/2;    u2 = u/a;
y = sign(x).*sqrt(-u1+sqrt(u1.^2 - u2)); 

f = erfinv(x); axis equal

figure(1);
plot(x, [y; f]); legend('Approx. erf(x)', 'erf(x)')

figure(2);
e = f-y;
plot(x, e);

MATLAB Plots:

Answer 5

Also quick and dirty: if less precision is allowed than I share my own approximation with the inverse hyperbolic tangent - the parameters are sought by monte carle simulation where all random values are between the range of 0.5 and 1.5:

p1 = 1.4872301551536515
p2 = 0.5739159012216655
p3 = 0.5803635928651558

atanh( p^( 1 / p3 ) ) / p2 )^( 1 / p1 )

This comes from the algebraic reordering of my erf function approximation with the hyperbolic tangent, where the RMSE error is 0.000367354 for x between 1 and 4:

tanh( x^p1 * p2 )^p3