msci Code

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include "peteys_tmpl_lib.h"
#include "supernewton.h"

#define ERFINV_XTOL 0.
#define ERFINV_YTOL 1e-15

#define ERFCINV_XTOL 0.
#define ERFCINV_YTOL 1e-50

#define ERFCINV_MAXITER 1000

//the inverse error function is normally not included with standard math
//libraries (it is absent from the GSL)
//this is my attempt to produce such a function
//while the current attempts do not work satisfactorily, they have
//revealed to me some limitations in the "supernewton" algorithm that
//need to be addressed
//a bit more polishing and it would probably work quite well...

void erfc_plus_deriv(double x, void *param, double *y, double *dydx) {
  double f;
  double y0;
  double dfdx;

  *y=gsl_sf_erfc(x)-*(double *) param; 
  *dydx=-M_2_SQRTPI*exp(-x*x);
}

void erf_plus_deriv(double x, void *param, double *y, double *dydx) {

  *y=gsl_sf_erf(x)-*(double *) param; 
  *dydx=M_2_SQRTPI*exp(-x*x);
}

double erfcinv2(double x) {
  double absx;
  double y;
  double y1, y2;
  long err;
  double sgnx;
  double xc;
  void *param;

  gsl_error_handler_t *old_handler;

  if (x==0.) return 0.;
  
  absx=fabs(x);
  sgnx=x/absx;

  old_handler=gsl_set_error_handler_off();
  if (absx < gsl_sf_erfc(25)) {
    return sgnx*INFINITY;
  } else if (absx < gsl_sf_erfc(2.5)) {
    param=&x;
    y1=sgnx*2.5;
    y2=sgnx*25;
    y=supernewton(&erfc_plus_deriv, param, y1, y2, ERFCINV_XTOL, ERFCINV_YTOL, ERFCINV_MAXITER, err);
  } else {
    xc=1-x;
    param=&xc;
    y1=0;
    y2=sgnx*2.5;
    y=supernewton(&erf_plus_deriv, param, y1, y2, ERFINV_XTOL, ERFINV_YTOL, ERFCINV_MAXITER, err);
  }

  gsl_set_error_handler(old_handler);

  printf("erfinv iterations=%ld\n", err);

  return y;

}

#define ERFCINV_NTABLE 20
double erfcinv_ytable[ERFCINV_NTABLE]={0.1000000000000000055511151,    0.1339997106006672877853703,     0.179559224410625883905368,    0.2406088410670414179381993,    0.3224151507094550339616035,    0.4320353688833750149811408,    0.5789261439962479771637049,    0.7757593575465746571495629,      1.03951529407000609062095,     1.392947485703483589958296,     1.866545599661939336399996,      2.50116570177648966932793,     3.351554802023644086261811,     4.491073735334452088352464,     6.018025808210742511050739,     8.064137166875866569171194,     10.80592046605450917695634,     14.47990215225132537568697,     19.40302697927656438992017,                             26};
double erfcinv_xtable[ERFCINV_NTABLE]={    0.8875370839817151580319887,    0.8496976572790083670483341,    0.7995457048771270613940487,    0.7336514856277412954810302,    0.6484159530328690301814731,     0.541206016476690532357452,     0.412943213103930006901976,    0.2726023114674728797801606,     0.141535585714363754128442,   0.04884694086378559702010804,  0.008298088947439162185726325, 0.0004044201743987975482610975, 2.139142245461042382694916e-06, 2.134509774098865909524717e-10, 1.727765284054945417466263e-17,  3.97453051075275160673825e-30, 1.009859790297311418956122e-52, 3.405357962804753120513636e-93, 9.139006961558066437665542e-166, 5.663192408856141018983002e-296};

double erfcinv(double x) {
  double absx;
  double y;
  double y1, y2;
  long err;
  long ind;
  double sgnx;
  double ytol;

  gsl_error_handler_t *old_handler;

  if (x==0.) return 0.;
  
  absx=fabs(x);
  sgnx=x/absx;

  void *param=&x;

  ind=bin_search(erfcinv_xtable, ERFCINV_NTABLE, absx);
  if (ind < 0) {
    y1=-erfcinv_ytable[0];
    y2=erfcinv_ytable[0];
    ytol=0;
  } else if (ind==ERFCINV_NTABLE) {
    return sgnx*INFINITY;
  } else {
    y1=erfcinv_ytable[ind]*sgnx;
    y2=erfcinv_ytable[ind+1]*sgnx;
    ytol=0;
  }

  //y1=-100.;
  //y2=100.;

  old_handler=gsl_set_error_handler_off();

  y=supernewton(&erfc_plus_deriv, param, y1, y2, ERFCINV_XTOL, ytol, ERFCINV_MAXITER, err);

  gsl_set_error_handler(old_handler);

  printf("erfinv iterations=%ld\n", err);
  if (err < 0) {
    return sgnx*INFINITY;
  }

  return y;

}

#define ERFINV_NTABLE 8
double erfinv_ytable[ERFINV_NTABLE+1]={0.1000000000000000055511151,     2.704865449790269771312978,     4.159675982375273584068509,     4.972183855244742822776516,     5.425967374602302939479159,     5.679404281882109550849691,     5.820948130945836851424247,     5.900000000000000355271368};
double erfinv_xtable[ERFINV_NTABLE]={0.1124629160182848974791625,    0.9998693644789686807428097,    0.9999999959630012646982777,    0.9999999999979600762145537,    0.9999999999999832356323282,    0.9999999999999990007992778,    0.9999999999999997779553951,    0.9999999999999998889776975};
double erfinv_tol[ERFINV_NTABLE]={1e-12, 1e-12, 1e-12, 1e-13, 1e-15, 1e-18,
		1e-20, 1e-22};

double erfinv(double x) {
  double absx;
  double y;
  double y1, y2;
  long err;
  long ind;
  double sgnx;
  double ytol;

  gsl_error_handler_t *old_handler;

  if (x==0.) return 0.;
  
  absx=fabs(x);
  sgnx=x/absx;

  void *param=&x;

  ind=bin_search(erfinv_xtable, ERFINV_NTABLE, absx);
  if (ind < 0) {
    y1=-erfinv_ytable[0];
    y2=erfinv_ytable[0];
    ytol=erfinv_tol[ind];
  } else if (ind==ERFINV_NTABLE) {
    return sgnx*INFINITY;
  } else {
    y1=erfinv_ytable[ind]*sgnx;
    y2=erfinv_ytable[ind+1]*sgnx;
    ytol=1e-12;
  }

  //y1=-100.;
  //y2=100.;

  old_handler=gsl_set_error_handler_off();

  y=supernewton(&erf_plus_deriv, param, y1, y2, ERFINV_XTOL, ytol, ERFCINV_MAXITER, err);

  gsl_set_error_handler(old_handler);

  printf("erfinv iterations=%ld\n", err);
  if (err < 0) {
    return sgnx*INFINITY;
  }

  return y;

}