/*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 #include "fcrt.h" #include "derfi.h" #include #include 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 */