Source Code 1992 March

home *** CD-ROM | disk | FTP | other *** search

/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2331 < prev next >

Wrap

Internet Message Format | 1990-12-28 | 54.9 KB

From: glenn@qed.physics.su.OZ.AU (Glenn Geers) Newsgroups: alt.sources Subject: 80386 alternative math lib v2.0 part02/06 Message-ID: <1990Dec16.210128.27937@metro.ucc.su.OZ.AU> Date: 16 Dec 90 21:01:28 GMT Submitted-by: root@trantor Archive-name: mathlib2.0/part02 ---- Cut Here and feed the following to sh ---- #!/bin/sh # this is mathlib.02 (part 2 of mathlib2.0) # do not concatenate these parts, unpack them in order with /bin/sh # file erf.c continued # if test ! -r _shar_seq_.tmp; then echo 'Please unpack part 1 first!' exit 1 fi (read Scheck if test "$Scheck" != 2; then echo Please unpack part "$Scheck" next! exit 1 else exit 0 fi ) < _shar_seq_.tmp || exit 1 if test ! -f _shar_wnt_.tmp; then echo 'x - still skipping erf.c' else echo 'x - continuing file erf.c' sed 's/^X//' << 'SHAR_EOF' >> 'erf.c' && X xden = ysq; X for (i = 0; i <= 3; i++) { X xnum = (xnum+P[i])*ysq; X xden = (xden+Q[i])*ysq; X } X result = ysq*(xnum+P[4])/(xden+Q[4]); X result = (SQRPI-result)/y; X if (jint != 2) { X i = (int)(y*SIXTEN); ysq = (double)i/SIXTEN; X result *= exp(-ysq*ysq)*exp(-(y-ysq)*(y+ysq)); X } X } X } X if (jint == 0) { /* Fix up for negative argument, erf, etc. */ X result = HALF-result; result += HALF; X if (x < ZERO) X result = -result; X } X else if (jint == 1) { X if (x < ZERO) X result = TWO-result; X } X else if (x < ZERO) { X if (x < XNEG) X result = XINF; X else { X i = (int)(x*SIXTEN); ysq = (double)i/SIXTEN; X y = exp(ysq*ysq)*exp((x-ysq)*(x+ysq)); X result = -result; result += y+y; X } X } X return result; } X /* X * This subprogram computes approximate values for erf(x). X * (see comments heading calerf()). X * X * Author/date: W. J. Cody, January 8, 1985 X */ double #if defined(__STDC__) || defined(__GNUC__) erf(double x) #else erf(x) double x; #endif { X return calerf(x,0); } X /* X * This subprogram computes approximate values for erfc(x). X * (see comments heading calerf()). X * X * Author/date: W. J. Cody, January 8, 1985 X */ double #if defined(__STDC__) || defined(__GNUC__) erfc(double x) #else erfc(x) double x; #endif { X return calerf(x,1); } X /* X * This subprogram computes approximate values for exp(x*x) * erfc(x). X * (see comments heading calerf()). X * X * Author/date: W. J. Cody, March 30, 1987 X */ double #if defined(__STDC__) || defined(__GNUC__) erfcx(double x) #else erfcx(x) double x; #endif { X return calerf(x,2); } SHAR_EOF echo 'File erf.c is complete' && chmod 0644 erf.c || echo 'restore of erf.c failed' Wc_c="`wc -c < 'erf.c'`" test 9681 -eq "$Wc_c" || echo 'erf.c: original size 9681, current size' "$Wc_c" rm -f _shar_wnt_.tmp fi # ============= nextafter.c ============== if test -f 'nextafter.c' -a X"$1" != X"-c"; then echo 'x - skipping nextafter.c (File already exists)' rm -f _shar_wnt_.tmp else > _shar_wnt_.tmp echo 'x - extracting nextafter.c (Text)' sed 's/^X//' << 'SHAR_EOF' > 'nextafter.c' && /* ** This file is part of the alternative 80386 math library and is ** covered by the GNU General Public license with my modification ** as noted in the README file that accompanied this file. ** ** Copyright 1990 G. Geers ** ** A mix of C and assembler - well I've got the functions so I might ** as well use them! ** */ X #include "fpumath.h" X asm(".align 4"); asm(".Lulp:"); asm(".double 1.1102230246251565e-16"); X asm(".align 4"); asm(".Lulpup:"); asm(".double 2.2204460492503131e-16"); X double nextafter(x, y) double x, y; { X asm("subl $8, %esp"); X X if (isnan(x) || isnan(y)) X return(quiet_nan(1.0)); X X if (isinf(x)) X if (y > x) X return(-max_normal()); X else X if (y < x) X return(max_normal()); X X if (x == 0.0) X if (y > 0.0) X return(min_subnormal()); X else X return(-min_subnormal()); X X if (isnormal(x)) { X if ((x == min_normal()) && y < x) X return(max_subnormal()); X X if ((x == max_normal()) && y > x) X return(infinity()); X X if ((x == -max_normal()) && y < x) X return(-infinity()); X X asm("movl 12(%ebp), %eax"); X asm("andl $0x7ff00000, %eax"); X asm("movl %eax, -12(%ebp)"); X asm("movl $0x0, -16(%ebp)"); X asm("fldl -16(%ebp)"); X X if (fabs(x) <= 2.0 && y < x) X asm("fldl .Lulp"); X else X asm("fldl .Lulpup"); X X asm("fmulp"); X X if (y > x) { X asm("fldl 8(%ebp)"); X asm("faddp"); X asm("leave"); X asm("ret"); X } X if (y < x) { X asm("fldl 8(%ebp)"); X asm("fsubp"); X asm("leave"); X asm("ret"); X } X } X else X if (issubnormal(x)) { X if ((x == max_subnormal()) && y > x) X return(min_normal()); X X if ((x == -max_subnormal()) && y < x) X return(-min_normal()); X X if (y > x) X return(x + min_subnormal()); X X if (y < x) X return(x - min_subnormal()); X } X X return(x); } SHAR_EOF chmod 0644 nextafter.c || echo 'restore of nextafter.c failed' Wc_c="`wc -c < 'nextafter.c'`" test 1725 -eq "$Wc_c" || echo 'nextafter.c: original size 1725, current size' "$Wc_c" rm -f _shar_wnt_.tmp fi # ============= j0f.c ============== if test -f 'j0f.c' -a X"$1" != X"-c"; then echo 'x - skipping j0f.c (File already exists)' rm -f _shar_wnt_.tmp else > _shar_wnt_.tmp echo 'x - extracting j0f.c (Text)' sed 's/^X//' << 'SHAR_EOF' > 'j0f.c' && #ifndef lint static char sccsid[] = "@(#)j0.c 1.2 (ucb.beef) 10/2/89"; #endif /* !defined(lint) */ /* X * This packet computes zero-order Bessel functions of the first and X * second kind (j0 and y0), for real arguments x, where 0 < x <= XMAX X * for y0, and |x| <= XMAX for j0. It contains three function-type X * subprograms, j0, y0 and caljy0. The calling statements for the X * primary entries are: X * X * y = j0(x) X * and X * y = y0(x), X * X * where the entry points correspond to the functions j0(x) and y0(x), X * respectively. The routine caljy0() is intended for internal packet X * use only, all computations within the packet being concentrated in X * this one routine. The function subprograms invoke caljy0 with X * the statement X * X * result = caljy0(x,jint); X * X * where the parameter usage is as follows: X * X * Function Parameters for caljy0 X * call x result jint X * X * j0(x) |x| <= XMAX j0(x) 0 X * y0(x) 0 < x <= XMAX y0(x) 1 X * X * The main computation uses unpublished minimax rational X * approximations for x <= 8.0, and an approximation from the X * book Computer Approximations by Hart, et. al., Wiley and Sons, X * New York, 1968, for arguments larger than 8.0 Part of this X * transportable packet is patterned after the machine-dependent X * FUNPACK program for j0(x), but cannot match that version for X * efficiency or accuracy. This version uses rational functions X * that are theoretically accurate to at least 18 significant decimal X * digits for x <= 8, and at least 18 decimal places for x > 8. The X * accuracy achieved depends on the arithmetic system, the compiler, X * the intrinsic functions, and proper selection of the machine- X * dependent constants. X * X ********************************************************************* X * X * Explanation of machine-dependent constants X * X * XINF = largest positive machine number X * XMAX = largest acceptable argument. The functions sin(), floor() X * and cos() must perform properly for fabs(x) <= XMAX. X * We recommend that XMAX be a small integer multiple of X * sqrt(1/eps), where eps is the smallest positive number X * such that 1+eps > 1. X * XSMALL = positive argument such that 1.0-(x/2)**2 = 1.0 X * to machine precision for all fabs(x) <= XSMALL. X * We recommend that XSMALL < sqrt(eps)/beta, where beta X * is the floating-point radix (usually 2 or 16). X * X * Approximate values for some important machines are X * X * eps XMAX XSMALL XINF X * X * CDC 7600 (S.P.) 7.11E-15 1.34E+08 2.98E-08 1.26E+322 X * CRAY-1 (S.P.) 7.11E-15 1.34E+08 2.98E-08 5.45E+2465 X * IBM PC (8087) (S.P.) 5.96E-08 8.19E+03 1.22E-04 3.40E+38 X * IBM PC (8087) (D.P.) 1.11D-16 2.68D+08 3.72D-09 1.79D+308 X * IBM 195 (D.P.) 2.22D-16 6.87D+09 9.09D-13 7.23D+75 X * UNIVAC 1108 (D.P.) 1.73D-18 4.30D+09 2.33D-10 8.98D+307 X * VAX 11/780 (D.P.) 1.39D-17 1.07D+09 9.31D-10 1.70D+38 X * X ********************************************************************* X ********************************************************************* X * X * Error Returns X * X * The program returns the value zero for x > XMAX, and returns X * -XINF when y0 is called with a negative or zero argument. X * X * X * Intrinsic functions required are: X * X * fabs, floor, cos, log, sin, sqrt X * X * X * Latest modification: June 2, 1989 X * X * Author: W. J. Cody X * Mathematics and Computer Science Division X * Argonne National Laboratory X * Argonne, IL 60439 X */ X X #include "fpumath.h" X X /* Machine-dependent constants */ #if defined(vax) || defined(tahoe) #define XMAX (double)1.07e9 #define XSMALL (double)9.31e-10 #define XINF (double)1.70e38 #else /* defined(vax) || defined(tahoe) */ #define XMAX (double)2.68e8 #define XSMALL (double)3.72e-9 #define XINF MAXFLOAT #endif /* defined(vax) || defined(tahoe) */ X /* Mathematical constants */ #define EIGHT (double)8 #define CONS (double)(-1.1593151565841244881e-1) /* ln(.5)+Euler's gamma */ #define FIVE5 (double)5.5 #define FOUR (double)4 #define ONE (double)1 #define ONEOV8 (double)0.125 #define PI2 (double)6.3661977236758134308e-1 #define P17 (double)1.716e-1 #define SIXTY4 (double)64 #define THREE (double)3 #define TWOPI (double)6.2831853071795864769e0 #define TWOPI1 (double)6.28125 #define TWOPI2 (double)1.9353071795864769253e-03 #define TWO56 (double)256 #define ZERO (double)0 X /* Zeroes of Bessel functions */ #define XJ0 (double)2.4048255576957727686e0 #define XJ1 (double)5.5200781102863106496e0 #define XY0 (double)8.9357696627916752158e-1 #define XY1 (double)3.9576784193148578684e0 #define XY2 (double)7.0860510603017726976e0 #define XJ01 (double)616 #define XJ02 (double)(-1.4244423042272313784e-3) #define XJ11 (double)1413 #define XJ12 (double)5.4686028631064959660e-4 #define XY01 (double)228.0 #define XY02 (double)2.9519662791675215849e-3 #define XY11 (double)1013 #define XY12 (double)6.4716931485786837568e-4 #define XY21 (double)1814 #define XY22 (double)1.1356030177269762362e-4 X /* X * Coefficients for rational approximation to ln(x/a) X */ static double PLG[] = { X -2.4562334077563243311e01, X 2.3642701335621505212e02, X -5.4989956895857911039e02, X 3.5687548468071500413e02, }; static double QLG[] = { X -3.5553900764052419184e01, X 1.9400230218539473193e02, X -3.3442903192607538956e02, X 1.7843774234035750207e02, }; X /* X * Coefficients for rational approximation of X * j0(x) / (x**2 - XJ0**2), XSMALL < |x| <= 4.0 X */ static double PJ0[] = { X 6.6302997904833794242e06, X -6.2140700423540120665e08, X 2.7282507878605942706e10, X -4.1298668500990866786e11, X -1.2117036164593528341e-01, X 1.0344222815443188943e02, X -3.6629814655107086448e04, }; static double QJ0[] = { X 4.5612696224219938200e05, X 1.3985097372263433271e08, X 2.6328198300859648632e10, X 2.3883787996332290397e12, X 9.3614022392337710626e02, }; X /* X * Coefficients for rational approximation of X * j0(x) / (x**2 - XJ1**2), 4.0 < |x| <= 8.0 X */ static double PJ1[] = { X 4.4176707025325087628e03, X 1.1725046279757103576e04, X 1.0341910641583726701e04, X -7.2879702464464618998e03, X -1.2254078161378989535e04, X -1.8319397969392084011e03, X 4.8591703355916499363e01, X 7.4321196680624245801e02, }; static double QJ1[] = { X 3.3307310774649071172e02, X -2.9458766545509337327e03, X 1.8680990008359188352e04, X -8.4055062591169562211e04, X 2.4599102262586308984e05, X -3.5783478026152301072e05, X -2.5258076240801555057e01, }; X /* X * Coefficients for rational approximation of X * (y0(x) - 2 LN(x/XY0) j0(x)) / (x**2 - XY0**2), X * XSMALL < |x| <= 3.0 X */ static double PY0[] = { X 1.0102532948020907590e04, X -2.1287548474401797963e06, X 2.0422274357376619816e08, X -8.3716255451260504098e09, X 1.0723538782003176831e11, X -1.8402381979244993524e01, }; static double QY0[] = { X 6.6475986689240190091e02, X 2.3889393209447253406e05, X 5.5662956624278251596e07, X 8.1617187777290363573e09, X 5.8873865738997033405e11, }; X /* X * Coefficients for rational approximation of X * (y0(x) - 2 LN(x/XY1) j0(x)) / (x**2 - XY1**2), X * 3.0 < |x| <= 5.5 X */ static double PY1[] = { X -1.4566865832663635920e04, X 4.6905288611678631510e06, X -6.9590439394619619534e08, X 4.3600098638603061642e10, X -5.5107435206722644429e11, X -2.2213976967566192242e13, X 1.7427031242901594547e01, }; static double QY1[] = { X 8.3030857612070288823e02, X 4.0669982352539552018e05, X 1.3960202770986831075e08, X 3.4015103849971240096e10, X 5.4266824419412347550e12, X 4.3386146580707264428e14, }; X /* X * Coefficients for rational approximation of X * (y0(x) - 2 LN(x/XY2) j0(x)) / (x**2 - XY2**2), X * 5.5 < |x| <= 8.0 X */ static double PY2[] = { X 2.1363534169313901632e04, X -1.0085539923498211426e07, X 2.1958827170518100757e09, X -1.9363051266772083678e11, X -1.2829912364088687306e11, X 6.7016641869173237784e14, X -8.0728726905150210443e15, X -1.7439661319197499338e01, }; static double QY2[] = { X 8.7903362168128450017e02, X 5.3924739209768057030e05, X 2.4727219475672302327e08, X 8.6926121104209825246e10, X 2.2598377924042897629e13, X 3.9272425569640309819e15, X 3.4563724628846457519e17, }; X /* X * Coefficients for Hart,s approximation, |x| > 8.0 X */ static double P0[] = { X 3.4806486443249270347e03, X 2.1170523380864944322e04, X 4.1345386639580765797e04, X 2.2779090197304684302e04, X 8.8961548424210455236e-01, X 1.5376201909008354296e02, }; static double Q0[] = { X 3.5028735138235608207e03, X 2.1215350561880115730e04, X 4.1370412495510416640e04, X 2.2779090197304684318e04, X 1.5711159858080893649e02, }; static double P1[] = { X -2.2300261666214198472e01, X -1.1183429920482737611e02, X -1.8591953644342993800e02, X -8.9226600200800094098e01, X -8.8033303048680751817e-03, X -1.2441026745835638459e00, }; static double Q1[] = { X 1.4887231232283756582e03, X 7.2642780169211018836e03, X 1.1951131543434613647e04, X 5.7105024128512061905e03, X 9.0593769594993125859e01, }; X static double #if defined(__STDC__) || defined(__GNUC__) caljy0(double x,int jint) #else caljy0(x,jint) double x; int jint; #endif { X int i; X double resj,down,up,xden,xnum,w,wsq,z,zsq; X X if (jint && x <= ZERO) /* Check for error conditions */ X return -XINF; #define ax x X else if ((ax = fabs(x)) > XMAX) X return ZERO; /* X * Calculate j0 or y0 for |x| > 8.0 X */ X if (ax > EIGHT) { X z = EIGHT/ax; X w = ax/TWOPI; X w = floor(w)+ONEOV8; X w = (ax-w*TWOPI1)-w*TWOPI2; X zsq = z*z; X xnum = P0[4]*zsq+P0[5]; X xden = zsq+Q0[4]; X up = P1[4]*zsq+P1[5]; X down = zsq+Q1[4]; X for (i = 0; i <= 3; i++) { X xnum = xnum*zsq+P0[i]; X xden = xden*zsq+Q0[i]; X up = up*zsq+P1[i]; X down = down*zsq+Q1[i]; X } #define r0 xnum #define r1 up X r0 = xnum/xden; X r1 = up/down; X return sqrt(PI2/ax)*(jint ? r0*sin(w)+z*r1*cos(w) : X r0*cos(w)-z*r1*sin(w)); #undef r1 #undef r0 X } X if (ax <= XSMALL) X return jint ? PI2*(log(ax)+CONS) : ONE; /* X * Calculate j0 for appropriate interval, preserving X * accuracy near the zero of j0 X */ X zsq = ax*ax; X if (ax <= FOUR) { X xnum = (PJ0[4]*zsq+PJ0[5])*zsq+PJ0[6]; X xden = zsq+QJ0[4]; X for (i = 0; i <= 3; i++) { X xnum = xnum*zsq+PJ0[i]; X xden = xden*zsq+QJ0[i]; X } #define prod resj X prod = ((ax-XJ01/TWO56)-XJ02)*(ax+XJ0); X } X else { X wsq = ONE-zsq/SIXTY4; X xnum = PJ1[6]*wsq+PJ1[7]; X xden = wsq+QJ1[6]; X for (i = 0; i <= 5; i++) { X xnum = xnum*wsq+PJ1[i]; X xden = xden*wsq+QJ1[i]; X } X prod = (ax+XJ1)*((ax-XJ11/TWO56)-XJ12); X } #define result resj X result = prod*xnum/xden; #undef prod X if (!jint) X return result; /* X * Calculate y0. First find resj = pi/2 ln(x/xn) j0(x), X * where xn is a zero of y0 X */ #define xy z X if (ax <= THREE) { X up = (ax-XY01/TWO56)-XY02; X xy = XY0; X } X else if (ax <= FIVE5) { X up = (ax-XY11/TWO56)-XY12; X xy = XY1; X } X else { X up = (ax-XY21/TWO56)-XY22; X xy = XY2; X } X down = ax+xy; X if (fabs(up) < P17*down) { X w = up/down; X wsq = w*w; X xnum = PLG[0]; X xden = wsq+QLG[0]; X for (i = 1; i <= 3; i++) { X xnum = xnum*wsq+PLG[i]; X xden = xden*wsq+QLG[i]; X } X resj = PI2*result*w*xnum/xden; X } X else X resj = PI2*result*log(ax/xy); #undef xy #undef result /* X * Now calculate y0 for appropriate interval, preserving X * accuracy near the zero of y0 X */ X if (ax <= THREE) { X xnum = PY0[5]*zsq+PY0[0]; X xden = zsq+QY0[0]; X for (i = 1; i <= 4; i++) { X xnum = xnum*zsq+PY0[i]; X xden = xden*zsq+QY0[i]; X } X } X else if (ax <= FIVE5) { #undef ax X xnum = PY1[6]*zsq+PY1[0]; X xden = zsq+QY1[0]; X for (i = 1; i <= 5; i++) { X xnum = xnum*zsq+PY1[i]; X xden = xden*zsq+QY1[i]; X } X } X else { X xnum = PY2[7]*zsq+PY2[0]; X xden = zsq+QY2[0]; X for (i = 1; i <= 6; i++) { X xnum = xnum*zsq+PY2[i]; X xden = xden*zsq+QY2[i]; X } X } X return resj+up*down*xnum/xden; } X /* X * This subprogram computes approximate values for Bessel functions X * of the first kind of order zero for arguments |x| <= XMAX X * (see comments heading caljy0). X */ float j0f(float x) { X return ((float)caljy0(x,0)); } X /* X * This subprogram computes approximate values for Bessel functions X * of the second kind of order zero for arguments 0 < x <= XMAX X * (see comments heading caljy0). X */ float y0f(float x) { X return ((float)caljy0(x,1)); } SHAR_EOF chmod 0644 j0f.c || echo 'restore of j0f.c failed' Wc_c="`wc -c < 'j0f.c'`" test 12458 -eq "$Wc_c" || echo 'j0f.c: original size 12458, current size' "$Wc_c" rm -f _shar_wnt_.tmp fi # ============= lgammaf.c ============== if test -f 'lgammaf.c' -a X"$1" != X"-c"; then echo 'x - skipping lgammaf.c (File already exists)' rm -f _shar_wnt_.tmp else > _shar_wnt_.tmp echo 'x - extracting lgammaf.c (Text)' sed 's/^X//' << 'SHAR_EOF' > 'lgammaf.c' && #ifndef lint static char sccsid[] = "@(#)lgamma.c 1.2 (ucb.beef) 3/1/89"; #endif /* !defined(lint) */ /* X * This routine calculates the log(GAMMA) function for a positive real X * argument x. Computation is based on an algorithm outlined in X * references 1 and 2. The program uses rational functions that X * theoretically approximate log(GAMMA) to at least 18 significant X * decimal digits. The approximation for x > 12 is from reference X * 3, while approximations for x < 12.0 are similar to those in X * reference 1, but are unpublished. The accuracy achieved depends X * on the arithmetic system, the compiler, the intrinsic functions, X * and proper selection of the machine-dependent constants. X * X ********************************************************************** X ********************************************************************** X * X * Explanation of machine-dependent constants X * X * beta - radix for the floating-point representation X * maxexp - the smallest positive power of beta that overflows X * XBIG - largest argument for which LN(GAMMA(x)) is representable X * in the machine, i.e., the solution to the equation X * LN(GAMMA(XBIG)) = beta**maxexp X * XINF - largest machine representable floating-point number; X * approximately beta**maxexp. X * EPS - The smallest positive floating-point number such that X * 1.0+EPS > 1.0 X * FRTBIG - Rough estimate of the fourth root of XBIG X * X * X * Approximate values for some important machines are: X * X * beta maxexp XBIG X * X * CRAY-1 (S.P.) 2 8191 9.62E+2461 X * Cyber 180/855 X * under NOS (S.P.) 2 1070 1.72E+319 X * IEEE (IBM/XT, X * SUN, etc.) (S.P.) 2 128 4.08E+36 X * IEEE (IBM/XT, X * SUN, etc.) (D.P.) 2 1024 2.55D+305 X * IBM 3033 (D.P.) 16 63 4.29D+73 X * VAX D-Format (D.P.) 2 127 2.05D+36 X * VAX G-Format (D.P.) 2 1023 1.28D+305 X * X * X * XINF EPS FRTBIG X * X * CRAY-1 (S.P.) 5.45E+2465 7.11E-15 3.13E+615 X * Cyber 180/855 X * under NOS (S.P.) 1.26E+322 3.55E-15 6.44E+79 X * IEEE (IBM/XT, X * SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.42E+9 X * IEEE (IBM/XT, X * SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.25D+76 X * IBM 3033 (D.P.) 7.23D+75 2.22D-16 2.56D+18 X * VAX D-Format (D.P.) 1.70D+38 1.39D-17 1.20D+9 X * VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.89D+76 X * X *************************************************************** X *************************************************************** X * X * Error returns X * X * The program returns the value XINF for x <= 0.0 or when X * overflow would occur. The computation is believed to X * be free of underflow and overflow. X * X * X * Intrinsic functions required are: X * X * log X * X * X * References: X * X * 1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for X * the Natural Logarithm of the Gamma Function,' Math. Comp. 21, X * 1967, pp. 198-203. X * X * 2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, X * 1969. X * X * 3) Hart, Et. Al., Computer Approximations, Wiley and sons, New X * York, 1968. X * X * X * Authors: W. J. Cody and L. Stoltz X * Argonne National Laboratory X * X * Latest modification: June 16, 1988 X */ X #include "fpumath.h" X X /* Machine-dependent constants */ #if defined(vax) || defined(tahoe) #define XBIG (double)2.05e36 #define XINF (double)1.70e38 #define EPS (double)1.39e-17 #define FRTBIG (double)1.20e9 #else /* defined(vax) || defined(tahoe) */ #define XBIG (double)2.55e305 #define XINF MAXFLOAT #define EPS (double)2.22e-16 #define FRTBIG (double)2.25e76 #endif /* defined(vax) || defined(tahoe) */ X /* Mathematical constants */ #define ONE (double)1 #define HALF (double)0.5 #define TWELVE (double)12 #define ZERO (double)0 #define FOUR (double)4 #define THRHAL (double)1.5 #define TWO (double)2 #define PNT68 (double)0.6796875 #define SQRTPI (double)0.9189385332046727417803297 X /* X * Numerator and denominator coefficients for rational minimax Approximation X * over (0.5,1.5). X */ static double D1 = -5.772156649015328605195174e-1; static double P1[] = { X 4.945235359296727046734888e0, X 2.018112620856775083915565e2, X 2.290838373831346393026739e3, X 1.131967205903380828685045e4, X 2.855724635671635335736389e4, X 3.848496228443793359990269e4, X 2.637748787624195437963534e4, X 7.225813979700288197698961e3, }; static double Q1[] = { X 6.748212550303777196073036e1, X 1.113332393857199323513008e3, X 7.738757056935398733233834e3, X 2.763987074403340708898585e4, X 5.499310206226157329794414e4, X 6.161122180066002127833352e4, X 3.635127591501940507276287e4, X 8.785536302431013170870835e3, }; X /* X * Numerator and denominator coefficients for rational minimax Approximation X * over (1.5,4.0). X */ static double D2 = 4.227843350984671393993777e-1; static double P2[] = { X 4.974607845568932035012064e0, X 5.424138599891070494101986e2, X 1.550693864978364947665077e4, X 1.847932904445632425417223e5, X 1.088204769468828767498470e6, X 3.338152967987029735917223e6, X 5.106661678927352456275255e6, X 3.074109054850539556250927e6, }; static double Q2[] = { X 1.830328399370592604055942e2, X 7.765049321445005871323047e3, X 1.331903827966074194402448e5, X 1.136705821321969608938755e6, X 5.267964117437946917577538e6, X 1.346701454311101692290052e7, X 1.782736530353274213975932e7, X 9.533095591844353613395747e6, }; X /* X * Numerator and denominator coefficients for rational minimax Approximation X * over (4.0,12.0). X */ static double D4 = 1.791759469228055000094023e0; static double P4[] = { X 1.474502166059939948905062e4, X 2.426813369486704502836312e6, X 1.214755574045093227939592e8, X 2.663432449630976949898078e9, X 2.940378956634553899906876e10, X 1.702665737765398868392998e11, X 4.926125793377430887588120e11, X 5.606251856223951465078242e11, }; static double Q4[] = { X 2.690530175870899333379843e3, X 6.393885654300092398984238e5, X 4.135599930241388052042842e7, X 1.120872109616147941376570e9, X 1.488613728678813811542398e10, X 1.016803586272438228077304e11, X 3.417476345507377132798597e11, X 4.463158187419713286462081e11, }; X /* X * Coefficients for minimax approximation over (12, INF). X */ static double C[] = { X -1.910444077728e-03, X 8.4171387781295e-04, X -5.952379913043012e-04, X 7.93650793500350248e-04, X -2.777777777777681622553e-03, X 8.333333333333333331554247e-02, X 5.7083835261e-03, }; X float lgammaf(float x) { X register i; X double res,corr,xden,xnum,dtmp; X #define y x X if (y > ZERO && y <= XBIG) { X if (y <= EPS) X res = -log(y); X else if (y <= THRHAL) { /* EPS < x <= 1.5 */ #define xm1 dtmp X if (y < PNT68) { X corr = -log(y); X xm1 = y; X } X else { X corr = ZERO; X xm1 = y-HALF; xm1 -= HALF; X } X if (y <= HALF || y >= PNT68) { X xden = ONE; X xnum = ZERO; X for (i = 0; i <= 7; i++) { X xnum = xnum*xm1+P1[i]; X xden = xden*xm1+Q1[i]; X } X res = xnum/xden; res = corr+xm1*(D1+xm1*res); #undef xm1 X } X else { #define xm2 dtmp X xm2 = y-HALF; xm2 -= HALF; X xden = ONE; X xnum = ZERO; X for (i = 0; i <= 7; i++) { X xnum = xnum*xm2+P2[i]; X xden = xden*xm2+Q2[i]; X } X res = xnum/xden; res = corr+xm2*(D2+xm2*res); X } X } X else if (y <= FOUR) { /* 1.5 < x <= 4.0 */ X xm2 = y-TWO; X xden = ONE; X xnum = ZERO; X for (i = 0; i <= 7; i++) { X xnum = xnum*xm2+P2[i]; X xden = xden*xm2+Q2[i]; X } X res = xnum/xden; res = xm2*(D2+xm2*res); #undef xm2 X } X else if (y <= TWELVE) { /* 4.0 < x <= 12.0 */ #define xm4 dtmp X xm4 = y-FOUR; X xden = -ONE; X xnum = ZERO; X for (i = 0; i <= 7; i++) { X xnum = xnum*xm4+P4[i]; X xden = xden*xm4+Q4[i]; X } X res = xnum/xden; res = D4+xm4*res; #undef xm4 X } X else { /* x >= 12.0 */ X res = ZERO; X if (y <= FRTBIG) { X res = C[6]; #define ysq dtmp X ysq = y*y; X for (i = 0; i <= 5; i++) X res = res/ysq+C[i]; #define ysq dtmp X } X res /= y; X corr = log(y); X res += SQRTPI; res -= HALF*corr; X res += y*(corr-ONE); X } #undef y X } X else /* Return for bad arguments */ X res = XINF; X return ((float)res); } SHAR_EOF chmod 0644 lgammaf.c || echo 'restore of lgammaf.c failed' Wc_c="`wc -c < 'lgammaf.c'`" test 8280 -eq "$Wc_c" || echo 'lgammaf.c: original size 8280, current size' "$Wc_c" rm -f _shar_wnt_.tmp fi # ============= gammaf.c ============== if test -f 'gammaf.c' -a X"$1" != X"-c"; then echo 'x - skipping gammaf.c (File already exists)' rm -f _shar_wnt_.tmp else > _shar_wnt_.tmp echo 'x - extracting gammaf.c (Text)' sed 's/^X//' << 'SHAR_EOF' > 'gammaf.c' && #ifndef lint static char sccsid[] = "@(#)gamma.c 1.6 (ucb.beef) 10/3/89"; #endif /* !defined(lint) */ /* X * This routine calculates the GAMMA function for a "double" argument X. X * Computation is based on an algorithm outlined in reference 1. X * The program uses rational functions that approximate the GAMMA X * function to at least 20 significant decimal digits. Coefficients X * for the approximation over the interval (1,2) are unpublished. X * Those for the approximation for X .GE. 12 are from reference 2. X * The accuracy achieved depends on the arithmetic system, the X * compiler, the intrinsic functions, and proper selection of the X * machine-dependent constants. X * X * X ********************************************************************* X ********************************************************************* X * X * Explanation of machine-dependent constants X * X * beta - radix for the floating-point representation X * maxexp - the smallest positive power of beta that overflows X * XBIG - the largest argument for which GAMMA(X) is representable X * in the machine, i.e., the solution to the equation X * GAMMA(XBIG) = beta**maxexp X * XINF - the largest machine representable floating-point number; X * approximately beta**maxexp X * EPS - the smallest positive floating-point number such that X * 1.0+EPS .GT. 1.0 X * XMININ - the smallest positive floating-point number such that X * 1/XMININ is machine representable X * X * Approximate values for some important machines are: X * X * beta maxexp XBIG X * X * CRAY-1 (S.P.) 2 8191 967.095 X * Cyber 180/855 X * under NOS (S.P.) 2 1070 177.980 X * IEEE (IBM/XT, X * SUN, etc.) (S.P.) 2 128 35.299 X * IEEE (IBM/XT, X * SUN, etc.) (D.P.) 2 1024 171.624 X * IBM 3033 (D.P.) 16 63 57.801 X * VAX D-Format (D.P.) 2 127 34.844 X * VAX G-Format (D.P.) 2 1023 171.668 X * X * XINF EPS XMININ X * X * CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 X * Cyber 180/855 X * under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 X * IEEE (IBM/XT, X * SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 X * IEEE (IBM/XT, X * SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 X * IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 X * VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39 X * VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.12D-308 X * X ********************************************************************* X ********************************************************************* X * X * Error returns X * X * The program returns the value XINF for singularities or X * when overflow would occur. The computation is believed X * to be free of underflow and overflow. X * X * X * Intrinsic functions required are: X * X * exp, floor, log, sin X * X * X * References: "An Overview of Software Development for Special X * Functions", W. J. Cody, Lecture Notes in Mathematics, X * 506, Numerical Analysis Dundee, 1975, G. A. Watson X * (ed.), Springer Verlag, Berlin, 1976. X * X * Computer Approximations, Hart, Et. Al., Wiley and X * sons, New York, 1968. X * X * Latest modification: May 30, 1989 X * X * Authors: W. J. Cody and L. Stoltz X * Applied Mathematics Division X * Argonne National Laboratory X * Argonne, IL 60439 X */ X #include "fpumath.h" X X /* Machine dependent parameters */ #if defined(vax) || defined(tahoe) #define XBIG (double)34.844e0 #define XMININ (double)5.88e-39 #define EPS (double)1.39e-17 #define XINF (double)1.70e38 #else /* defined(vax) || defined(tahoe) */ #define XBIG (double)171.624e0 #define XMININ (double)2.23e-308 #define EPS (double)2.22e-16 #define XINF MAXFLOAT #endif /* defined(vax) || defined(tahoe) */ X /* Mathematical constants */ #define ONE (double)1 #define HALF (double)0.5 #define TWELVE (double)12 #define ZERO (double)0 #define SQRTPI (double)0.9189385332046727417803297e0 #define PI (double)3.1415926535897932384626434e0 X typedef int logical; #define FALSE (logical)0 #define TRUE (logical)1 X /* X * Numerator and denominator coefficients for rational minimax X * approximation over (1,2). X */ static double P[] = { X -1.71618513886549492533811e0, X 2.47656508055759199108314e1, X -3.79804256470945635097577e2, X 6.29331155312818442661052e2, X 8.66966202790413211295064e2, X -3.14512729688483675254357e4, X -3.61444134186911729807069e4, X 6.64561438202405440627855e4, }; static double Q[] = { X -3.08402300119738975254353e1, X 3.15350626979604161529144e2, X -1.01515636749021914166146e3, X -3.10777167157231109440444e3, X 2.25381184209801510330112e4, X 4.75584627752788110767815e3, X -1.34659959864969306392456e5, X -1.15132259675553483497211e5, }; X /* X * Coefficients for minimax approximation over (12, INF). X */ static double C[] = { X -1.910444077728e-03, X 8.4171387781295e-04, X -5.952379913043012e-04, X 7.93650793500350248e-04, X -2.777777777777681622553e-03, X 8.333333333333333331554247e-02, X 5.7083835261e-03, }; X float gammaf(float x) { X register i,n; X logical parity; X double fact,res,xden,xnum,dtmp1,dtmp2; X X parity = FALSE; X fact = ONE; X n = 0; #define y x X if (y <= ZERO) { /* Arg. is negative */ X y = -x; #define y1 dtmp1 X y1 = floor(y); X res = y-y1; X if (res != ZERO) { X parity = floor(y1/2)*2 != y1; X fact = -PI/sin(PI*res); X y += ONE; X } X else X return XINF; X } X /* Arg. is positive */ X if (y < EPS) { /* Arg. < EPS */ X if (y >= XMININ) X res = ONE/y; X else X return XINF; X } X else if (y < TWELVE) { #define y1 dtmp1 #define z dtmp2 X y1 = y; X if (y < ONE) { /* 0.0 < arg. < 1.0 */ X z = y; X y += ONE; X } X else { /* 1.0 < arg. < 12.0 */ X X n = (int)y; n--; X y -= (double)n; /* reduce arg. if needed */ X z = y-ONE; X } X /* Evaluate approximation for 1.0 < arg. < 2.0 */ X xnum = ZERO; X xden = ONE; X for (i = 0; i <= 7; i++) { X xnum = (xnum+P[i])*z; X xden = xden*z+Q[i]; X } X res = xnum/xden+ONE; X if (y1 < y) /* Adjust res for 0.0 < arg. < 1.0 */ X res /= y1; X else if (y1 > y) { /* Adjust res for 2.0 < arg. < 12.0 */ X for (i = 0; i <= n-1; i++) { X res *= y; X y += ONE; X } X } #undef z #undef y1 X } X else { /* Evaluate for arg. >= 12.0 */ X if (y <= XBIG) { #define ysq dtmp1 #define sum dtmp2 X ysq = y*y; X sum = C[6]; X for (i = 0; i <= 5; i++) X sum = sum/ysq+C[i]; X sum = sum/y-y; sum += SQRTPI; X sum += (y-HALF)*log(y); X res = exp(sum); #undef sum #undef ysq X } X else X return XINF; X } X /* Final adjustments and return */ X if (parity == TRUE) X res = -res; X if (fact != ONE) X res = fact/res; X return ((float)res); #undef y } SHAR_EOF chmod 0644 gammaf.c || echo 'restore of gammaf.c failed' Wc_c="`wc -c < 'gammaf.c'`" test 6996 -eq "$Wc_c" || echo 'gammaf.c: original size 6996, current size' "$Wc_c" rm -f _shar_wnt_.tmp fi # ============= j1f.c ============== if test -f 'j1f.c' -a X"$1" != X"-c"; then echo 'x - skipping j1f.c (File already exists)' rm -f _shar_wnt_.tmp else > _shar_wnt_.tmp echo 'x - extracting j1f.c (Text)' sed 's/^X//' << 'SHAR_EOF' > 'j1f.c' && #ifndef lint static char sccsid[] = "@(#)j1.c 1.1 (ucb.beef) 10/1/89"; #endif /* !defined(lint) */ /* X * This packet computes first-order Bessel functions of the first and X * second kind (j1 and y1), for real arguments x, where 0 < x <= XMAX X * for y1, and |x| <= XMAX for j1. It contains three function-type X * subprograms, j1, y1 and caljy1. The calling statements for the X * primary entries are: X * X * y = j1(x) X * and X * y = y1(x), X * X * where the entry points correspond to the functions j1(x) and y1(x), X * respectively. The routine caljy1() is intended for internal packet X * use only, all computations within the packet being concentrated in X * this one routine. The function subprograms invoke caljy1 with X * the statement X * X * result = caljy1(x,jint); X * X * where the parameter usage is as follows: X * X * Function Parameters for caljy1 X * call x result jint X * X * j1(x) |x| <= XMAX j1(x) 0 X * y1(x) 0 < x <= XMAX y1(x) 1 X * X * The main computation uses unpublished minimax rational X * approximations for x <= 8.0, and an approximation from the X * book Computer Approximations by Hart, et. al., Wiley and Sons, X * New York, 1968, for arguments larger than 8.0 Part of this X * transportable packet is patterned after the machine-dependent X * FUNPACK program for j1(x), but cannot match that version for X * efficiency or accuracy. This version uses rational functions X * that are theoretically accurate to at least 18 significant decimal X * digits for x <= 8, and at least 18 decimal places for x > 8. The X * accuracy achieved depends on the arithmetic system, the compiler, X * the intrinsic functions, and proper selection of the machine- X * dependent constants. X * X ********************************************************************* X * X * Explanation of machine-dependent constants X * X * XINF = largest positive machine number X * XMAX = largest acceptable argument. The functions floor(), sin() X * and cos() must perform properly for fabs(x) <= XMAX. X * We recommend that XMAX be a small integer multiple of X * sqrt(1/eps), where eps is the smallest positive number X * such that 1+eps > 1. X * XSMALL = positive argument such that 1.0-(1/2)(x/2)**2 = 1.0 X * to machine precision for all fabs(x) <= XSMALL. X * We recommend that XSMALL < sqrt(eps)/beta, where beta X * is the floating-point radix (usually 2 or 16). X * X * Approximate values for some important machines are X * X * eps XMAX XSMALL XINF X * X * CDC 7600 (S.P.) 7.11E-15 1.34E+08 2.98E-08 1.26E+322 X * CRAY-1 (S.P.) 7.11E-15 1.34E+08 2.98E-08 5.45E+2465 X * IBM PC (8087) (S.P.) 5.96E-08 8.19E+03 1.22E-04 3.40E+38 X * IBM PC (8087) (D.P.) 1.11e-16 2.68e08 3.72e-09 1.79e308 X * IBM 195 (D.P.) 2.22e-16 6.87e09 9.09e-13 7.23e75 X * UNIVAC 1108 (D.P.) 1.73e-18 4.30e09 2.33e-10 8.98e307 X * VAX 11/780 (D.P.) 1.39e-17 1.07e09 9.31e-10 1.70e38 X * X ********************************************************************* X ********************************************************************* X * X * Error Returns X * X * The program returns the value zero for x > XMAX, and returns X * -XINF when BESLY1 is called with a negative or zero argument. X * X * X * Intrinsic functions required are: X * X * cos, fabs, floor, log, sin, sqrt X * X * X * Author: w. J. Cody X * Mathematics and Computer Science Division X * Argonne National Laboratory X * Argonne, IL 60439 X * X * Latest modification: November 10, 1987 X */ X #include "fpumath.h" X X /* Machine-dependent constants */ #if defined(vax) || defined(tahoe) #define XMAX (double)1.07e9 #define XSMALL (double)9.31e-10 #define XINF (double)1.7e38 #else /* defined(vax) || defined(tahoe) */ #define XMAX (double)2.68e08 #define XSMALL (double)3.72e-09 #define XINF MAXFLOAT #endif /* defined(vax) || defined(tahoe) */ X /* Mathematical constants */ #define EIGHT (double)8 #define FOUR (double)4 #define HALF (double)0.5 #define THROV8 (double)0.375 #define PI2 (double)6.3661977236758134308e-1 #define P17 (double)1.716e-1 #define TWOPI (double)6.2831853071795864769e0 #define ZERO (double)0 #define TWOPI1 (double)6.28125 #define TWOPI2 (double)1.9353071795864769253e-03 #define TWO56 (double)256 #define RTPI2 (double)7.9788456080286535588e-1 X /* Zeroes of Bessel functions */ #define XJ0 (double)3.8317059702075123156e0 #define XJ1 (double)7.0155866698156187535e0 #define XY0 (double)2.1971413260310170351e0 #define XY1 (double)5.4296810407941351328e0 #define XJ01 (double)981 #define XJ02 (double)(-3.2527979248768438556e-4) #define XJ11 (double)1796 #define XJ12 (double)(-3.8330184381246462950e-5) #define XY01 (double)562 #define XY02 (double)1.8288260310170351490e-3 #define XY11 (double)1390 #define XY12 (double)(-6.4592058648672279948e-6) X /* X * Coefficients for rational approximation to ln(x/a) X */ static double PLG[] = { X -2.4562334077563243311e01, X 2.3642701335621505212e02, X -5.4989956895857911039e02, X 3.5687548468071500413e02, }; static double QLG[] = { X -3.5553900764052419184e01, X 1.9400230218539473193e02, X -3.3442903192607538956e02, X 1.7843774234035750207e02, }; X /* X * Coefficients for rational approximation of X * j1(x) / (x * (x**2 - XJ0**2)), XSMALL < |x| <= 4.0 X */ static double PJ0[] = { X 9.8062904098958257677e05, X -1.1548696764841276794e08, X 6.6781041261492395835e09, X -1.4258509801366645672e11, X -4.4615792982775076130e03, X 1.0650724020080236441e01, X -1.0767857011487300348e-02, }; static double QJ0[] = { X 5.9117614494174794095e05, X 2.0228375140097033958e08, X 4.2091902282580133541e10, X 4.1868604460820175290e12, X 1.0742272239517380498e03, }; X /* X * Coefficients for rational approximation of X * j1(x) / (x * (x**2 - XJ1**2)), 4.0 < |x| <= 8.0 X */ static double PJ1[] = { X 4.6179191852758252280e00, X -7.1329006872560947377e03, X 4.5039658105749078904e06, X -1.4437717718363239107e09, X 2.3569285397217157313e11, X -1.6324168293282543629e13, X 1.1357022719979468624e14, X 1.0051899717115285432e15, }; static double QJ1[] = { X 1.1267125065029138050e06, X 6.4872502899596389593e08, X 2.7622777286244082666e11, X 8.4899346165481429307e13, X 1.7128800897135812012e16, X 1.7253905888447681194e18, X 1.3886978985861357615e03, }; X /* X * Coefficients for rational approximation of X * (y1(x) - 2 LN(x/XY0) j1(x)) / (x**2 - XY0**2), X * XSMALL < |x| <= 4.0 X */ static double PY0[] = { X 2.2157953222280260820e05, X -5.9157479997408395984e07, X 7.2144548214502560419e09, X -3.7595974497819597599e11, X 5.4708611716525426053e12, X 4.0535726612579544093e13, X -3.1714424660046133456e02, }; static double QY0[] = { X 8.2079908168393867438e02, X 3.8136470753052572164e05, X 1.2250435122182963220e08, X 2.7800352738690585613e10, X 4.1272286200406461981e12, X 3.0737873921079286084e14, }; X /* X * Coefficients for rational approximation of X * (y1(x) - 2 LN(x/XY1) j1(x)) / (x**2 - XY1**2), X * .0 < |x| <= 8.0 X */ static double PY1[] = { X 1.9153806858264202986e06, X -1.1957961912070617006e09, X 3.7453673962438488783e11, X -5.9530713129741981618e13, X 4.0686275289804744814e15, X -2.3638408497043134724e16, X -5.6808094574724204577e18, X 1.1514276357909013326e19, X -1.2337180442012953128e03, }; static double QY1[] = { X 1.2855164849321609336e03, X 1.0453748201934079734e06, X 6.3550318087088919566e08, X 3.0221766852960403645e11, X 1.1187010065856971027e14, X 3.0837179548112881950e16, X 5.6968198822857178911e18, X 5.3321844313316185697e20, }; X /* X * Coefficients for Hart,s approximation, |x| > 8.0 X */ static double P0[] = { X -1.0982405543459346727e05, X -1.5235293511811373833e06, X -6.6033732483649391093e06, X -9.9422465050776411957e06, X -4.4357578167941278571e06, X -1.6116166443246101165e03, }; static double Q0[] = { X -1.0726385991103820119e05, X -1.5118095066341608816e06, X -6.5853394797230870728e06, X -9.9341243899345856590e06, X -4.4357578167941278568e06, X -1.4550094401904961825e03, }; static double P1[] = { X 1.7063754290207680021e03, X 1.8494262873223866797e04, X 6.6178836581270835179e04, X 8.5145160675335701966e04, X 3.3220913409857223519e04, X 3.5265133846636032186e01, }; static double Q1[] = { X 3.7890229745772202641e04, X 4.0029443582266975117e05, X 1.4194606696037208929e06, X 1.8194580422439972989e06, X 7.0871281941028743574e05, X 8.6383677696049909675e02, }; X static double #if defined(__STDC__) || defined(__GNUC__) caljy1(double x, int jint) #else caljy1(x,jint) double x; int jint; #endif { X int i; X double ax,resj,down,up,xden,xnum,w,wsq,z,zsq; X X ax = fabs(x); /* Check for error conditions */ X if (jint && (x <= ZERO || x < HALF && ax*XINF < PI2)) X return -XINF; X else if (ax > XMAX) X return ZERO; /* X * Calculate j1 or y1 for |x| > 8.0 X */ X if (ax > EIGHT) { X z = EIGHT/ax; X w = floor(ax/TWOPI)+THROV8; X w = (ax-w*TWOPI1)-w*TWOPI2; X zsq = z*z; X xnum = P0[5]; X xden = zsq+Q0[5]; X up = P1[5]; X down = zsq+Q1[5]; X for (i = 0; i <= 4; i++) { X xnum = xnum*zsq+P0[i]; X xden = xden*zsq+Q0[i]; X up = up*zsq+P1[i]; X down = down*zsq+Q1[i]; X } #define r0 xnum #define r1 up X r0 = xnum/xden; X r1 = up/down; X return RTPI2/sqrt(ax)*(jint ? r0*sin(w)+z*r1*cos(w) : X (x < ZERO ? z*r1*sin(w)-r0*cos(w) : X r0*cos(w)-z*r1*sin(w))); #undef r1 #undef r0 X } X else if (ax <= XSMALL) X return jint ? -PI2/ax : x*HALF; /* X * Calculate j1 for appropriate interval, preserving X * accuracy near the zero of j1 X */ X zsq = ax*ax; X if (ax <= FOUR) { X xnum = (PJ0[6]*zsq+PJ0[5])*zsq+PJ0[4]; X xden = zsq+QJ0[4]; X for (i = 0; i <= 3; i++) { X xnum = xnum*zsq+PJ0[i]; X xden = xden*zsq+QJ0[i]; X } #define prod resj X prod = x*((ax-XJ01/TWO56)-XJ02)*(ax+XJ0); X } X else { X xnum = PJ1[0]; X xden = (zsq+QJ1[6])*zsq+QJ1[0]; X for (i = 1; i <= 5; i++) { X xnum = xnum*zsq+PJ1[i]; X xden = xden*zsq+QJ1[i]; X } X xnum = xnum*(ax-EIGHT)*(ax+EIGHT)+PJ1[6]; X xnum = xnum*(ax-FOUR)*(ax+FOUR)+PJ1[7]; X prod = x*((ax-XJ11/TWO56)-XJ12)*(ax+XJ1); X } #define result resj X result = prod*(xnum/xden); #undef prod X if (!jint) X return result; /* X * Calculate y1. First find resj = pi/2 ln(x/xn) j1(x), X * where xn is a zero of y1 X */ #define xy z X if (ax <= FOUR) { X up = (ax-XY01/TWO56)-XY02; X xy = XY0; X } X else { X up = (ax-XY11/TWO56)-XY12; X xy = XY1; X } X down = ax+xy; X if (fabs(up) < P17*down) { X w = up/down; X wsq = w*w; X xnum = PLG[0]; X xden = wsq+QLG[0]; X for (i = 1; i <= 3; i++) { X xnum = xnum*wsq+PLG[i]; X xden = xden*wsq+QLG[i]; X } X resj = PI2*result*w*xnum/xden; X } X else X resj = PI2*result*log(ax/xy); #undef xy #undef result /* X * Now calculate y1 for appropriate interval, preserving X * accuracy near the zero of y1 X */ X if (ax <= FOUR) { X xnum = PY0[6]*zsq+PY0[0]; X xden = zsq+QY0[0]; X for (i = 1; i <= 5; i++) { X xnum = xnum*zsq+PY0[i]; X xden = xden*zsq+QY0[i]; X } X } X else { X xnum = PY1[8]*zsq+PY1[0]; X xden = zsq+QY1[0]; X for (i = 1; i <= 7; i++) { X xnum = xnum*zsq+PY1[i]; X xden = xden*zsq+QY1[i]; X } X } X up *= down; up /= ax; up *= xnum; up /= xden; up += resj; X return up; } X /* X * This subprogram computes approximate values for Bessel functions X * of the first kind of order zero for arguments |x| <= XMAX X * (see comments heading caljy1). X */ float j1f(float x) { X return ((float)caljy1(x,0)); } X /* X * This subprogram computes approximate values for Bessel functions X * of the second kind of order zero for arguments 0 < x <= XMAX X * (see comments heading caljy1). X */ float y1f(float x) { X return ((float)caljy1(x,1)); } SHAR_EOF chmod 0644 j1f.c || echo 'restore of j1f.c failed' Wc_c="`wc -c < 'j1f.c'`" test 11769 -eq "$Wc_c" || echo 'j1f.c: original size 11769, current size' "$Wc_c" rm -f _shar_wnt_.tmp fi # ============= erff.c ============== if test -f 'erff.c' -a X"$1" != X"-c"; then echo 'x - skipping erff.c (File already exists)' rm -f _shar_wnt_.tmp else > _shar_wnt_.tmp echo 'x - extracting erff.c (Text)' sed 's/^X//' << 'SHAR_EOF' > 'erff.c' && #ifndef lint static char sccsid[] = "@(#)erf.c 1.3 (ucb.beef) 10/2/89"; #endif /* !defined(lint) */ /* X * This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) X * for a "float" argument x. It contains three subprograms: X * erff(), erfcf(), and erfcxf() and * one "static" subprogram, X * calerf. The calling statements for the primary entries are: X * X * y = erf(x), X * X * y = erfc(x), X * and X * y = erfcx(x). X * X * The routine calerf() is intended for internal packet use only, X * all computations within the packet being concentrated in this X * routine. The 3 primary subprograms invoke calerf with the X * statement X * X * result = calerf(arg,jint); X * X * where the parameter usage is as follows X * X * Function Parameters for calerf X * call arg result jint X * X * erf(arg) ANY "double" ARGUMENT erf(arg) 0 X * erfc(arg) fabs(arg) < XBIG erfc(arg) 1 X * erfcx(arg) XNEG < arg < XMAX erfcx(arg) 2 X * X * The main computation evaluates near-minimax approximations X * from "Rational Chebyshev approximations for the error function" X * by W. J. Cody, Math. Comp., 1969, PP. 631-638. This X * transportable program uses rational functions that theoretically X * approximate erf(x) and erfc(x) to at least 18 significant X * decimal digits. The accuracy achieved depends on the arithmetic X * system, the compiler, the intrinsic functions, and proper X * selection of the machine-dependent constants. X * X ******************************************************************** X ******************************************************************** X * X * Explanation of machine-dependent constants X * X * XMIN = the smallest positive floating-point number. X * XINF = the largest positive finite floating-point number. X * XNEG = the largest negative argument acceptable to erfcx; X * the negative of the solution to the equation X * 2*exp(x*x) = XINF. X * XSMALL = argument below which erf(x) may be represented by X * 2*x/sqrt(pi) and above which x*x will not underflow. X * A conservative value is the largest machine number x X * such that 1.0 + x = 1.0 to machine precision. X * XBIG = largest argument acceptable to erfc; solution to X * the equation: W(x) * (1-0.5/x**2) = XMIN, where X * W(x) = exp(-x*x)/[x*sqrt(pi)]. X * XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to X * machine precision. A conservative value is X * 1/[2*sqrt(XSMALL)] X * XMAX = largest acceptable argument to erfcx; the minimum X * of XINF and 1/[sqrt(pi)*XMIN]. X * X * Approximate values for some important machines are: X * X * XMIN XINF XNEG XSMALL X * X * CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 X * CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 X * IEEE (IBM/XT, X * SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 X * IEEE (IBM/XT, X * SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 X * IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 X * UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 X * VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 X * VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 X * X * X * XBIG XHUGE XMAX X * X * CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293 X * CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 X * IEEE (IBM/XT, X * SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 X * IEEE (IBM/XT, X * SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 X * IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 X * UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 X * VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 X * VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 X * X ******************************************************************** X ******************************************************************** X * X * Error returns X * X * The program returns erfc = 0 for arg .GE. XBIG; X * X * erfcx = XINF for arg .LT. XNEG; X * and X * erfcx = 0 for arg .GE. XMAX. X * X * X * Intrinsic functions required are: X * X * fabs, exp X * X * X * Author: W. J. Cody X * Mathematics and Computer Science Division X * Argonne National Laboratory X * Argonne, IL 60439 X * X * Latest modification: June 16, 1988 X */ X #include "fpumath.h" X X /* Machine-dependent constants */ #if defined(vax) || defined(tahoe) #define XINF (double)1.70e38 #define XNEG (double)(-9.345e0) #define XSMALL (double)1.39e-17 #define XBIG (double)9.269e0 #define XHUGE (double)1.90e8 #define XMAX (double)1.70e38 #else /* defined(vax) || defined(tahoe) */ #define XINF MAXFLOAT #define XNEG (double)(-26.628e0) #define XSMALL (double)1.11e-16 #define XBIG (double)26.543e0 #define XHUGE (double)6.71e7 #define XMAX (double)2.53e307 #endif /* defined(vax) || defined(tahoe) */ X /* Mathematical constants */ #define FOUR (double)4 #define ONE (double)1 #define HALF (double)0.5 #define TWO (double)2 #define ZERO (double)0 #define SQRPI (double)5.6418958354775628695e-1 #define THRESH (double)0.46875 #define SIXTEN (double)16 X /* X * Coefficients for approximation to erf in first interval X */ static double A[] = { X 3.16112374387056560e00, X 1.13864154151050156e02, X 3.77485237685302021e02, X 3.20937758913846947e03, X 1.85777706184603153e-1, }; static double B[] = { X 2.36012909523441209e01, X 2.44024637934444173e02, X 1.28261652607737228e03, X 2.84423683343917062e03, }; X /* X * Coefficients for approximation to erfc in second interval X */ static double C[] = { X 5.64188496988670089e-1, X 8.88314979438837594e0, X 6.61191906371416295e01, X 2.98635138197400131e02, X 8.81952221241769090e02, X 1.71204761263407058e03, X 2.05107837782607147e03, X 1.23033935479799725e03, X 2.15311535474403846e-8, }; static double D[] = { X 1.57449261107098347e01, X 1.17693950891312499e02, X 5.37181101862009858e02, X 1.62138957456669019e03, X 3.29079923573345963e03, X 4.36261909014324716e03, X 3.43936767414372164e03, X 1.23033935480374942e03, }; X /* X * Coefficients for approximation to erfc in third interval X */ static double P[] = { X 3.05326634961232344e-1, X 3.60344899949804439e-1, X 1.25781726111229246e-1, X 1.60837851487422766e-2, X 6.58749161529837803e-4, X 1.63153871373020978e-2, }; static double Q[] = { X 2.56852019228982242e00, X 1.87295284992346047e00, X 5.27905102951428412e-1, X 6.05183413124413191e-2, X 2.33520497626869185e-3, }; X static double #if defined(__STDC__) || defined(__GNUC__) calerf(double x,int jint) #else calerf(x,jint) double x; int jint; #endif { X register i,skip; X double y,ysq,xnum,xden,result; X X y = fabs(x); X if (y <= THRESH) { /* Evaluate erf for |x| <= 0.46875 */ X ysq = y > XSMALL ? y*y : ZERO; X xnum = A[4]*ysq; X xden = ysq; X for (i = 0; i <= 2; i++) { X xnum = (xnum+A[i])*ysq; X xden = (xden+B[i])*ysq; X } X result = x*(xnum+A[3])/(xden+B[3]); X if (jint) X result = ONE-result; X if (jint == 2) X result *= exp(ysq); X return result; X } X else if (y <= FOUR) { /* Evaluate erfc for 0.46875 <= |x| <= 4.0 */ X ysq = y*y; X xnum = C[8]*y; X xden = y; X for (i = 0; i <= 6; i++) { X xnum = (xnum+C[i])*y; X xden = (xden+D[i])*y; X } X result = (xnum+C[7])/(xden+D[7]); X if (jint != 2) { X i = (int)(y*SIXTEN); ysq = (double)i/SIXTEN; X result *= exp(-ysq*ysq)*exp(-(y-ysq)*(y+ysq)); X } X } X else { /* Evaluate erfc for |x| > 4.0 */ X result = ZERO; X skip = 0; X if (y >= XBIG) { X if (jint != 2 || y >= XMAX) X skip++; X else if (y >= XHUGE) { X result = SQRPI/y; SHAR_EOF true || echo 'restore of erff.c failed' fi echo 'End of mathlib2.0 part 2' echo 'File erff.c is continued in part 3' echo 3 > _shar_seq_.tmp exit 0 -- Glenn Geers | "So when it's over, we're back to people. Department of Theoretical Physics | Just to prove that human touch can have The University of Sydney | no equal." Sydney NSW 2006 Australia | - Basia Trzetrzelewska, 'Prime Time TV'