Category : C Source Code
Archive   : JPL-C2.ZIP
Filename : INVERF.C

/* 1.1 10-08-85 (inverf.c)
************************************************************************
* Robert C. Tausworthe *
* Jet Propulsion Laboratory *
* Pasadena, CA 91009 1985 *
************************************************************************/

#include "defs.h"
#include "stdtyp.h"
#include "errno.h"
#include "mathcons.h"

/*----------------------------------------------------------------------*

This routine uses the approximations given in J. M. Blair,
et. al., "Rational Chebyshev Approximations for the Inverse Error
Function," Math Comp. Vol. 30, No. 136, Oct. 1976, pp. 827-830, with
the appendices given in separate microfiche.

*----------------------------------------------------------------------*/

LOCAL int dP[4] = {7, 8, 11, 8}; /* degrees of P(x) */
LOCAL int dQ[4] = {7, 8, 10, 9}; /* degrees of Q(x) */

LOCAL double P[4][12] = /* numerator polynomials */
{
/* P[0][0] */ -.30866886527764497339e2, /* Table 19 of Ref. below */
/* P[0][1] */ 0.20652786402942339589e3, /* |x| <= 0.75 */
/* P[0][2] */ -.52856770835093823310e3,
/* P[0][3] */ 0.64880509544036249088e3,
/* P[0][4] */ -.39205509901971391289e3,
/* P[0][5] */ 0.10706278097770074402e3,
/* P[0][6] */ -.10303488455439678272e2,
/* P[0][7] */ 0.15641510821923860000e0,
/* P[0][8] */ 0.0,
/* P[0][9] */ 0.0,
/* P[0][10] */ 0.0,
/* P[0][11] */ 0.0,

/* P[1][0] */ 0.94897362808681080020e-2, /* Table 39 of Ref. below */
/* P[1][1] */ -.24758242362823355485e0, /* 0.75 <= x <= 0.9375 */
/* P[1][2] */ 0.25349389220714893916e1,
/* P[1][3] */ -.12954198980646771502e2,
/* P[1][4] */ 0.34810057749357500873e2,
/* P[1][5] */ -.47644367129787181802e2,
/* P[1][6] */ 0.29631331505876308122e2,
/* P[1][7] */ -.64200071507209448654e1,
/* P[1][8] */ 0.21489185007307062000e0,
/* P[1][9] */ 0.0,
/* P[1][10] */ 0.0,
/* P[1][11] */ 0.0,

/* P[2][0] */ 0.19953288964537210824e-5, /* Table 60 of Ref. below */
/* P[2][1] */ 0.21273702631785953343e-3, /* 0.9375 <= x <= 1 - 10^(-100) */
/* P[2][2] */ 0.58975595952407247651e-2,
/* P[2][3] */ 0.59481901452735809123e-1,
/* P[2][4] */ 0.31328289030932667506e0,
/* P[2][5] */ 0.13630199956442260990e1,
/* P[2][6] */ 0.34152815205652930673e1,
/* P[2][7] */ 0.30184181468933606839e1,
/* P[2][8] */ 0.20842433546207339433e1,
/* P[2][9] */ 0.85545635026743499993e0,
/* P[2][10] */ 0.40273918408712893132e-2,
/* P[2][11] */ -.15196139115744716810e-3,

/* P[3][0] */ .45344689563209398449e-11, /* Table 82 of Ref. below */
/* P[3][1] */ .46156006321345332510e-8, /* 1 - 10^(-100) <= x */
/* P[3][2] */ .12964481560643197452e-5, /* <= 1 - 10^(-10000) */
/* P[3][3] */ .13714329569665128933e-3,
/* P[3][4] */ .60537914739162189689e-2,
/* P[3][5] */ .11279046353630280004e0,
/* P[3][6] */ .82810030904462690215e0,
/* P[3][7] */ .19507620287580568829e1,
/* P[3][8] */ .69952990607058154857e0,
/* P[3][9] */ 0.0,
/* P[3][10] */ 0.0,
/* P[3][11] */ 0.0
};

LOCAL double Q[4][11] = /* denominator polynomials */
{
/* Q[0][0] */ -.28460290173882339383e2, /* Table 19 */
/* Q[0][1] */ 0.20518924149238530630e3,
/* Q[0][2] */ -.57617125090127638064e3,
/* Q[0][3] */ 0.79669388170563770334e3,
/* Q[0][4] */ -.56509305564023424022e3,
/* Q[0][5] */ 0.19450320483954087700e3,
/* Q[0][6] */ -.27371852306002662877e2,
/* Q[0][7] */ 1.0,
/* Q[0][8] */ 0.0,
/* Q[0][9] */ 0.0,
/* Q[0][10] */ 0.0,

/* Q[1][0] */ 0.67544512778850945940e-2, /* Table 39 */
/* Q[1][1] */ -.18611650627372178511e0,
/* Q[1][2] */ 0.20369295047216351160e1,
/* Q[1][3] */ -.11315360624238054876e2,
/* Q[1][4] */ 0.33880176779595142684e2,
/* Q[1][5] */ -.53715373448862143348e2,
/* Q[1][6] */ 0.41409991778428888715e2,
/* Q[1][7] */ -.12831383833953226499e2,
/* Q[1][8] */ 1.0,
/* Q[1][9] */ 0.0,
/* Q[1][10] */ 0.0,

/* Q[2][0] */ .19953210379374212953e-5, /* Table 60 */
/* Q[2][1] */ .21274156963404084598e-3,
/* Q[2][2] */ .59037062023731354671e-2,
/* Q[2][3] */ .59959150110861092334e-1,
/* Q[2][4] */ .32318083080817836442e0,
/* Q[2][5] */ .14378337804749344527e1,
/* Q[2][6] */ .37644571508257969664e1,
/* Q[2][7] */ .44081435698143841173e1,
/* Q[2][8] */ .42508710497182804606e1,
/* Q[2][9] */ .22127469427969785343e1,
/* Q[2][10] */ 1.0,

/* Q[3][0] */ .45344687377088206782e-11, /* Table 82 */
/* Q[3][1] */ .46156017600933592558e-8,
/* Q[3][2] */ .12964671850944981712e-5,
/* Q[3][3] */ .13715891988350205065e-3,
/* Q[3][4] */ .60574830550097140404e-2,
/* Q[3][5] */ .11311889334355782064e0,
/* Q[3][6] */ .84001814918178042918e0,
/* Q[3][7] */ .21238242087454993541e1,
/* Q[3][8] */ .15771922386662040545e1,
/* Q[3][9] */ 1.0,
/* Q[3][10] */ 0.0
};

/************************************************************************/
double
inverf(x) /* inverse of erf(x), which see. */

/*----------------------------------------------------------------------*/
double x;
{
int interval, sign;
double y, inverfc(), ratfun();
LOCAL double bias[2] = {0.5625, 0.87890625};

sign = (x < 0.0 ? 1 : 0);
if (sign)
x = -x;
if (x >= 1.0)
{ errno = EDOM;
return (sign ? -INFINITY : INFINITY);
}
if (x > 0.9375)
{ y = inverfc(1.0 - x);
return (sign ? -y : y);
}
if (x < 0.75)
interval = 0;
else interval = 1;
y = x * ratfun(x*x - bias[interval], P[interval], Q[interval],
dP[interval], dQ[interval]);
return (sign ? -y : y);
}

#define BOUNDARY 0.065901028 /* 1 / sqrt(100 * log(10)) */

/*\p*********************************************************************/
double
inverfc(x) /* returns double inverse of erfc(y) */

/*----------------------------------------------------------------------*/
double x;
{
int sign, interval;
double y, log(), fabs(), inverf(), sqrt(), ratfun();

y = fabs(1.0 - x);
if (y >= 1.0)
{ errno = EDOM;
return (x > 2 ? -INFINITY : INFINITY);
}
if (x > 1.0)
{ x = 2.0 - x;
sign = 1;
} else sign = 0;
if (x IS 0.0)
{ errno = ERANGE;
return (sign ? -INFINITY : INFINITY);
}
if (x > 0.0625)
{ y = inverf(1.0 - x);
return (sign ? -y : y);
}
y = 1.0 / sqrt(-log(x));
if (y > BOUNDARY)
interval = 2;
else interval = 3;
y = ratfun(y, P[interval], Q[interval], dP[interval], dQ[interval])/y;
return (sign ? -y : y);
}