80386 alternative math lib v2.0 part02/06
Glenn Geers
glenn at qed.physics.su.OZ.AU
Mon Dec 17 08:01:28 AEST 1990
Submitted-by: root at 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'
More information about the Alt.sources
mailing list