/* See {erfinv.h } */
/* Last edited on 2023-12-15 04:12:04 by stolfi */

#define _GNU_SOURCE
#include <math.h>

#include <affirm.h>
#include <erfinv.h>

double erfinv(double y)
  {
    /* Based on a Pascal version posted by David Heffernan to StackOverflow on 2011-05-12. */
    
    double a0 =  0.886226899, a1 = -1.645349621, a2 =  0.914624893, a3 = -0.140543331;
    double b0 = -2.118377725, b1 =  1.442710462, b2 = -0.329097515, b3 =  0.012229801;
    double c0 = -1.970840454, c1 = -1.624906493, c2 =  3.429567803, c3 =  1.641345311;
    double d0 =  3.543889200, d1 =  1.637067800;

    double y0 = 0.7;

    double x, z;

    demand((y >= -1.0) && (y <= 1.0), "{erfinv}: argument out of range");
    
    if (fabs(y) == 1) { return y*INF; }
    
    /* Rational function approximation: */
    if (y < -y0)
      { z = sqrt(-log((1.0+y)/2.0));
        x = -(((c3*z+c2)*z+c1)*z+c0)/((d1*z+d0)*z+1.0);
      }
    else if (y > y0)
      { z = sqrt(-log((1.0-y)/2.0));
        x = (((c3*z+c2)*z+c1)*z+c0)/((d1*z+d0)*z+1.0);
      }
    else
      { z = y*y;
        x = y*(((a3*z+a2)*z+a1)*z+a0)/((((b3*z+b3)*z+b1)*z+b0)*z+1.0);
      }
      
    /* Polish x to full accuracy with Newton iterations: */
    x = x - (erf(x) - y) / (2.0/sqrt(M_PI) * exp(-x*x));
    x = x - (erf(x) - y) / (2.0/sqrt(M_PI) * exp(-x*x));

    return x;
  }