/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:55 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "derf.h" #include #include /* PARAMETER translations */ #define FAC 1.12837916709551257389615890312154516e0 /* end of PARAMETER translations */ double /*FUNCTION*/ derf( double x) { double derf_v, u, y; static double sqeps, xbig; static LOGICAL32 first = TRUE; /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1999-12-27 DERF Krogh Added SAVE XMAX in DERFCE. *>> 1997-06-05 DERF Snyder Added SAVE EPS3 in DERFC *>> 1996-04-27 DERF Krogh Changes to use .C. and C%%. *>> 1996-03-30 DERF Krogh Added external statements. *>> 1995-11-28 DERF Krogh Converted SFTRAN to Fortran *>> 1995-10-24 DERF Krogh Removed blanks in numbers for C conversion. *>> 1994-10-19 DERF Krogh Changes to use M77CON *>> 1991-01-14 DERF CLL Changed to use generic names: abs and sign. *>> 1990-11-29 CLL *>> 1989-08-16 DERF CLL Changes for Cray *>> 1989-06-30 DERF CLL More digits in sqrt(PI), S.P. & D.P. more alike. *>> 1986-03-18 DERF Lawson Initial code. *--D replaces "?": ?ERF, ?ERFC, ?ERFCE, ?INITS, ?CSEVL, ?ERM1, ?ERF1, *--& ?ERFC1, ?ERFE1, ?ERFE2, ?ERFSP * * JULY 1977 EDITION. W.FULLERTON, C3, * LOS ALAMOS SCIENTIFIC LAB. * * REORGANIZATION OF FULLERTON'S DERF & DERFC * C.L.LAWSON & S.CHAN, JUNE 2 ,1983, JPL. * When DERFC is called with 0.5 .lt. X .le. 1.0 the original code * computed erfc from erfce, computing erfce from Cody's rational * formula whose coeffs (PS() and QS()) are only given to 18 decimal * digits. For Cray D.P. of 29D precision we can do better by * computing DERFC as 1 - erf in this interval, since the coeffs for * DERF are good up to 30D. * With this change DERFC will also be good up to 30 or 29D for all * X up to XMAX, where it would underflow. * We compute DERFCE only for X .ge. 0. DERFCE is computed from * DERF in [0, 1] and from DERFC in [1, XMAX] and so is good to * about 30D there. For X > XMAX we use the Cody rational * approximations for DERFCE which have coeffs good only to 18D. * Added variable EPS3. * ------------------------------------------------------------------ *-- Begin mask code changes * MACHINE DEPENDENT VALUES ARE SET ON THE FIRST ENTRY * TO THIS CODE. EXAMPLES OF THESE VALUES FOLLOW: * * SYSTEM NTERF NTERFC NTERC2 SQEPS XBIG XMAX * ------ ----- ------ ------ ----- ---- ---- * IEEE S.P. 7 9 10 .34E-3 4.01 9.18 * IEEE D.P. 12 25 24 .15E-7 6.01 26.53 * UNIVAC S.P. 8 11 11 .12E-3 4.26 9.29 * UNIVAC D.P. 14 30 27 .13E-8 6.40 26.58 * Cray S.P. 11 21 21 .12E-6 5.66 75.30 * Cray D.P. 19 55 45 .10E-13 8.04 74.89 * * The above values are derived from the following values: * * SYSTEM R1MACH(1) R1MACH(3) D1MACH(1) D1MACH(3) * ------ --------- --------- --------- --------- * IEEE 2**(-126) 2**(-24) 2**(-1022) 2**(-53) * VAX 2**(-128) 2**(-27) 2**(-128) 2**(-56) * UNIVAC 2**(-129) 2**(-27) 2**(-1025) 2**(-60) * Cray 2**(-8190) 2**(-47) 2**(-8100) 2**(-94) * *-- End mask code changes * ------------------------------------------------------------------ */ /* FAC = 2/sqrt(PI) */ /* ------------------------------------------------------------------ * */ if (first) derfsp( &first, &u, &sqeps, &xbig, &y ); y = fabs( x ); if (y <= sqeps) { u = FAC*x; } else if (y <= 1.0e0) { u = derf1( x ); } else if (y < xbig) { u = sign( 0.5e0 + (0.5e0 - derfc1( y )*exp( -y*y )/y), x ); } else { u = sign( 1.0e0, x ); } derf_v = u; return( derf_v ); } /* end of function */ /* ------------------------------------------------------------------ */ double /*FUNCTION*/ derfc( double x) { double derfc_v, u, y; static double eps3, sqeps, xbig, xmax; static LOGICAL32 first = TRUE; /* FAC = 2/sqrt(PI) */ if (first) derfsp( &first, &eps3, &sqeps, &xbig, &xmax ); if (x <= -xbig) { u = 2.0e0; } else { y = fabs( x ); if (y <= sqeps) { u = 0.5e0 + (0.5e0 - FAC*x); } else if (y <= 1.0e0) { if (x < 0.5e0) { /* *** Here X is in [-1,-SQEPS] or [SQEPS,0.5] * */ u = 0.5e0 + (0.5e0 - derf1( x )); } else { /* . *** Here X is in [0.5,1.0] * . Following test on 1.0D-19 added 8/89 to make better use * . of 28 D precision available on a Cray. The stored * . approximation for erfce is only good to about 18 D. * . On machines having precision greater than 18 D we get * . a more accurate value for erfc by computing U = erf * . and then 0.5D0 - (0.5D0 - U). * */ if (eps3 > 1.0e-19) { u = exp( -SQ(x) )*derfe1( y ); } else { u = 0.5e0 + (0.5e0 - derf1( x )); } } } else if (y <= xmax) { u = derfc1( y )*exp( -y*y )/y; if (x < 0.0e0) u = 2.0e0 - u; } else { u = 0.0e0; derm1( "DERFC", 1, 0, "X SO BIG DERFC UNDERFLOWS", "X" , x, '.' ); } } derfc_v = u; return( derfc_v ); } /* end of function */ /* ------------------------------------------------------------------ */ /* PARAMETER translations */ #define CUT2 4.0e0 #define CUT3 1.0e16 #define PII 0.56418958354775628694807945156077258e0 /* end of PARAMETER translations */ double /*FUNCTION*/ derfce( double x) { double derfce_v, u, y; static double xmax; static LOGICAL32 first = TRUE; /* PII = 1/sqrt(PI) */ /* CUT2 = 4.0 * CUT3 = 10.**16 WILL BE SATISFACTORY FOR ALL * MACHINES HAVING RELATIVE PRECISION LESS * THAN 32 DECIMAL PLACES AND OVERFLOW LIMIT * GREATER THAN 10.**32 */ if (first) derfsp( &first, &u, &y, &y, &xmax ); y = x; if (x < 0.0e0) { derm1( "DERFCE", 1, 0, "X MUST BE .GE. ZERO", "X", x, '.' ); u = 0.0e0; } else if (x < 1.0e0) { u = exp( y*y )*(0.5e0 + (0.5e0 - derf1( x ))); } else if (x < xmax) { u = derfc1( y )/x; } else if (x < CUT2) { u = derfe1( y ); } else if (x < CUT3) { u = derfe2( y ); } else { u = PII/x; } derfce_v = u; return( derfce_v ); } /* end of function */ /* PARAMETER translations */ #define SQRTPI 1.77245385090551602729816748334114518e0 /* end of PARAMETER translations */ void /*FUNCTION*/ derfsp( LOGICAL32 *first, double *eps3, double *sqeps, double *xbig, double *xmax) { long int _l0; static double eps3i, sqepsi, xbigi, xmaxi; static LOGICAL32 firsti = TRUE; /* PROCEDURE (SET PARAMETERS) */ if (firsti) { firsti = FALSE; eps3i = DBL_EPSILON/FLT_RADIX; sqepsi = sqrt( 2.0e0*eps3i ); xbigi = sqrt( -log( SQRTPI*eps3i ) ); xmaxi = sqrt( -log( SQRTPI*DBL_MIN ) ); xmaxi += -0.5e0*log( xmaxi )/xmaxi - 0.01e0; } *eps3 = eps3i; *sqeps = sqepsi; *xbig = xbigi; *xmax = xmaxi; /* PRINT*,'SQEPS = ',SQEPS, SQEPSI * PRINT*,'XBIG = ',XBIG, XBIGI * PRINT*,'XMAX = ',XMAX, XMAXI */ *first = FALSE; return; } /* end of function */ double /*FUNCTION*/ derf1( double x) { long int _l0; double derf1_v; static double erfcs[21]={-.49046121234691808039984544033376e-1, -.14226120510371364237824741899631e0,.10035582187599795575754676712933e-1, -.57687646997674847650827025509167e-3,.27419931252196061034422160791471e-4, -.11043175507344507604135381295905e-5,.38488755420345036949961311498174e-7, -.11808582533875466969631751801581e-8,.32334215826050909646402930953354e-10, -.79910159470045487581607374708595e-12,.17990725113961455611967245486634e-13, -.37186354878186926382316828209493e-15,.71035990037142529711689908394666e-17, -.12612455119155225832495424853333e-18,.20916406941769294369170500266666e-20, -.32539731029314072982364160000000e-22,.47668672097976748332373333333333e-24, -.65980120782851343155199999999999e-26,.86550114699637626197333333333333e-28, -.10788925177498064213333333333333e-29,.12811883993017002666666666666666e-31}; static long nterf = 0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Erfcs = &erfcs[0] - 1; /* end of OFFSET VECTORS */ /* PROCEDURE (DERF(X) FOR ABS(X) .LE. 1) */ /* SERIES FOR ERF ON THE INTERVAL 0. TO 1.00000D+00 * WITH WEIGHTED ERROR 1.28D-32 * LOG WEIGHTED ERROR 31.89 * SIGNIFICANT FIGURES REQUIRED 31.05 * DECIMAL PLACES REQUIRED 32.55 * *++ Save data by elements if ~.C. */ if (nterf == 0) dinits( erfcs, 21, 0.1e0*DBL_EPSILON/FLT_RADIX, &nterf ); derf1_v = x + x*dcsevl( 2.0e0*x*x - 1.0e0, erfcs, nterf ); return( derf1_v ); } /* end of function */ double /*FUNCTION*/ derfc1( double y) { long int _l0; static long int nterfc; double derfc1_v, eta, ysq; static double erc2cs[49]={-.6960134660230950112739150826197e-1, -.4110133936262089348982212084666e-1,.3914495866689626881561143705244e-2, -.4906395650548979161280935450774e-3,.7157479001377036380760894141825e-4, -.1153071634131232833808232847912e-4,.1994670590201997635052314867709e-5, -.3642666471599222873936118430711e-6,.6944372610005012589931277214633e-7, -.1371220902104366019534605141210e-7,.2788389661007137131963860348087e-8, -.5814164724331161551864791050316e-9,.1238920491752753181180168817950e-9, -.2690639145306743432390424937889e-10,.5942614350847910982444709683840e-11, -.1332386735758119579287754420570e-11,.3028046806177132017173697243304e-12, -.6966648814941032588795867588954e-13,.1620854541053922969812893227628e-13, -.3809934465250491999876913057729e-14,.9040487815978831149368971012975e-15, -.2164006195089607347809812047003e-15,.5222102233995854984607980244172e-16, -.1269729602364555336372415527780e-16,.3109145504276197583836227412951e-17, -.7663762920320385524009566714811e-18,.1900819251362745202536929733290e-18, -.4742207279069039545225655999965e-19,.1189649200076528382880683078451e-19, -.3000035590325780256845271313066e-20,.7602993453043246173019385277098e-21, -.1935909447606872881569811049130e-21,.4951399124773337881000042386773e-22, -.1271807481336371879608621989888e-22,.3280049600469513043315841652053e-23, -.8492320176822896568924792422399e-24,.2206917892807560223519879987199e-24, -.5755617245696528498312819507199e-25,.1506191533639234250354144051199e-25, -.3954502959018796953104285695999e-26,.1041529704151500979984645051733e-26, -.2751487795278765079450178901333e-27,.7290058205497557408997703680000e-28, -.1936939645915947804077501098666e-28,.5160357112051487298370054826666e-29, -.1378419322193094099389644800000e-29,.3691326793107069042251093333333e-30, -.9909389590624365420653226666666e-31,.2666491705195388413323946666666e-31}; static double erfccs[59]={.715179310202924774503697709496e-1,-.265324343376067157558893386681e-1, .171115397792085588332699194606e-2,-.163751663458517884163746404749e-3, .198712935005520364995974806758e-4,-.284371241276655508750175183152e-5, .460616130896313036969379968464e-6,-.822775302587920842057766536366e-7, .159214187277090112989358340826e-7,-.329507136225284321486631665072e-8, .722343976040055546581261153890e-9,-.166485581339872959344695966886e-9, .401039258823766482077671768814e-10,-.100481621442573113272170176283e-10, .260827591330033380859341009439e-11,-.699111056040402486557697812476e-12, .192949233326170708624205749803e-12,-.547013118875433106490125085271e-13, .158966330976269744839084032762e-13,-.472689398019755483920369584290e-14, .143587337678498478672873997840e-14,-.444951056181735839417250062829e-15, .140481088476823343737305537466e-15,-.451381838776421089625963281623e-16, .147452154104513307787018713262e-16,-.489262140694577615436841552532e-17, .164761214141064673895301522827e-17,-.562681717632940809299928521323e-18, .194744338223207851429197867821e-18,-.682630564294842072956664144723e-19, .242198888729864924018301125438e-19,-.869341413350307042563800861857e-20, .315518034622808557122363401262e-20,-.115737232404960874261239486742e-20, .428894716160565394623737097442e-21,-.160503074205761685005737770964e-21, .606329875745380264495069923027e-22,-.231140425169795849098840801367e-22, .888877854066188552554702955697e-23,-.344726057665137652230718495566e-23, .134786546020696506827582774181e-23,-.531179407112502173645873201807e-24, .210934105861978316828954734537e-24,-.843836558792378911598133256738e-25, .339998252494520890627359576337e-25,-.137945238807324209002238377110e-25, .563449031183325261513392634811e-26,-.231649043447706544823427752700e-26, .958446284460181015263158381226e-27,-.399072288033010972624224850193e-27, .167212922594447736017228709669e-27,-.704599152276601385638803782587e-28, .297976840286420635412357989444e-28,-.126252246646061929722422632994e-28, .539543870454248793985299653154e-29,-.238099288253145918675346190062e-29, .109905283010276157359726683750e-29,-.486771374164496572732518677435e-30, .152587726411035756763200828211e-30}; static long nterc2 = 0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Erc2cs = &erc2cs[0] - 1; double *const Erfccs = &erfccs[0] - 1; /* end of OFFSET VECTORS */ /* *(exp(-Y*Y) / Y) = PROCEDURE (DERFC(Y) FOR 1 .LT. Y .LT. XMAX) * *(1 / X) = PROCEDURE (DERFCE(X) FOR 1 .LT. X .LT. XMAX) */ /* SERIES FOR ERC2 ON THE INTERVAL 2.50000D-01 TO 1.00000D+00 * WITH WEIGHTED ERROR 2.67D-32 * LOG WEIGHTED ERROR 31.57 * SIGNIFICANT FIGURES REQUIRED 30.31 * DECIMAL PLACES REQUIRED 32.42 * *++ Save data by elements if ~.C. */ /* SERIES FOR ERFC ON THE INTERVAL 0. TO 2.50000D-01 * WITH WEIGHTED ERROR 1.53D-31 * LOG WEIGHTED ERROR 30.82 * SIGNIFICANT FIGURES REQUIRED 29.47 * DECIMAL PLACES REQUIRED 31.70 * *++ Save data by elements if ~.C. */ if (nterc2 == 0) { eta = 0.1e0*DBL_EPSILON/FLT_RADIX; dinits( erc2cs, 49, eta, &nterc2 ); dinits( erfccs, 59, eta, &nterfc ); } ysq = y*y; if (ysq <= 4.0e0) { derfc1_v = 0.5e0 + dcsevl( (8.0e0/ysq - 5.0e0)/3.0e0, erc2cs, nterc2 ); } else { derfc1_v = 0.5e0 + dcsevl( 8.0e0/ysq - 1.0e0, erfccs, nterfc ); } return( derfc1_v ); } /* end of function */ double /*FUNCTION*/ derfe1( double y) { double derfe1_v, psum, qsum; static double ps[9]={1.00000000000036828e0,1.87051017604560834e0, 1.74642369370058320e0,1.02438464807598001e0,4.07413180167223764e-1, 1.11870870991098165e-1,2.07045775788719818e-2,2.37133372752999036e-3, 1.29992515945788642e-4}; static double qs[10]={1.00000000000000000e0,2.99888934314798253e0, 4.13030795287321183e0,3.43830153103630866e0,1.91273588328781533e0, 7.40352738163508723e-1,2.00387662412610424e-1,3.68131014202168126e-2, 4.20307996290648223e-3,2.30405728794132537e-4}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Ps = &ps[0] - 1; double *const Qs = &qs[0] - 1; /* end of OFFSET VECTORS */ /* PROCEDURE (DERFCE(Y) FOR 15/32 .LE. Y .LE. 4) */ psum = Ps[1] + y*(Ps[2] + y*(Ps[3] + y*(Ps[4] + y*(Ps[5] + y*(Ps[6] + y*(Ps[7] + y*(Ps[8] + y*Ps[9]))))))); qsum = Qs[1] + y*(Qs[2] + y*(Qs[3] + y*(Qs[4] + y*(Qs[5] + y*(Qs[6] + y*(Qs[7] + y*(Qs[8] + y*(Qs[9] + y*Qs[10])))))))); derfe1_v = psum/qsum; return( derfe1_v ); } /* end of function */ double /*FUNCTION*/ derfe2( double y) { double derfe2_v, psum, qsum, ysq; static double pl[6]={-2.82094791773876911e-1,-6.88752606824337911e0, -5.38632485753818598e1,-1.54309751654255924e2,-1.30749393763753716e2, -6.98670450907125305e0}; static double ql[6]={1.00000000000000000e0,2.59156442057350244e1, 2.26063711030173368e2,8.02050727433854173e2,1.09991209268230208e3, 4.28227932949102556e2}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Pl = &pl[0] - 1; double *const Ql = &ql[0] - 1; /* end of OFFSET VECTORS */ /* PROCEDURE (DERFCE(Y) FOR 4 .LE. Y .LE. 10**16) */ /* PII = 1/sqrt(PI) */ ysq = 1.e0/y/y; psum = Pl[1] + ysq*(Pl[2] + ysq*(Pl[3] + ysq*(Pl[4] + ysq*(Pl[5] + ysq*Pl[6])))); qsum = Ql[1] + ysq*(Ql[2] + ysq*(Ql[3] + ysq*(Ql[4] + ysq*(Ql[5] + ysq*Ql[6])))); derfe2_v = (PII + ysq*psum/qsum)/y; return( derfe2_v ); } /* end of function */