/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:41 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "derfi.h"
#include <float.h>
#include <stdlib.h>
double /*FUNCTION*/ derfi(
double x)
{
long int j;
double arg, derfi_v, fsign, s;
static double d[6]={-1.548813042373261659512742e0,2.565490123147816151928163e0,
-.5594576313298323225436913e0,2.287915716263357638965891e0,-9.199992358830151031278420e0,
2.794990820124599493768426e0};
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const D = &d[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*--D replaces "?": ?ERFI, ?ERFCI, ?ERFIX, ?ERM1
*>> 1998-11-01 DERFI Krogh Removed some equivalence for "mangle".
*>> 1996-06-18 DERFI Krogh Changes to use .C. and C%%. J not changed.
*>> 1996-03-30 DERFI Krogh Added external statements.
*>> 1995-11-28 DERFI Krogh Removed multiple entries.
*>> 1995-11-03 DERFI Krogh Removed blanks in numbers for C conversion.
*>> 1994-10-20 DERFI Krogh Changes to use M77CON
*>> 1994-04-20 DERFI CLL Edited type stmts to make DP & SP files similar
*>> 1987-10-29 DERFI Snyder Initial code.
*
* For -1.0 .LT. X .LT. 1.0 calculate the inverse of the error
* function. That is, X = ERF(ERFI).
*
* For 0.0 .LT. X .LT. 2.0 calculate the inverse of the
* complementary error function. that is, X = ERFC(ERFCI). This
* calculation is carried out by invoking the alternate entry *ERFCI.
*
* If X is out of range, program execution is terminated by calling
* the error message processor.
*
* This subprogram uses approximations due to A. Strecok from
* Mathematics of Computation 22, (1968) pp 144-158.
* */
/* ***** Parameters *****************************************
*
* MAX... is the position in C of the last coefficient of a Chebyshev
* polynomial expansion.
* MIN... is the position in C of the first coefficient of a Chebyshev
* polynomial expansion.
* NC is the upper dimension of the array of coefficients.
* NDELTA is the number of coefficients of the Chebyshev polynomial
* expansion used to approximate R(X) in the range
* 0.9975 .LT. X .LE. 1-5.0D-16
* NLAMDA is the number of coefficients of the Chebyshev polynomial
* expansion used to approximate R(X) in the range
* 0.8 .LT. X .LE. 0.9975.
* NMU is the number of coefficients of the Chebyshev polynomial
* expansion used to approximate R(X) in the range
* 5.0D-16 .GT. 1-X .GE. 1.D-300.
* NXI is the number of coefficients of the Chebyshev polynomial
* expansion used to approximate DERFCI(X)/X in the
* range 0.0 .LE. X .LE. 0.8.
*
*
* ***** External References ********************************
*
* D1MACH Provides the round-off level. Used to calculate the number
* of coefficients to retain in each Chebyshev expansion.
* DERM1 Prints an error message and stops if X .LE. -1.0 or
* X .GE. 1.0 (ERFI) or X .LE. 0.0 or X .GE. 2.0 (ERFCI).
* LOG Calculates the natural logarithm.
* SQRT Calculates the square root.
*
*
* ***** Local Variables ***********************************
*
* ARG If ERFI or ERFCI is being approximated by a Chebyshev
* expansion then ARG is the argument of ERFI or the argument
* that would be used if ERFCI(X) were computed as ERFC(1-X),
* that is, ARG = X if ERFI is being computed, or ARG = 1-X if
* ERFCI is being computed. If ERFI or ERFCI is being computed
* using the initial approximation ERFI=SQRT(-LOG((1-X)*(1+X))),
* then ARG is that initial approximation.
* C contains the coefficients of polynomial expansions. They are
* stored in C in the order DELTA(0..37), LAMDA(0..26),
* MU(0..25), XI(0..38).
* D are used to scale the argument of the Chebyshev polynomial
* expansion in the range 1.D-300 .LT. 1-X .LT. 0.2.
* DELTA are coefficients of the Chebyshev polynomial expansion of R(X)
* for 0.9975 .LT. X .LE. 1-5.0D-16.
* FIRST is a logical SAVE variable indicating whether it is necessary
* to calculate the number of coefficients to use for each
* Chebyshev expansion.
* FSIGN is X or 1.0 - X. It is used to remember the sign to be
* assigned to the function value.
* I, J are used as indices.
* IMIN is the minimum index of a coefficient in the Chebyshev
* polynomial expansion to be used.
* JIX is an array containing MINXI, MAXXI, MINLAM, MAXLAM, MINDEL,
* MAXDEL, MINMU, MAXMU in locations -1..6
* LAMDA are coefficients of the Chebyshev polynomial expansion of R(X)
* for 0.8 .LT. X .LE. 0.9975.
* MU are coefficients of the Chebyshev polynomial expansion of R(X)
* for 5.0D-16 .GT. 1-X .GE. 1.D-300.
* S2 is 2.0 * S.
* S is the argument of the Chebyshev polynomial expansion.
* W1..W3 are adjacent elements of the recurrence used to evaluate the
* Chebyshev polynomial expansion.
* XI are coefficients of the Chebyshev polynomial expansion of
* ERFC(X)/X for 0.0 .LE. X .LE. 0.8.
* */
fsign = x;
arg = fabs( x );
if (arg < 0.0e0 || arg >= 1.0e0)
{
derm1( "DERFI", 1, 2, "Argument out of range", "X", x, '.' );
/* In case the error level is shifted to zero by the caller: */
derfi_v = 0.0e0;
return( derfi_v );
}
if (arg == 0.0e0)
{
derfi_v = 0.0e0;
return( derfi_v );
}
if (arg <= 0.8e0)
{
s = 3.125e0*arg*arg - 1.0e0;
j = -1;
}
else
{
if (arg <= 0.9975e0)
{
j = 1;
}
else
{
j = 3;
}
arg = sqrt( -log( (1.0e0 - arg)*(1.0e0 + arg) ) );
s = D[j]*arg + D[j + 1];
}
derfi_v = sign( arg*derfix( s, j ), fsign );
return( derfi_v );
} /* end of function */
/* ***** entry ERFCI ****************************************
* */
double /*FUNCTION*/ derfci(
double x)
{
long int j;
double arg, derfci_v, fsign, s;
static double d[6]={-1.548813042373261659512742e0,2.565490123147816151928163e0,
-.5594576313298323225436913e0,2.287915716263357638965891e0,-9.199992358830151031278420e0,
2.794990820124599493768426e0};
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const D = &d[0] - 1;
/* end of OFFSET VECTORS */
/* Calculate the inverse of the complementary error function.
* */
/* Decide which approximation to use, and calculate the argument of
* the Chebyshev polynomial expansion.
* */
if (x <= 0.0e0 || x >= 2.0e0)
{
derm1( "DERFCI", 1, 2, "Argument out of range", "X", x, '.' );
/* In case the error level is shifted to zero by the caller: */
derfci_v = 0.0e0;
}
if (x == 1.0e0)
{
derfci_v = 0.0e0;
return( derfci_v );
}
fsign = 1.0e0 - x;
arg = fabs( fsign );
if (arg <= 0.8e0)
{
s = 3.125e0*arg*arg - 1.0e0;
j = -1;
}
else
{
arg = 2.0e0 - x;
if (x < 1.0e0)
{
s = x;
}
else
{
s = arg;
}
arg = sqrt( -log( x*arg ) );
if (s < 5.0e-16)
{
j = 5;
s = D[5]/sqrt( arg ) + D[6];
}
else
{
if (s >= 0.0025e0)
{
j = 1;
}
else if (s >= 5.0e-16)
{
j = 3;
}
s = D[j]*arg + D[j + 1];
}
}
derfci_v = sign( arg*derfix( s, j ), fsign );
return( derfci_v );
} /* end of function */
/* PARAMETER translations */
#define MAXDEL (MINDEL + NDELTA)
#define MAXLAM (MINLAM + NLAMDA)
#define MAXMU (MINMU + NMU)
#define MAXXI (MINXI + NXI)
#define MINDEL 0
#define MINLAM (MAXDEL + 1)
#define MINMU (MAXLAM + 1)
#define MINXI (MAXMU + 1)
#define NC MAXXI
#define NDELTA 37
#define NLAMDA 26
#define NMU 25
#define NXI 38
/* end of PARAMETER translations */
double /*FUNCTION*/ derfix(
double s,
long j)
{
long int _l0, i, imin;
double derfix_v, s2, w1, w2, w3;
static double c[NC-(0)+1]={.9566797090204925274526373e0,-.0231070043090649036999908e0,
-.0043742360975084077333218e0,-.0005765034226511854809364e0,-.0000109610223070923931242e0,
.0000251085470246442787982e0,.0000105623360679477511955e0,.0000027544123300306391503e0,
.0000004324844983283380689e0,-.0000000205303366552086916e0,-.0000000438915366654316784e0,
-.0000000176840095080881795e0,-.0000000039912890280463420e0,-.0000000001869324124559212e0,
.0000000002729227396746077e0,.0000000001328172131565497e0,.0000000000318342484482286e0,
.0000000000016700607751926e0,-.0000000000020364649611537e0,-.0000000000009648468127965e0,
-.0000000000002195672778128e0,-.0000000000000095689813014e0,.0000000000000137032572230e0,
.0000000000000062538505417e0,.0000000000000014584615266e0,.0000000000000001078123993e0,
-.0000000000000000709229988e0,-.0000000000000000391411775e0,-.0000000000000000111659209e0,
-.0000000000000000015770366e0,.0000000000000000002853149e0,.0000000000000000002716662e0,
.0000000000000000000957770e0,.0000000000000000000176835e0,-.0000000000000000000009828e0,
-.0000000000000000000020464e0,-.0000000000000000000008020e0,-.0000000000000000000001650e0,
.9121588034175537733059200e0,-.0162662818676636958546661e0,.0004335564729494453650589e0,
.0002144385700744592065205e0,.0000026257510757648130176e0,-.0000030210910501037969912e0,
-.0000000124060618367572157e0,.0000000624066092999917380e0,-.0000000005401247900957858e0,
-.0000000014232078975315910e0,.0000000000343840281955305e0,.0000000000335848703900138e0,
-.0000000000014584288516512e0,-.0000000000008102174258833e0,.0000000000000525324085874e0,
.0000000000000197115408612e0,-.0000000000000017494333828e0,-.0000000000000004800596619e0,
.0000000000000000557302987e0,.0000000000000000116326054e0,-.0000000000000000017262489e0,
-.0000000000000000002784973e0,.0000000000000000000524481e0,.0000000000000000000065270e0,
-.0000000000000000000015707e0,-.0000000000000000000001475e0,.0000000000000000000000450e0,
.9885750640661893136460358e0,.0108577051845994776160281e0,-.0017511651027627952494825e0,
.0000211969932065633437984e0,.0000156648714042435087911e0,-.0000005190416869103124261e0,
-.0000000371357897426717780e0,.0000000012174308662357429e0,-.0000000001768115526613442e0,
-.0000000000119372182556161e0,.0000000000003802505358299e0,-.0000000000000660188322362e0,
-.0000000000000087917055170e0,-.0000000000000003506869329e0,-.0000000000000000697221497e0,
-.0000000000000000109567941e0,-.0000000000000000011536390e0,-.0000000000000000001326235e0,
-.0000000000000000000263938e0,.0000000000000000000005341e0,-.0000000000000000000022610e0,
.0000000000000000000009552e0,-.0000000000000000000005250e0,.0000000000000000000002487e0,
-.0000000000000000000001134e0,.0000000000000000000000420e0,.9928853766189408231495800e0,
.1204675161431044864647846e0,.0160781993420999447257039e0,.0026867044371623158279591e0,
.0004996347302357262947170e0,.0000988982185991204409911e0,.0000203918127639944337340e0,
.0000043272716177354218758e0,.0000009380814128593406758e0,.0000002067347208683427411e0,
.0000000461596991054300078e0,.0000000104166797027146217e0,.0000000023715009995921222e0,
.0000000005439284068471390e0,.0000000001255489864097987e0,.0000000000291381803663201e0,
.0000000000067949421808797e0,.0000000000015912343331469e0,.0000000000003740250585245e0,
.0000000000000882087762421e0,.0000000000000208650897725e0,.0000000000000049488041039e0,
.0000000000000011766394740e0,.0000000000000002803855725e0,.0000000000000000669506638e0,
.0000000000000000160165495e0,.0000000000000000038382583e0,.0000000000000000009212851e0,
.0000000000000000002214615e0,.0000000000000000000533091e0,.0000000000000000000128488e0,
.0000000000000000000031006e0,.0000000000000000000007491e0,.0000000000000000000001812e0,
.0000000000000000000000439e0,.0000000000000000000000106e0,.0000000000000000000000026e0,
.0000000000000000000000006e0,.0000000000000000000000002e0};
static LOGICAL32 first = TRUE;
static long jix[6-(-1)+1]={MINXI,MAXXI,MINLAM,MAXLAM,MINDEL,MAXDEL,
MINMU,MAXMU};
/* Subroutine where most of calculations are done. */
/* ***** Equivalence Statements *****************************
*
* Equivalence statements connecting arrays DELTA, LAMDA, MU, XI removed
* by FTK to make "mangle" work. All references to these arrays have
* been replaced by references to C.
*
* ***** Data Statements ************************************
*
* DELTA(J), J = 0..NDELTA, then
* LAMDA(J), J = 0..NLAMDA, then
* MU(J), J = 0..NMU, then
* XI(J), J = 0..NXI
*
*++ With first index 0, save data by elements if ~.C. */
/* ***** Procedures *****************************************
*
* Decide which approximation to use, and calculate the argument of
* the Chebyshev polynomial expansion.
*
*
* If this is the first call, calculate the degree of each expansion.
* */
if (first)
{
first = FALSE;
s2 = 0.5e0*DBL_EPSILON/FLT_RADIX;
for (imin = -1; imin <= 5; imin += 2)
{
for (i = jix[imin-(-1)]; i <= jix[imin + 1-(-1)]; i++)
{
if (fabs( c[i] ) < s2)
{
jix[imin + 1-(-1)] = i;
goto L_120;
}
}
L_120:
;
}
}
/* Evaluate the Chebyshev polynomial expansion.
* */
s2 = s + s;
w1 = 0.0e0;
w2 = 0.0e0;
imin = jix[j-(-1)];
i = jix[j + 1-(-1)];
L_200:
w3 = w2;
w2 = w1;
w1 = (s2*w2 - w3) + c[i];
i -= 1;
if (i > imin)
goto L_200;
derfix_v = (s*w1 - w2) + c[imin];
return( derfix_v );
} /* end of function */