#! /bin/sh
echo ########
echo # editorial remark, 11 Nov 1995, ehg@research.bell-labs.com
echo # See also http://netlib.bell-labs.com/netlib/toms/322.Z.
echo # Matthew Belmonte <mkb4@Cornell.edu> points out that
echo # line 114 here: p *= zk + w * z * (zk - 1.0)/(z-1.0);
echo # is wrong, and therefore unfairly disparages Tolman's improvement.
echo # Gary Perlman <perlman@cis.ohio-state.edu> agrees, but says even
echo # the slow version is plenty fast enough, and since this version
echo # works correctly, he is inclined to leave well enough alone.
echo ########
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
#	z.c
#	chisq.c
#	f.c
# This archive created: Wed May 27 02:19:35 1987
# By:	Gary Perlman (Wang Institute, Tyngsboro, MA 01879 USA)
export PATH; PATH=/bin:/usr/bin:$PATH
echo shar: "extracting 'z.c'" '(3123 characters)'
if test -f 'z.c'
then
	echo shar: "will not over-write existing file 'z.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'z.c'
	X/*HEADER
	X	Module:       z.c
	X	Purpose:      compute approximations to normal z distribution probabilities
	X	Programmer:   Gary Perlman
	X	Organization: Wang Institute, Tyngsboro, MA 01879
	X	Tester:       compile with -DZTEST to include main program
	X	Copyright:    none
	X	Tabstops:     4
	X*/
	X
	X/*LINTLIBRARY*/
	Xstatic char sccsfid[] = "@(#) z.c 5.1 (|stat) 12/26/85";
	X#include	<math.h>
	X
	X#define	Z_EPSILON      0.000001       /* accuracy of critz approximation */
	X#define	Z_MAX          6.0            /* maximum meaningful z value */
	X
	Xdouble	poz ();
	Xdouble	critz ();
	X
	X#ifdef	ZTEST
	Xmain ()
	X	{
	X	double	z;
	X	printf ("%4s  %10s  %10s  %10s\n",
	X		"z", "poz(z)", "poz(-z)", "z'");
	X	for (z = 0.0; z <= Z_MAX; z += .01)
	X		{
	X		printf ("%4.2f  %10.6f  %10.6f  %10.6f\n",
	X			z, poz (z), poz (-z), critz (poz (z)));
	X		}
	X	}
	X#endif	ZTEST
	X
	X/*FUNCTION poz: probability of normal z value */
	X/*ALGORITHM
	X	Adapted from a polynomial approximation in:
	X		Ibbetson D, Algorithm 209
	X		Collected Algorithms of the CACM 1963 p. 616
	X	Note:
	X		This routine has six digit accuracy, so it is only useful for absolute
	X		z values < 6.  For z values >= to 6.0, poz() returns 0.0.
	X*/
	Xdouble            /*VAR returns cumulative probability from -oo to z */
	Xpoz (z)
	Xdouble	z;        /*VAR normal z value */
	X	{
	X	double	y, x, w;
	X	
	X	if (z == 0.0)
	X		x = 0.0;
	X	else
	X		{
	X		y = 0.5 * fabs (z);
	X		if (y >= (Z_MAX * 0.5))
	X			x = 1.0;
	X		else if (y < 1.0)
	X			{
	X			w = y*y;
	X			x = ((((((((0.000124818987 * w
	X				-0.001075204047) * w +0.005198775019) * w
	X				-0.019198292004) * w +0.059054035642) * w
	X				-0.151968751364) * w +0.319152932694) * w
	X				-0.531923007300) * w +0.797884560593) * y * 2.0;
	X			}
	X		else
	X			{
	X			y -= 2.0;
	X			x = (((((((((((((-0.000045255659 * y
	X				+0.000152529290) * y -0.000019538132) * y
	X				-0.000676904986) * y +0.001390604284) * y
	X				-0.000794620820) * y -0.002034254874) * y
	X				+0.006549791214) * y -0.010557625006) * y
	X				+0.011630447319) * y -0.009279453341) * y
	X				+0.005353579108) * y -0.002141268741) * y
	X				+0.000535310849) * y +0.999936657524;
	X			}
	X		}
	X	return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5));
	X	}
	X
	X/*FUNCTION critz: compute critical z value to produce given probability */
	X/*ALGORITHM
	X	Begin with upper and lower limits for z values (maxz and minz)
	X	set to extremes.  Choose a z value (zval) between the extremes.
	X	Compute the probability of the z value.  Set minz or maxz, based
	X	on whether the probability is less than or greater than the
	X	desired p.  Continue adjusting the extremes until they are
	X	within Z_EPSILON of each other.
	X*/
	Xdouble        /*VAR returns z such that fabs (poz(p) - z) <= .000001 */
	Xcritz (p)
	Xdouble	p;    /*VAR critical probability level */
	X	{
	X	double	minz = -Z_MAX;    /* minimum of range of z */
	X	double	maxz = Z_MAX;     /* maximum of range of z */
	X	double	zval = 0.0;       /* computed/returned z value */
	X	double	poz (), pval;     /* prob (z) function, pval := poz (zval) */
	X	
	X	if (p <= 0.0 || p >= 1.0)
	X		return (0.0);
	X	
	X	while (maxz - minz > Z_EPSILON)
	X		{
	X		pval = poz (zval);
	X		if (pval > p)
	X			maxz = zval;
	X		else
	X			minz = zval;
	X		zval = (maxz + minz) * 0.5;
	X		}
	X	return (zval);
	X	}
SHAR_EOF
if test 3123 -ne "`wc -c < 'z.c'`"
then
	echo shar: "error transmitting 'z.c'" '(should have been 3123 characters)'
fi
fi
echo shar: "extracting 'chisq.c'" '(3003 characters)'
if test -f 'chisq.c'
then
	echo shar: "will not over-write existing file 'chisq.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'chisq.c'
	X/*
	X	Module:       chisq.c
	X	Purpose:      compute approximations to chisquare distribution probabilities
	X	Contents:     pochisq(), critchi()
	X	Uses:         poz() in z.c (Algorithm 209)
	X	Programmer:   Gary Perlman
	X	Organization: Wang Institute, Tyngsboro, MA 01879
	X	Tester:       compile with -DCHISQTEST to include main program
	X	Copyright:    none
	X	Tabstops:     4
	X*/
	X
	X/*LINTLIBRARY*/
	X#include <math.h>
	Xstatic char sccsfid[] = "@(#) chisq.c 5.2 (|stat) 12/08/86";
	X
	X#define	CHI_EPSILON     0.000001    /* accuracy of critchi approximation */
	X#define	CHI_MAX     99999.0         /* maximum chi square value */
	X
	X#define	LOG_SQRT_PI     0.5723649429247000870717135 /* log (sqrt (pi)) */
	X#define	I_SQRT_PI       0.5641895835477562869480795 /* 1 / sqrt (pi) */
	X#define	BIGX           20.0         /* max value to represent exp (x) */
	X#define	ex(x)             (((x) < -BIGX) ? 0.0 : exp (x))
	X
	Xdouble	pochisq ();
	Xdouble	critchi ();
	X
	X#ifdef	CHISQTEST
	Xdouble	Prob[] = { .10, .05, .01, .005, .001, -1.0 };
	Xmain ()
	X	{
	X	int 	df;
	X	int 	p;
	X	printf ("%-4s ", "df");
	X	for (p = 0; Prob[p] > 0.0; p++)
	X		printf ("%8.3f ", Prob[p]);
	X	putchar ('\n');
	X	for (df = 1; df < 30; df++)
	X		{
	X		printf ("%4d ", df);
	X		for (p = 0; Prob[p] > 0.0; p++)
	X			printf ("%8.3f ", critchi (Prob[p], df));
	X		putchar ('\n');
	X		}
	X	}
	X#endif	CHISQTEST
	X
	X/*FUNCTION pochisq: probability of chi sqaure value */
	X/*ALGORITHM Compute probability of chi square value.
	X	Adapted from:
	X		Hill, I. D. and Pike, M. C.  Algorithm 299
	X		Collected Algorithms for the CACM 1967 p. 243
	X	Updated for rounding errors based on remark in
	X		ACM TOMS June 1985, page 185
	X*/
	Xdouble
	Xpochisq (x, df)
	Xdouble	x;       /* obtained chi-square value */
	Xint 	df;      /* degrees of freedom */
	X	{
	X	double	a, y, s;
	X	double	e, c, z;
	X	double	poz ();   /* computes probability of normal z score */
	X	int 	even;     /* true if df is an even number */
	X	
	X	if (x <= 0.0 || df < 1)
	X		return (1.0);
	X	
	X	a = 0.5 * x;
	X	even = (2*(df/2)) == df;
	X	if (df > 1)
	X		y = ex (-a);
	X	s = (even ? y : (2.0 * poz (-sqrt (x))));
	X	if (df > 2)
	X		{
	X		x = 0.5 * (df - 1.0);
	X		z = (even ? 1.0 : 0.5);
	X		if (a > BIGX)
	X			{
	X			e = (even ? 0.0 : LOG_SQRT_PI);
	X			c = log (a);
	X			while (z <= x)
	X				{
	X				e = log (z) + e;
	X				s += ex (c*z-a-e);
	X				z += 1.0;
	X				}
	X			return (s);
	X			}
	X		else
	X			{
	X			e = (even ? 1.0 : (I_SQRT_PI / sqrt (a)));
	X			c = 0.0;
	X			while (z <= x)
	X				{
	X				e = e * (a / z);
	X				c = c + e;
	X				z += 1.0;
	X				}
	X			return (c * y + s);
	X			}
	X		}
	X	else
	X		return (s);
	X	}
	X
	X/*FUNCTION critchi: compute critical chi square value to produce given p */
	Xdouble
	Xcritchi (p, df)
	Xdouble	p;
	Xint 	df;
	X	{
	X	double	minchisq = 0.0;
	X	double	maxchisq = CHI_MAX;
	X	double	chisqval;
	X	
	X	if (p <= 0.0)
	X		return (maxchisq);
	X	else if (p >= 1.0)
	X		return (0.0);
	X	
	X	chisqval = df / sqrt (p);    /* fair first value */
	X	while (maxchisq - minchisq > CHI_EPSILON)
	X		{
	X		if (pochisq (chisqval, df) < p)
	X			maxchisq = chisqval;
	X		else
	X			minchisq = chisqval;
	X		chisqval = (maxchisq + minchisq) * 0.5;
	X		}
	X	return (chisqval);
	X	}
SHAR_EOF
if test 3003 -ne "`wc -c < 'chisq.c'`"
then
	echo shar: "error transmitting 'chisq.c'" '(should have been 3003 characters)'
fi
fi
echo shar: "extracting 'f.c'" '(3847 characters)'
if test -f 'f.c'
then
	echo shar: "will not over-write existing file 'f.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'f.c'
	X/*
	X	Module:       f.c
	X	Purpose:      compute approximations to F distribution probabilities
	X	Contents:     pof(), critf()
	X	Programmer:   Gary Perlman
	X	Organization: Wang Institute, Tyngsboro, MA 01879
	X	Tester:       compile with -DFTEST to include main program
	X	Copyright:    none
	X	Tabstops:     4
	X*/
	X
	X/*LINTLIBRARY*/
	X#include <math.h>
	Xstatic char sccsfid[] = "@(#) f.c 5.2 (|stat) 12/26/85";
	X
	X#ifndef	I_PI        /* 1 / pi */
	X#define	I_PI        0.3183098861837906715377675
	X#endif
	X#define	F_EPSILON     0.000001       /* accuracy of critf approximation */
	X#define	F_MAX      9999.0            /* maximum F ratio */
	X
	Xdouble	pof ();
	Xdouble	critf ();
	X
	X#ifdef	FTEST
	X
	Xint 	DFs[] = { 1, 2, 5, 10, 20, 40, 60, 120, -1 };
	Xdouble	Prob[] = { .10, .05, .01, .005, .001, -1.0 };
	X
	Xmain ()
	X	{
	X	int 	dfnum;
	X	int 	dfdenom;
	X	int 	p;
	X	
	X	for (p = 0; Prob[p] > 0.0; p++)
	X		{
	X		printf ("alpha = %g                      df1\n", Prob[p]);
	X		printf ("%-4s\\", "df2");
	X		for (dfnum = 0; DFs[dfnum] > 0; dfnum++)
	X			printf ("%7d ", DFs[dfnum]);
	X		putchar ('\n');
	X		for (dfdenom = 0; DFs[dfdenom] > 0; dfdenom++)
	X			{
	X			printf ("%4d ", DFs[dfdenom]);
	X			for (dfnum = 0; DFs[dfnum] > 0; dfnum++)
	X				printf ("%7.2f ", critf (Prob[p], DFs[dfnum], DFs[dfdenom]));
	X			putchar ('\n');
	X			}
	X		putchar ('\n');
	X		}
	X	}
	X#endif	FTEST
	X
	X/*FUNCTION pof: probability of F */
	X/*ALGORITHM Compute probability of F ratio.
	X	Adapted from Collected Algorithms of the CACM
	X	Algorithm 322
	X	Egon Dorrer
	X*/
	Xdouble
	Xpof (F, df1, df2)
	Xdouble	F;
	Xint 	df1, df2;
	X	{
	X	int	i, j;
	X	int	a, b;
	X	double	w, y, z, d, p;
	X	
	X	if (F < F_EPSILON || df1 < 1 || df2 < 1)
	X		return (1.0);
	X	a = df1%2 ? 1 : 2;
	X	b = df2%2 ? 1 : 2;
	X	w = (F * df1) / df2;
	X	z = 1.0 / (1.0 + w);
	X	if (a == 1)
	X		if (b == 1)
	X			{
	X			p = sqrt (w);
	X			y = I_PI; /* 1 / 3.14159 */
	X			d = y * z / p;
	X			p = 2.0 * y * atan (p);
	X			}
	X		else
	X			{
	X			p = sqrt (w * z);
	X			d = 0.5 * p * z / w;
	X			}
	X	else if (b == 1)
	X		{
	X		p = sqrt (z);
	X		d = 0.5 * z * p;
	X		p = 1.0 - p;
	X		}
	X	else
	X		{
	X		d = z * z;
	X		p = w * z;
	X		}
	X	y = 2.0 * w / z;
	X#ifdef	REMARK /* speedup modification suggested by Tolman (wrong answer!) */
	X	if (a == 1)
	X		for (j = b + 2; j <= df2; j += 2)
	X			{
	X			d *= (1.0 + a / (j - 2.0)) * z;
	X			p += d * y / (j - 1.0);
	X			}
	X	else
	X		{
	X		double	zk = 1.0;
	X		for (j = (df2 - 1) / 2; j; j--)
	X			zk *= z;
	X		d *= zk * df2/b;
	X		p *= zk + w * z * (zk - 1.0)/(z-1.0);
	X		}
	X#else /* original version */
	X	for (j = b + 2; j <= df2; j += 2)
	X		{
	X		d *= (1.0 + a / (j - 2.0)) * z;
	X		p = (a == 1 ? p + d * y / (j - 1.0) : (p + w) * z);
	X		}
	X#endif	REMARK
	X	y = w * z;
	X	z = 2.0 / z;
	X	b = df2 - 2;
	X	for (i = a + 2; i <= df1; i += 2)
	X		{
	X		j = i + b;
	X		d *= y * j / (i - 2.0);
	X		p -= z * d / j;
	X		}
	X	/* correction for approximation errors suggested in certification */
	X	if (p < 0.0)
	X		p = 0.0;
	X	else if (p > 1.0)
	X		p = 1.0;
	X	return (1.0-p);
	X	}
	X
	X/*FUNCTION critf: compute critical F value t produce given probability */
	X/*ALGORITHM
	X	Begin with upper and lower limits for F values (maxf and minf)
	X	set to extremes.  Choose an f value (fval) between the extremes.
	X	Compute the probability of the f value.  Set minf or maxf, based
	X	on whether the probability is less than or greater than the
	X	desired p.  Continue adjusting the extremes until they are
	X	within F_EPSILON of each other.
	X*/
	Xdouble
	Xcritf (p, df1, df2)
	Xdouble	p;
	Xint 	df1;
	Xint 	df2;
	X	{
	X	double	fval;
	X	double	fabs ();          /* floating absolute value */
	X	double	maxf = F_MAX;     /* maximum F ratio */
	X	double	minf = 0.0;       /* minimum F ratio */
	X	
	X	if (p <= 0.0 || p >= 1.0)
	X		return (0.0);
	X	
	X	fval = 1.0 / p;             /* the smaller the p, the larger the F */
	X	
	X	while (fabs (maxf - minf) > F_EPSILON)
	X		{
	X		if (pof (fval, df1, df2) < p)     /* F too large */
	X			maxf = fval;
	X		else                              /* F too small */
	X			minf = fval;
	X		fval = (maxf + minf) * 0.5;
	X		}
	X	
	X	return (fval);
	X	}
SHAR_EOF
if test 3847 -ne "`wc -c < 'f.c'`"
then
	echo shar: "error transmitting 'f.c'" '(should have been 3847 characters)'
fi
fi
exit 0
#	End of shell archive

moved to a/perlman