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