/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:03 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dfrenl.h" #include #include /* File: dfrenl.[for|f|c] * Contains procedures: DFRENC(), DFRENF(), DFRENG(), DFRENS() * and private low-level procedure: DFREN1(). * Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2000-12-01 DFRENL Krogh Removed parameter PID2. *>> 2000-09-07 DFRENL Krogh Changed text on x to one on ax. *>> 1996-01-08 DFRENL WV Snyder Use DCSPXX for cos(Pi/2 x**2), etc. *>> 1995-11-03 DFRENL Krogh Removed blanks in numbers for C conversion. *>> 1994-11-02 DFRENL Krogh Changes to use M77CON *>> 1994-10-18 DFRENL WV Snyder More specializing instructions *>> 1993-02-25 DFRENL CLL. Edited to eliminate ENTRY and EQUIVALENCE. *>> 1992-09-15 DFRENL WV Snyder Specializing instructions *>> 1992-04-13 DFRENL WV Snyder Declare DFRENF, DFRENG, DFRENS *>> 1992-03-18 DFRENL WV Snyder Move declarations for coefficient arrays *>> 1992-01-24 DFRENL WV Snyder Original code *--D replaces "?": ?FRENC, ?FREN1, ?FRENF, ?FRENG, ?FRENS, ?FRENL *--& ?CSPXX, ?SNPXX * Subprograms in this file compute the Fresnel Cosine and Sine * integrals C(x) and S(x), and the auxiliary functions f(x) and g(x), * for any X: * DFRENC(X) for Fresnel integral C(X) * DFRENS(X) for Fresnel integral S(x) * DFRENF(X) for Fresnel integral auxiliary function f(x) * DFRENG(X) for Fresnel integral auxiliary function g(x). * * Developed by W. V. Snyder, Jet Propulsion Laboratory, 24 January 1992. * * Ref: W. J. Cody, "Chebyshev Approximations for the Fresnel Integrals", * Mathematics of Computation, 1968, pp 450-453 plus Microfiche Suppl. * W. V. Snyder, "Algorithm 723: Fresnel Integrals," ACM Trans. Math. * Softw. 19, 4 (December 1993) 452-456. * Accuracies of highest order formulae, where E is relative error: * * Range Function -log10(E) Function -log10(E) * |X|<=1.2 C(x) 16.24 S(x) 17.26 * 1.2<|X|<=1.6 C(x) 17.47 S(x) 18.66 * 1.6<|X|<=1.9 f(x) 17.13 g(x) 16.25 * 1.9<|X|<=2.4 f(x) 16.64 g(x) 15.65 * 2.4<|X| f(x) 16.89 g(x) 15.58 * * Refer to Cody for accuracy of other approximations. * * ================================================================== */ double /*FUNCTION*/ dfrenc( double x) { double dfrenc_v; dfrenc_v = dfren1( 1, x ); return( dfrenc_v ); } /* end of function */ /* ================================================================== */ double /*FUNCTION*/ dfrenf( double x) { double dfrenf_v; dfrenf_v = dfren1( 3, x ); return( dfrenf_v ); } /* end of function */ /* ================================================================== */ double /*FUNCTION*/ dfreng( double x) { double dfreng_v; dfreng_v = dfren1( 4, x ); return( dfreng_v ); } /* end of function */ /* ================================================================== */ double /*FUNCTION*/ dfrens( double x) { double dfrens_v; dfrens_v = dfren1( 2, x ); return( dfrens_v ); } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define RPI 0.3183098861837906715377675267450287240689e0 #define RPISQ (RPI*RPI) /* end of PARAMETER translations */ double /*FUNCTION*/ dfren1( long mode, double x) { LOGICAL32 wantc, wantf, wantg; long int _l0; double ax, cx, dfren1_v, sx, x4; static double largef, largeg, largex; static double bigx = -1.0e0; static double c = 0.0e0; static double f = 0.5e0; static double g = 0.5e0; static double s = 0.0e0; static double xsave = 0.0e0; static LOGICAL32 havec = TRUE; static LOGICAL32 havef = TRUE; static LOGICAL32 haveg = TRUE; static LOGICAL32 haves = TRUE; static LOGICAL32 needc[8]={TRUE,FALSE,TRUE,TRUE,TRUE,FALSE,FALSE, FALSE}; static LOGICAL32 needs[8]={FALSE,TRUE,TRUE,TRUE,FALSE,TRUE,FALSE, FALSE}; static LOGICAL32 needf[8]={FALSE,FALSE,TRUE,FALSE,TRUE,TRUE,TRUE, FALSE}; static LOGICAL32 needg[8]={FALSE,FALSE,FALSE,TRUE,TRUE,TRUE,FALSE, TRUE}; static double pc1[4-(0)+1]={9.999999999999999421e-1,-1.994608988261842706e-1, 1.761939525434914045e-2,-5.280796513726226960e-4,5.477113856826871660e-6}; static double qc1[4-(1)+1]={4.727921120104532689e-2,1.099572150256418851e-3, 1.552378852769941331e-5,1.189389014228757184e-7}; static double pc2[5-(0)+1]={1.00000000000111043640e0,-2.07073360335323894245e-1, 1.91870279431746926505e-2,-6.71376034694922109230e-4,1.02365435056105864908e-5, -5.68293310121870728343e-8}; static double qc2[5-(1)+1]={3.96667496952323433510e-2,7.88905245052359907842e-4, 1.01344630866749406081e-5,8.77945377892369265356e-8,4.41701374065009620393e-10}; static double ps1[4-(0)+1]={5.2359877559829887021e-1,-7.0748991514452302596e-2, 3.8778212346368287939e-3,-8.4555728435277680591e-5,6.7174846662514086196e-7}; static double qs1[4-(1)+1]={4.1122315114238422205e-2,8.1709194215213447204e-4, 9.6269087593903403370e-6,5.9528122767840998345e-8}; static double ps2[5-(0)+1]={5.23598775598344165913e-1,-7.37766914010191323867e-2, 4.30730526504366510217e-3,-1.09540023911434994566e-4,1.28531043742724820610e-6, -5.76765815593088804567e-9}; static double qs2[5-(1)+1]={3.53398342767472162540e-2,6.18224620195473216538e-4, 6.87086265718620117905e-6,5.03090581246612375866e-8,2.05539124458579596075e-10}; static double pf1[5-(0)+1]={3.1830975293580985290e-1,1.2226000551672961219e1, 1.2924886131901657025e2,4.3886367156695547655e2,4.1466722177958961672e2, 5.6771463664185116454e1}; static double qf1[5-(1)+1]={3.8713003365583442831e1,4.1674359830705629745e2, 1.4740030733966610568e3,1.5371675584895759916e3,2.9113088788847831515e2}; static double pf2[5-(0)+1]={3.183098818220169217e-1,1.958839410219691002e1, 3.398371349269842400e2,1.930076407867157531e3,3.091451615744296552e3, 7.177032493651399590e2}; static double qf2[5-(1)+1]={6.184271381728873709e1,1.085350675006501251e3, 6.337471558511437898e3,1.093342489888087888e4,3.361216991805511494e3}; static double pf3[6-(0)+1]={-9.675460329952532343e-2,-2.431275407194161683e1, -1.947621998306889176e3,-6.059852197160773639e4,-7.076806952837779823e5, -2.417656749061154155e6,-7.834914590078317336e5}; static double qf3[6-(1)+1]={2.548289012949732752e2,2.099761536857815105e4, 6.924122509827708985e5,9.178823229918143780e6,4.292733255630186679e7, 4.803294784260528342e7}; static double pg1[5-(0)+1]={1.013206188102747985e-1,4.445338275505123778e0, 5.311228134809894481e1,1.991828186789025318e2,1.962320379716626191e2, 2.054214324985006303e1}; static double qg1[5-(1)+1]={4.539250196736893605e1,5.835905757164290666e2, 2.544731331818221034e3,3.481121478565452837e3,1.013794833960028555e3}; static double pg2[5-(0)+1]={1.01321161761804586e-1,7.11205001789782823e0, 1.40959617911315524e2,9.08311749529593938e2,1.59268006085353864e3, 3.13330163068755950e2}; static double qg2[5-(1)+1]={7.17128596939302198e1,1.49051922797329229e3, 1.06729678030580897e4,2.41315567213369742e4,1.15149832376260604e4}; static double pg3[6-(0)+1]={-1.53989733819769316e-1,-4.31710157823357568e1, -3.87754141746378493e3,-1.35678867813756347e5,-1.77758950838029676e6, -6.66907061668636416e6,-1.72590224654836845e6}; static double qg3[6-(1)+1]={2.86733194975899483e2,2.69183180396242536e4, 1.02878693056687506e6,1.62095600500231646e7,9.38695862531635179e7, 1.40622441123580005e8}; /* OFFSET Vectors w/subscript range: 1 to dimension */ LOGICAL32 *const Needc = &needc[0] - 1; LOGICAL32 *const Needf = &needf[0] - 1; LOGICAL32 *const Needg = &needg[0] - 1; LOGICAL32 *const Needs = &needs[0] - 1; /* end of OFFSET VECTORS */ /* MODE = 1 means compute C. * MODE = 2 means compute S. * MODE = 3 means compute F. * MODE = 4 means compute G. * ------------------------------------------------------------------ * Internal variables. * * RPI is the reciprocal of PI. * RPISQ is the reciprocal of PI squared. * AX is abs(x). * BIGX is 1/sqrt(round-off). If X > BIGX then to the working * precision x**2 is an integer (which we assume to be a multiple * of four), so cos(pi/2 * x**2) = 1, and sin(pi/2 * x**2) = 0. * C and S are values of C(x) and S(x), respectively. * CX and SX are cos(pi/2 * ax**2) and sin(pi/2 * ax**2), respectively. * F and G are used to compute f(x) and g(x) when X > 1.6. * HAVEC, HAVEF, HAVEG, HAVES are logical variables that indicate * whether the values stored in C, F, G and S correspond to the * value stored in X. HAVEF indicates we have both F and G when * XSAVE .le. 1.6, and HAVEC indicates we have both C and S when * XSAVE .gt. 1.6. * LARGEF is 1/(pi * underflow). If X > LARGEF then f ~ 0. * LARGEG is cbrt(1/(pi**2 * underflow)). If X > LARGEG then g ~ 0. * LARGEX is 1/sqrt(sqrt(underflow)). If X > LARGEX then f ~ 1/(pi * x) * and g ~ 1/(pi**2 * x**3). * MODE indicates the function to be computed: 1 = C(x), 2 = S(x), * 3 = f(x), 4 = g(x). * NEEDC, NEEDF, NEEDG, NEEDS are arrays indexed by MODE (MODE+4 when * X .gt. 1.6) that indicate what functions are needed. * WANTC indicates whether C and S must be computed from F and G. * WANTF and WANTG indicate we computed F and G on the present call. * XSAVE is the most recently provided value of X. * X4 is either X ** 4 or (1.0/X) ** 4. * If you change the order of approximation, you must change the * declarations and DATA statements for the coefficient arrays, * and the executable statements that evaluate the approximations. * ------------------------------------------------------------------ */ /*++ default digits=16 *++ code for digits <= 3 is inactive * double precision PC1(0:1), QC1(1:1) *++ code for digits > 3 & digits <= 6 is inactive * double precision PC1(0:2), QC1(1:2) *++ code for digits > 6 & digits <= 11 is inactive * double precision PC1(0:3), QC1(1:3) *++ code for digits > 11 is active */ /*++ end *++ code for digits <= 2 is inactive * double precision PC2(0:1), QC2(1:1) *++ code for digits > 2 & digits <= 5 is inactive * double precision PC2(0:2), QC2(1:2) *++ code for digits > 5 & digits <= 8 is inactive * double precision PC2(0:3), QC2(1:3) *++ code for digits > 8 & digits <= 12 is inactive * double precision PC2(0:4), QC2(1:4) *++ code for digits > 12 is active */ /*++ end *++ code for digits <= 3 is inactive * double precision PS1(0:1), QS1(1:1) *++ code for digits > 3 & digits <= 7 is inactive * double precision PS1(0:2), QS1(1:2) *++ code for digits > 7 & digits <= 12 is inactive * double precision PS1(0:3), QS1(1:3) *++ code for digits > 12 is active */ /*++ end *++ code for digits <= 2 is inactive * double precision PS2(0:1), QS2(1:1) *++ code for digits > 2 & digits <= 5 is inactive * double precision PS2(0:2), QS2(1:2) *++ code for digits > 5 & digits <= 9 is inactive * double precision PS2(0:3), QS2(1:3) *++ code for digits > 9 & digits <= 14 is inactive * double precision PS2(0:4), QS2(1:4) *++ code for digits > 14 is active */ /*++ end *++ code for digits <= 5 is inactive * double precision PF1(0:1), QF1(1:1) *++ code for digits > 5 & digits <= 8 is inactive * double precision PF1(0:2), QF1(1:2) *++ code for digits > 8 & digits <= 11 is inactive * double precision PF1(0:3), QF1(1:3) *++ code for digits > 11 & digits <= 14 is inactive * double precision PF1(0:4), QF1(1:4) *++ code for digits > 14 is active */ /*++ end *++ code for digits <= 5 is inactive * double precision PF2(0:1), QF2(1:1) *++ code for digits > 5 & digits <= 8 is inactive * double precision PF2(0:2), QF2(1:2) *++ code for digits > 8 & digits <= 11 is inactive * double precision PF2(0:3), QF2(1:3) *++ code for digits > 11 & digits <= 14 is inactive * double precision PF2(0:4), QF2(1:4) *++ code for digits > 14 is active */ /*++ end *++ code for digits <= 3 is inactive * double precision PF3(0:0) *++ code for digits > 3 & digits <= 6 is inactive * double precision PF3(0:1), QF3(1:1) *++ code for digits > 6 & digits <= 9 is inactive * double precision PF3(0:2), QF3(1:2) *++ code for digits > 9 & digits <= 11 is inactive * double precision PF3(0:3), QF3(1:3) *++ code for digits > 11 & digits <= 13 is inactive * double precision PF3(0:4), QF3(1:4) *++ code for digits > 13 & digits <= 15 is inactive * double precision PF3(0:5), QF3(1:5) *++ code for digits > 15 is active */ /*++ end *++ code for digits <= 4 is inactive * double precision PG1(0:1), QG1(1:1) *++ code for digits > 4 & digits <= 7 is inactive * double precision PG1(0:2), QG1(1:2) *++ code for digits > 7 & digits <= 10 is inactive * double precision PG1(0:3), QG1(1:3) *++ code for digits > 10 & digits <= 13 is inactive * double precision PG1(0:4), QG1(1:4) *++ code for digits > 13 is active */ /*++ end *++ code for digits <= 4 is inactive * double precision PG2(0:1), QG2(1:1) *++ code for digits > 4 & digits <= 7 is inactive * double precision PG2(0:2), QG2(1:2) *++ code for digits > 7 & digits <= 10 is inactive * double precision PG2(0:3), QG2(1:3) *++ code for digits > 10 & digits <= 13 is inactive * double precision PG2(0:4), QG2(1:4) *++ code for digits > 13 is active */ /*++ end *++ code for digits <= 3 is inactive * double precision PG3(0:0) *++ code for digits > 3 & digits <= 5 is inactive * double precision PG3(0:1), QG3(r:1) *++ code for digits > 5 & digits <= 8 is inactive * double precision PG3(0:2), QG3(1:2) *++ code for digits > 8 & digits <= 10 is inactive * double precision PG3(0:3), QG3(1:3) *++ code for digits > 10 & digits <= 12 is inactive * double precision PG3(0:4), QG3(1:4) *++ code for digits = 13 is inactive * double precision PG3(0:5), QG3(1:5) *++ code for digits > 13 is active */ /*++ end * */ /* C(x) S(x) f(x) g(x) C(x) S(x) f(x) g(x) */ /* Coefficients for C(x), |X| <= 1.2 *++ code for digits <= 3 is inactive * data pc1 / 1.00053d0, -1.12353d-1/, qc1 /1.38937d-1/ *++ code for digits > 3 & digits <= 6 is inactive * data pc1 / 9.99999896d-1, -1.63090954d-1, 1.06388604d-2/ * data qc1 / 8.36467414d-2, 3.10155884d-3/ *++ code for digits > 6 & digits <= 11 is inactive * data pc1 / 1.0000000000042d0, -1.8651127631106d-1, * * 1.5065663274457d-2, -3.1058074693185d-4/ * data qc1 / 6.0228833908128d-2, 1.7410300558198d-3, * * .6308617899210d-5/ *++ code for digits > 11 is active */ /*++ end * * Coefficients for C(x), 1.2 < |X| <= 1.6 *++ code for digits <= 2 is inactive * data pc2 / 1.3139d0, -5.8827d-2/, qc2 / 4.7324d-1/ *++ code for digits > 2 & digits <= 5 is inactive * data pc2 / 9.9790890d-1, -1.5391895d-1, 9.9817933d-3/ * data qc2 / 8.9878647d-2, 5.6003071d-3/ *++ code for digits > 5 & digits <= 8 is inactive * data pc2 / 1.00000440109d0, -1.83413851908d-1, * * 1.45776170828d-2,-2.78980705270d-4/ * data qc2 / 6.33347445637d-2, 2.01245738369d-3, * * 4.03949004646d-5/ *++ code for digits > 8 & digits <= 12 is inactive * data pc2 / 9.999999966496876d-1, -1.980300987022688d-1, * * 1.735518748450023d-2, -5.069442297935788d-4, * * 5.059962254678234d-6/ * data qc2 / 4.871000308997918d-2, 1.188406974004084d-3, * * 1.824527635843850d-5, 1.678314257874103d-7/ *++ code for digits > 12 is active */ /*++ end * * Coefficients for S(x), |X| <= 1.2 *++ code for digits <= 3 is inactive * data ps1 / 5.23677d-1, -4.67900d-2/, qs1 / 8.81622d-2/ *++ code for digits > 3 & digits <= 7 is inactive * data ps1 / 5.235987665d-1,-5.837763961d-2, 2.165920196d-3/ * data qs1 / 6.474944765d-2, 1.713287588d-3/ *++ code for digits > 7 & digits <= 12 is inactive * data ps1 / 5.23598775598566d-1, -6.59149581139046d-2, * * 3.21501649828293d-3, -4.88704436240178d-5/ * data qs1 / 5.03546388670085d-2, 1.17835980356588d-3, * * 1.37089875826980d-5/ *++ code for digits > 12 is active */ /*++ end * * coefficients for S(x), 1.2 < |X| <= 1.6 *++ code for digits <= 2 is inactive * data ps2 / 5.4766d-1, -3.7151d-2/, qs2 / 1.4559d-1/ *++ code for digits > 2 & digits <= 5 is inactive * data ps2 / 5.2343994d-1, -5.4828347d-2, 1.9020881d-3/ * data qs2 / 7.1104704d-2, 2.5596274d-3/ *++ code for digits > 5 & digits <= 9 is inactive * data ps2 / 5.23599040498d-1, -6.46593392426d-2, * * 3.08030794361d-3, -4.40800854418d-5/ * data qs2 / 5.27536681685d-2, 1.34311026821d-3, * * 1.90476612849d-5/ *++ code for digits > 9 & digits <= 14 is inactive * data ps2 / 5.2359877543509178d-1,-7.0149076634833662d-2, * * 3.8031581605987038d-3,-8.0964948714408156d-5, * * 6.1908080210052772d-7/ * data qs2 / 4.2268067370395487d-2, 8.7642753831073237d-4, * * 1.1088542889789282d-5, 7.8725829545478464d-8/ *++ code for digits > 14 is active */ /*++ end * * coefficients for f(x), 1.6 < |X| <= 1.9 *++ code for digits <= 5 is inactive * data pf1 / 3.1803519d-1, 4.2352332d-1/, qf1 / 1.6015403d0/ *++ code for digits > 5 & digits <= 8 is inactive * data pf1 / 3.18285252094d-1, 2.02860275713d0, * * 1.08108136048d0 / * data qf1 / 6.67171548135d0, 4.52997553972d0 / *++ code for digits > 8 & digits <= 11 is inactive * data pf1 / 3.1830646311448d-1, 4.6077249266684d0, * * 1.1914578438074d1, 3.5555073665830d0 / * data qf1 / 1.4778464938936d1, 4.0902971163705d1, * * 1.6130432784819d1 / *++ code for digits > 11 & digits <= 14 is inactive * data pf1 / 3.1830926850490599d-1, 8.0358812280394156d0, * * 4.8034065557792487d1, 6.9853426160102065d1, * * 1.3530423554038784d1 / * data qf1 / 2.5549161843579529d1, 1.5761100558012250d2, * * 2.4956199380517229d2, 6.5563064008391554d1 / *++ code for digits > 14 is active */ /*++ end * * coefficients for f(x), 1.9 < |X| <= 2.4 *++ code for digits <= 5 is inactive * data pf2 / 3.1825112d-1, 5.8395951d-1/, qf2 / 2.1243944d0 / *++ code for digits > 5 & digits <= 8 is inactive * data pf2 / 3.1830699932d-1, 2.9993457087d0, 2.3956340010d0/ * data qf2 / 9.7254517756d0, 9.4814077696d0 / *++ code for digits > 8 & digits <= 11 is inactive * data pf2 / 3.1830963944521d-1, 7.0770431878327d0, * * 2.8756083262903d1, 1.3583685742326d1 / * data qf2 / 2.2536994405207d1, 9.6127568469278d1, * * 5.7407028004031d1 / *++ code for digits > 11 & digits <= 14 is inactive * data pf2 / 3.183098568640159d-1, 1.265129469683175d1, * * 1.219467616498339d2, 2.910033655512762d2, * * 9.278397828631516d1/ * data qf2 / 4.004915302781009d1, 3.942059697951583d2, * * 1.001368403495691d3, 4.142676224222433d2/ *++ code for digits > 14 is active */ /*++ end * * coefficients for f(x), 2.4 < |X| *++ code for digits <= 3 is inactive * data pf3 /-8.97969d-2/ *++ code for digits > 3 & digits <= 6 is inactive * data pf3 / -9.67122165d-2, -3.73565920d-1/ * data qf3 / 7.29466572d0 / *++ code for digits > 6 & digits <= 9 is inactive * data pf3 /-9.67541443648d-2, -2.26566293818d0, * * -3.53429821084d0 / * data qf3 / 2.69596939977d1, 9.72883957469d1 / *++ code for digits > 9 & digits <= 11 is inactive * data pf3 /-9.6754596389025d-2, -5.5254943840897d0, * * -5.8606282086171d1, -5.2235560918394d1 / * data qf3 / 6.0654483020979d1, 7.8528841711294d2, * * 1.8592498578831d3 / *++ code for digits > 11 & digits <= 13 is inactive * data pf3 /-9.675460316952504d-2, -1.023876428129288d1, * * -2.712579634037998d2, -1.766119345282127d3, * * -1.043464426656267d3 / * data qf3 / 1.093682244053459d2, 3.155843461920557d3, * * 2.625634316204420d4, 4.578252057246393d4 / *++ code for digits > 13 & digits <= 15 is inactive * data pf3 /-9.67546032967090380d-2,-1.64797712841245767d1, * * -8.16343401784374598d2, -1.34922028171857248d4, * * -6.13547113614699772d4, -2.61294753225141779d4/ * data qf3 / 1.73871690673649114d2, 9.01827596231524147d3, * * 1.65946462621853184d5, 1.00105478900791339d6, * * 1.37012364817225972d6 / *++ code for digits > 15 is active */ /*++ end * * coefficients for g(x), 1.6 < |X| <= 1.9 *++ code for digits <= 4 is inactive * data pg1 / 1.007011d-1, 1.457780d-1/, qg1 / 2.700253d0 / *++ code for digits > 4 & digits <= 7 is inactive * data pg1 / 1.012500483d-1, 7.735207446d-1, 3.878282770d-1/ * data qg1 / 9.101956178d0, 1.002234185d1 / *++ code for digits > 7 & digits <= 10 is inactive * data pg1 / 1.013095436817d-1, 1.731937984173d0, * * 5.036588245265d0, 1.290707246507d0 / * data qg1 / 1.860077430076d1, 6.901951935725d1, * * 4.308918659989d1 / *++ code for digits > 10 & digits <= 13 is inactive * data pg1 / 1.013188096509180d-1, 2.966220554725899d0, * * 2.015435299505393d1, 3.155605679387908d1, * * 4.916012830646366d0 / * data qg1 / 3.079178736724045d1, 2.362874318980473d2, * * 4.953861258248338d2, 2.026144493403599d2 / *++ code for digits > 13 is active */ /*++ end * * coefficients for g(x), 1.9 < |X| <= 2.4 *++ code for digits <= 4 is inactive * data pg2 / 1.011711d-1, 2.239630d-1/, qg2 / 3.609237d0 / *++ code for digits > 4 & digits <= 7 is inactive * data pg2 / 1.013115463d-1, 1.182021191d0, 9.886315969d-1/ * data qg2 / 1.317219285d1, 2.101104994d1 / *++ code for digits > 7 & digits <= 10 is inactive * data pg2 / 1.013202025653d-1, 2.690767013770d0, * * 1.283624749271d1, 5.791337587723d0 / * data qg2 / 2.807457005500d1, 1.598704313522d2, * * 1.525630385045d2 / *++ code for digits > 10 & digits <= 13 is inactive * data pg2 / 1.013210509409046d-1, 4.682769769757399d0, * * 5.241351613346472d1, 1.411735944550041d2, * * 4.019756012788710d1 / * data qg2 / 4.773652966704634d1, 5.802036967947208d2, * * 1.947179327244406d3, 1.266041236960445d3 / *++ code for digits > 13 is active */ /*++ end * * coefficients for g(x), 2.4 < |X| *++ code for digits <= 3 is inactive * data pg3 /-1.3549d-1/ *++ code for digits > 3 & digits <= 5 is inactive * data pg3 /-1.5382824d-1, -6.4547464d-1/, qg3 / 1.0290328d1/ *++ code for digits > 5 & digits <= 8 is inactive * data pg3 /-1.5398760302d-1, -4.2772857997d0, * * -6.5940181141d0 / * data qg3 / 3.4149986477d1, 1.7071708816d2 / *++ code for digits > 8 & digits <= 10 is inactive * data pg3 /-1.539896971616d-1, -1.024638314446d1, * * -1.252507402120d2, -1.029189676144d2 / * data qg3 / 7.292229303754d1, 1.186528673248d3, * * 3.844430908476d3 / *++ code for digits > 10 & digits <= 12 is inactive * data pg3 /-1.53989733057971d-1, -1.86483223831639d1, * * -5.66882778026550d2, -4.16714647017489d3, * * -2.14678074364341d3 / * data qg3 / 1.27484298075273d2, 4.40258858159927d3, * * 4.57583830694684d4, 1.08207883281275d5 / *++ code for digits = 13 is inactive * data pg3 /-1.539897338019247d-1, -2.959078231855258d1, * * -1.661867087183632d3, -3.112796454920657d4, * * -1.573536286819766d5, -5.574884113746041d4 / * data qg3 / 1.985439813544440d2, 1.196693134508007d4, * * 2.625573864512749d5, 1.969055829311119d6, * * 3.629081333313312d6 / *++ code for digits > 13 is active */ /*++ end * ------------------------------------------------------------------ */ if (bigx < 0.0e0) { bigx = 1.0e0/DBL_EPSILON; largef = RPI/DBL_MIN; largeg = pow(RPI*largef,1.0e0/3.0e0); largex = 1.0e0/sqrt( sqrt( DBL_MIN ) ); } if (x != xsave) { havec = FALSE; havef = FALSE; haveg = FALSE; haves = FALSE; } ax = fabs( x ); if (ax <= 1.6e0) { x4 = powi(ax,4); if (Needc[mode] && !havec) { if (ax <= 1.2e0) { /*++ code for digits <= 3 is inactive * c = x * (pc1(1)*x4+pc1(0)) / (qc1(1)*x4+1.0d0) *++ code for digits > 3 & digits <= 6 is inactive * c = x * ((pc1(2)*x4+pc1(1))*x4+pc1(0)) * 1 / ((qc1(2)*x4+qc1(1))*x4+1.0d0) *++ code for digits > 6 & digits <= 11 is inactive * c = x * (((pc1(3)*x4+pc1(2))*x4+pc1(1))*x4+pc1(0)) * 1 / (((qc1(3)*x4+qc1(2))*x4+qc1(1))*x4+1.0d0) *++ code for digits > 11 is active */ c = x*((((pc1[4]*x4 + pc1[3])*x4 + pc1[2])*x4 + pc1[1])* x4 + pc1[0])/((((qc1[3]*x4 + qc1[2])*x4 + qc1[1])* x4 + qc1[0])*x4 + 1.0e0); /*++ end */ } else { /*++ code for digits <= 2 is inactive * c = x * (pc2(1)*x4+pc2(0)) / (qc2(1)*x4+1.0d0) *++ code for digits > 2 & digits <= 5 is inactive * c = x * ((pc2(2)*x4+pc2(1))*x4+pc2(0)) * 1 / ((qc2(2)*x4+qc2(1))*x4+1.0d0) *++ code for digits > 5 & digits <= 8 is inactive * c = x * (((pc2(3)*x4+pc2(2))*x4+pc2(1))*x4+pc2(0)) * 1 / (((qc2(3)*x4+qc2(2))*x4+qc2(1))*x4+1.0d0) *++ code for digits > 8 & digits <= 12 is inactive * c = x * ((((pc2(4)*x4+pc2(3))*x4+pc2(2))*x4+pc2(1))*x4+ * 1 pc2(0)) * 2 / ((((qc2(4)*x4+qc2(3))*x4+qc2(2))*x4+qc2(1))*x4+ * 3 1.0d0) *++ code for digits > 12 is active */ c = x*(((((pc2[5]*x4 + pc2[4])*x4 + pc2[3])*x4 + pc2[2])* x4 + pc2[1])*x4 + pc2[0])/(((((qc2[4]*x4 + qc2[3])* x4 + qc2[2])*x4 + qc2[1])*x4 + qc2[0])*x4 + 1.0e0); /*++ end */ } havec = TRUE; } if (Needs[mode] && !haves) { if (ax <= 1.2e0) { /*++ code for digits <= 3 is inactive * s = x**3*(ps1(1)*x4+ps1(0)) / (qs1(1)*x4+1.0d0) *++ code for digits > 3 & digits <= 7 is inactive * s = x**3*((ps1(2)*x4+ps1(1))*x4+ps1(0)) * 1 / ((qs1(2)*x4+qs1(1))*x4+1.0d0) *++ code for digits > 7 & digits <= 12 is inactive * s = x**3*(((ps1(3)*x4+ps1(2))*x4+ps1(1))*x4+ps1(0)) * 1 / (((qs1(3)*x4+qs1(2))*x4+qs1(1))*x4+1.0d0) *++ code for digits > 12 is active */ s = CUBE(x)*((((ps1[4]*x4 + ps1[3])*x4 + ps1[2])*x4 + ps1[1])*x4 + ps1[0])/((((qs1[3]*x4 + qs1[2])*x4 + qs1[1])*x4 + qs1[0])*x4 + 1.0e0); /*++ end */ } else { /*++ code for digits <= 2 is inactive * s = x**3*(ps2(1)*x4+ps2(0)) / (qs2(1)*x4+1.0d0) *++ code for digits > 2 & digits <= 5 is inactive * s = x**3*((ps2(2)*x4+ps2(1))*x4+ps2(0)) * 1 / ((qs2(2)*x4+qs2(1))*x4+1.0d0) *++ code for digits > 5 & digits <= 9 is inactive * s = x**3*(((ps2(3)*x4+ps2(2))*x4+ps2(1))*x4+ps2(0)) * 1 / (((qs2(3)*x4+qs2(2))*x4+qs2(1))*x4+1.0d0) *++ code for digits > 9 & digits <= 14 is inactive * s = x**3*((((ps2(4)*x4+ps2(3))*x4+ps2(2))*x4+ps2(1))*x4+ * 1 ps2(0)) * 2 / ((((qs2(4)*x4+qs2(3))*x4+qs2(2))*x4+qs2(1))*x4+ * 3 1.0d0) *++ code for digits > 14 is active */ s = CUBE(x)*(((((ps2[5]*x4 + ps2[4])*x4 + ps2[3])* x4 + ps2[2])*x4 + ps2[1])*x4 + ps2[0])/(((((qs2[4]* x4 + qs2[3])*x4 + qs2[2])*x4 + qs2[1])*x4 + qs2[0])* x4 + 1.0e0); /*++ end */ } haves = TRUE; } if ((Needf[mode] || Needg[mode]) && !havef) { cx = dcspxx( ax ); sx = dsnpxx( ax ); f = (0.5e0 - s)*cx - (0.5e0 - c)*sx; g = (0.5e0 - c)*cx + (0.5e0 - s)*sx; havef = TRUE; } } else { if (ax <= largex) { x4 = powi(1.0e0/ax,4); wantf = Needf[mode + 4] && !havef; if (wantf) { if (ax <= 1.9e0) { /*++ code for digits <= 5 is inactive * f = (pf1(1)*x4+pf1(0)) / ((qf1(1)*x4+1.0d0) * ax) *++ code for digits > 5 & digits <= 8 is inactive * f = ((pf1(2)*x4+pf1(1))*x4+pf1(0)) * 1 / (((qf1(2)*x4+qf1(1))*x4+1.0d0) * ax) *++ code for digits > 8 & digits <= 11 is inactive * f = (((pf1(3)*x4+pf1(2))*x4+pf1(1))*x4+pf1(0)) * 1 / ((((qf1(3)*x4+qf1(2))*x4+qf1(1))*x4+1.0d0) * ax) *++ code for digits > 11 & digits <= 14 is inactive * f = ((((pf1(4)*x4+pf1(3))*x4+pf1(2))*x4+pf1(1))*x4+ * 1 pf1(0)) * 2 / (((((qf1(4)*x4+qf1(3))*x4+qf1(2))*x4+qf1(1))*x4+ * 3 1.0d0) * ax) *++ code for digits > 14 is active */ f = (((((pf1[5]*x4 + pf1[4])*x4 + pf1[3])*x4 + pf1[2])*x4 + pf1[1])*x4 + pf1[0])/((((((qf1[4]* x4 + qf1[3])*x4 + qf1[2])*x4 + qf1[1])*x4 + qf1[0])* x4 + 1.0e0)*ax); /*++ end */ } else if (ax <= 2.4) { /*++ code for digits <= 5 is inactive * f = (pf2(1)*x4+pf2(0)) / ((qf2(1))*x4+1.0d0) * ax) *++ code for digits > 5 & digits <= 8 is inactive * f = ((pf2(2)*x4+pf2(1))*x4+pf2(0)) * 1 / (((qf2(2)*x4+qf2(1))*x4+1.0d0) * ax) *++ code for digits > 8 & digits <= 11 is inactive * f = (((pf2(3)*x4+pf2(2))*x4+pf2(1))*x4+pf2(0)) * 1 / ((((qf2(3)*x4+qf2(2))*x4+qf2(1))*x4+1.0d0) * ax) *++ code for digits > 11 & digits <= 14 is inactive * f = ((((pf2(4)*x4+pf2(3))*x4+pf2(2))*x4+pf2(1))*x4+ * 1 pf2(0)) * 2 / (((((qf2(4)*x4+qf2(3))*x4+qf2(2))*x4+qf2(1))*x4+ * 3 1.0d0) * ax) *++ code for digits > 14 is active */ f = (((((pf2[5]*x4 + pf2[4])*x4 + pf2[3])*x4 + pf2[2])*x4 + pf2[1])*x4 + pf2[0])/((((((qf2[4]* x4 + qf2[3])*x4 + qf2[2])*x4 + qf2[1])*x4 + qf2[0])* x4 + 1.0e0)*ax); /*++ end */ } else { /*++ code for digits <= 3 is inactive * f = (rpi + x4*pf3(0)) / ax *++ code for digits > 3 & digits <= 6 is inactive * f = (rpi + x4*(pf3(1)*x4+pf3(0)) / (qf3(1)*x4+1.0d0)) * / ax *++ code for digits > 6 & digits <= 9 is inactive * f = (rpi + * 1 x4*((pf3(2)*x4+pf3(1))*x4+pf3(0)) * 2 / ((qf3(2)*x4+qf3(1))*x4+1.0d0)) / ax *++ code for digits > 9 & digits <= 11 is inactive * f = (rpi + * 1 x4*(((pf3(3)*x4+pf3(2))*x4+pf3(1))*x4+pf3(0)) * 2 / (((qf3(3)*x4+qf3(2))*x4+qf3(1))*x4+1.0d0)) / ax *++ code for digits > 11 & digits <= 13 is inactive * f = (rpi + * 1 x4*((((pf3(4)*x4+pf3(3))*x4+pf3(2))*x4+pf3(1))*x4+ * 2 pf3(0)) * 3 / ((((qf3(4)*x4+qf3(3))*x4+qf3(2))*x4+qf3(1))*x4+ * 4 1.0d0)) / ax *++ code for digits > 13 & digits <= 15 is inactive * f = (rpi + * 1 x4*(((((pf3(5)*x4+pf3(4))*x4+pf3(3))*x4+pf3(2))*x4+ * 2 pf3(1))*x4+pf3(0)) * 3 / (((((qf3(5)*x4+qf3(4))*x4+qf3(3))*x4+qf3(2))*x4+ * 4 qf3(1))*x4+1.0d0)) / ax *++ code for digits > 15 is active */ f = (RPI + x4*((((((pf3[6]*x4 + pf3[5])*x4 + pf3[4])* x4 + pf3[3])*x4 + pf3[2])*x4 + pf3[1])*x4 + pf3[0])/ ((((((qf3[5]*x4 + qf3[4])*x4 + qf3[3])*x4 + qf3[2])* x4 + qf3[1])*x4 + qf3[0])*x4 + 1.0e0))/ax; /*++ end */ } havef = TRUE; } wantg = Needg[mode + 4] && !haveg; if (wantg) { if (ax <= 1.9e0) { /*++ code for digits <=4 is inactive * g = (pg1(1)*x4+pg1(0)) / ((qg1(1)*x4+1.0d0) * ax**3) *++ code for digits > 4 & digits <= 7 is inactive * g = ((pg1(2)*x4+pg1(1))*x4+pg1(0)) * 1 / (((qg1(2)*x4+qg1(1))*x4+1.0d0) * ax**3) *++ code for digits > 7 & digits <= 10 is inactive * g = (((pg1(3)*x4+pg1(2))*x4+pg1(1))*x4+pg1(0)) * 1 / ((((qg1(3)*x4+qg1(2))*x4+qg1(1))*x4+1.0d0) *ax**3) *++ code for digits > 10 & digits <= 13 is inactive * g = ((((pg1(4)*x4+pg1(3))*x4+pg1(2))*x4+pg1(1))*x4+ * 1 pg1(0)) * 2 / (((((qg1(4)*x4+qg1(3))*x4+qg1(2))*x4+qg1(1))*x4+ * 3 1.0d0) * ax**3) *++ code for digits > 13 is active */ g = (((((pg1[5]*x4 + pg1[4])*x4 + pg1[3])*x4 + pg1[2])*x4 + pg1[1])*x4 + pg1[0])/((((((qg1[4]* x4 + qg1[3])*x4 + qg1[2])*x4 + qg1[1])*x4 + qg1[0])* x4 + 1.0e0)*CUBE(ax)); /*++ end */ } else if (ax <= 2.4e0) { /*++ code for digits <=4 is inactive * g = (pg2(1)*x4+pg2(0)) / ((qg2(1)*x4+1.0d0) * ax**3) *++ code for digits > 4 & digits <= 7 is inactive * g = ((pg2(2)*x4+pg2(1))*x4+pg2(0)) * 1 / (((qg2(2)*x4+qg2(1))*x4+1.0d0) * ax**3) *++ code for digits > 7 & digits <= 10 is inactive * g = (((pg2(3)*x4+pg2(2))*x4+pg2(1))*x4+pg2(0)) * 1 / ((((qg2(3)*x4+qg2(2))*x4+qg2(1))*x4+1.0d0) *ax**3) *++ code for digits > 10 & digits <= 13 is inactive * g = ((((pg2(4)*x4+pg2(3))*x4+pg2(2))*x4+pg2(1))*x4+ * 1 pg2(0)) * 2 / (((((qg2(4)*x4+qg2(3))*x4+qg2(2))*x4+qg2(1))*x4+ * 3 1.0d0) * ax**3) *++ code for digits > 13 is active */ g = (((((pg2[5]*x4 + pg2[4])*x4 + pg2[3])*x4 + pg2[2])*x4 + pg2[1])*x4 + pg2[0])/((((((qg2[4]* x4 + qg2[3])*x4 + qg2[2])*x4 + qg2[1])*x4 + qg2[0])* x4 + 1.0e0)*CUBE(ax)); /*++ end */ } else { /*++ code for digits <=3 is inactive * g = (rpisq + x4*pg3(0)) / ax**3 *++ code for digits > 3 & digits <= 5 is inactive * g = (rpisq + x4*(pg3(1)*x4+pg3(0)) * 1 / (qg3(1)*x4+1.0d0)) / ax**3 *++ code for digits > 5 & digits <= 8 is inactive * g = (rpisq + x4*((pg3(2)*x4+pg3(1))*x4+pg3(0)) * 1 / ((qg3(2)*x4+qg3(1))*x4+1.0d0)) / ax**3 *++ code for digits > 8 & digits <= 10 is inactive * g = (rpisq + * 1 x4*(((pg3(3)*x4+pg3(2))*x4+pg3(1))*x4+pg3(0)) * 2 / (((qg3(3)*x4+qg3(2))*x4+qg3(1))*x4+1.0d0)) * 3 / ax**3 *++ code for digits > 10 & digits <= 12 is inactive * g = (rpisq + * 1 x4*((((pg3(4)*x4+pg3(3))*x4+pg3(2))*x4+pg3(1))*x4+ * 2 pg3(0)) * 3 / ((((qg3(4)*x4+qg3(3))*x4+qg3(2))*x4+qg3(1))*x4+ * 4 1.0d0)) / ax**3 *++ code for digits = 13 is inactive * g = (rpisq + * 1 x4*(((((pg3(5)*x4+pg3(4))*x4+pg3(3))*x4+pg3(2))*x4+ * 2 pg3(1))*x4+pg3(0)) * 3 / (((((qg3(5)*x4+qg3(4))*x4+qg3(3))*x4+qg3(2))*x4+ * 4 qg3(1))*x4+1.0d0)) / ax**3 *++ code for digits > 13 is active */ g = (RPISQ + x4*((((((pg3[6]*x4 + pg3[5])*x4 + pg3[4])*x4 + pg3[3])*x4 + pg3[2])*x4 + pg3[1])* x4 + pg3[0])/((((((qg3[5]*x4 + qg3[4])*x4 + qg3[3])* x4 + qg3[2])*x4 + qg3[1])*x4 + qg3[0])*x4 + 1.0e0))/ CUBE(ax); /*++ end */ } haveg = TRUE; } } else { wantf = Needf[mode]; if (wantf) { if (x <= largef) { f = RPI/ax; } else { f = 0.0e0; } } wantg = Needg[mode]; if (wantg) { if (x <= largeg) { g = RPISQ/CUBE(ax); } else { g = 0.0e0; } } } wantc = (Needc[mode + 4] || Needs[mode + 4]) && !havec; if (wantc || x < 0.0e0) { cx = dcspxx( ax ); sx = dsnpxx( ax ); if (wantc) { c = 0.5e0 + f*sx - g*cx; s = 0.5e0 - f*cx - g*sx; if (x < 0.0) { c = -c; s = -s; } havec = TRUE; } if (x < 0.0) { /* We COULD do the following before the preceeding, and then * not put in a test in the preceeding for x .lt. 0, but * even though the results are mathematically identical, we * would have some cancellation above if we did so. */ if (wantg) g = cx + sx - g; if (wantf) f = cx - sx - f; } } } xsave = x; switch (mode) { case 1: goto L_1; case 2: goto L_2; case 3: goto L_3; case 4: goto L_4; } L_1: dfren1_v = c; return( dfren1_v ); L_2: dfren1_v = s; return( dfren1_v ); L_3: dfren1_v = f; return( dfren1_v ); L_4: dfren1_v = g; return( dfren1_v ); } /* end of function */