Initial revision

lang/cem/libcc/math/LIST

@@ -0,0 +1,22 @@
+tail_m.a
+asin.c
+atan2.c
+atan.c
+ceil.c
+cosh.c
+fabs.c
+gamma.c
+hypot.c
+jn.c
+j0.c
+j1.c
+log10.c
+pow.c
+log.c
+sin.c
+sinh.c
+sqrt.c
+tan.c
+tanh.c
+exp.c
+floor.c

lang/cem/libcc/math/asin.c

@@ -0,0 +1,53 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+
+static double
+asin_acos(x, cosfl)
+	double x;
+{
+	int negative = x < 0;
+	extern double sqrt(), atan();
+
+	if (negative) {
+		x = -x;
+	}
+	if (x > 1) {
+		errno = EDOM;
+		return 0;
+	}
+	if (x == 1) {
+		x = M_PI_2;
+	}
+	else x = atan(x/sqrt(1-x*x));
+	if (negative) x = -x;
+	if (cosfl) {
+		return M_PI_2 - x;
+	}
+	return x;
+}
+
+double
+asin(x)
+	double x;
+{
+	return asin_acos(x, 0);
+}
+
+double
+acos(x)
+	double x;
+{
+	return asin_acos(x, 1);
+}

lang/cem/libcc/math/atan.c

@@ -0,0 +1,103 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+double
+atan(x)
+	double x;
+{
+	/*	The interval [0, infinity) is treated as follows:
+		Define partition points Xi
+			X0 = 0
+			X1 = tan(pi/16)
+			X2 = tan(3pi/16)
+			X3 = tan(5pi/16)
+			X4 = tan(7pi/16)
+			X5 = infinity
+		and evaluation nodes xi
+			x2 = tan(2pi/16)
+			x3 = tan(4pi/16)
+			x4 = tan(6pi/16)
+			x5 = infinity
+		An argument x in [Xn-1, Xn] is now reduced to an argument
+		t in [-X1, X1] by the following formulas:
+			
+			t = 1/xn - (1/(xn*xn) + 1)/((1/xn) + x)
+
+			arctan(x) = arctan(xi) + arctan(t)
+
+		For the interval [0, p/16] an approximation is used:
+			arctan(x) = x * P(x*x)/Q(x*x)
+	*/
+	static struct precomputed {
+		double X;		/* partition point */
+		double arctan;		/* arctan of evaluation node */
+		double one_o_x;		/* 1 / xn */
+		double one_o_xsq_p_1;	/* 1 / (xn*xn) + 1 */
+	} prec[5] = {
+		{ 0.19891236737965800691159762264467622,
+		  0.0,
+		  0.0,		/* these don't matter */
+		  0.0 } ,
+		{ 0.66817863791929891999775768652308076, /* tan(3pi/16)	*/
+		  M_PI_8,
+		  2.41421356237309504880168872420969808,
+		  6.82842712474619009760337744841939616 },
+		{ 1.49660576266548901760113513494247691, /* tan(5pi/16) */
+		  M_PI_4,
+		  1.0,
+		  2.0 },
+		{ 5.02733949212584810451497507106407238, /* tan(7pi/16) */
+		  M_3PI_8,
+		  0.41421356237309504880168872420969808,
+		  1.17157287525380998659662255158060384 },
+		{ MAXDOUBLE,
+		  M_PI_2,
+		  0.0,
+		  1.0 }};
+
+	/*	Hart & Cheney # 5037 */
+
+	static double p[5] = {
+		0.7698297257888171026986294745e+03,
+		0.1557282793158363491416585283e+04,
+		0.1033384651675161628243434662e+04,
+		0.2485841954911840502660889866e+03,
+		0.1566564964979791769948970100e+02
+	};
+
+	static double q[6] = {
+		0.7698297257888171026986294911e+03,
+		0.1813892701754635858982709369e+04,
+		0.1484049607102276827437401170e+04,
+		0.4904645326203706217748848797e+03,
+		0.5593479839280348664778328000e+02,
+		0.1000000000000000000000000000e+01
+	};
+
+	int negative = x < 0.0;
+	register struct precomputed *pr = prec;
+
+	if (negative) {
+		x = -x;
+	}
+	while (x > pr->X) pr++;
+	if (pr != prec) {
+		x = pr->arctan +
+			atan(pr->one_o_x - pr->one_o_xsq_p_1/(pr->one_o_x + x));
+	}
+	else {
+		double xsq = x*x;
+
+		x = x * POLYNOM4(xsq, p)/POLYNOM5(xsq, q);
+	}
+	return negative ? -x : x;
+}

lang/cem/libcc/math/atan2.c

@@ -0,0 +1,46 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+atan2(y, x)
+	double x, y;
+{
+	extern double atan();
+	double absx, absy, val;
+
+	if (x == 0 && y == 0) {
+		errno = EDOM;
+		return 0;
+	}
+	absy = y < 0 ? -y : y;
+	absx = x < 0 ? -x : x;
+	if (absy - absx == absy) {
+		/* x negligible compared to y */
+		return y < 0 ? -M_PI_2 : M_PI_2;
+	}
+	if (absx - absy == absx) {
+		/* y negligible compared to x */
+		val = 0.0;
+	}
+	else	val = atan(y/x);
+	if (x > 0) {
+		/* first or fourth quadrant; already correct */
+		return val;
+	}
+	if (y < 0) {
+		/* third quadrant */
+		return val - M_PI;
+	}
+	return val + M_PI;
+}

lang/cem/libcc/math/ceil.c

@@ -0,0 +1,21 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+double
+ceil(x)
+	double x;
+{
+	extern double modf();
+	double val;
+
+	return modf(x, &val) > 0 ? val + 1.0 : val ;
+	/*	this also works if modf always returns a positive
+		fractional part
+	*/
+}

lang/cem/libcc/math/cosh.c

@@ -0,0 +1,38 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+cosh(x)
+	double x;
+{
+	extern double exp();
+
+	if (x < 0) {
+		x = -x;
+	}
+	if (x > M_LN_MAX_D) {
+		/* exp(x) would overflow */
+		if (x >= M_LN_MAX_D + M_LN2) {
+			/* not representable */
+			x = HUGE;
+			errno = ERANGE;
+		}
+		else	x = exp (x - M_LN2);
+	}
+	else {
+		double expx = exp(x);
+		x = 0.5 * (expx + 1.0/expx);
+	}
+	return x;
+}

lang/cem/libcc/math/exp.c

@@ -0,0 +1,67 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+exp(x)
+	double x;
+{
+	/*	2**x = (Q(x*x)+x*P(x*x))/(Q(x*x)-x*P(x*x)) for x in [0,0.5] */
+	/*	Hart & Cheney #1069 */
+
+	static double p[3] = {
+		 0.2080384346694663001443843411e+07,
+		 0.3028697169744036299076048876e+05,
+		 0.6061485330061080841615584556e+02
+	};
+
+	static double q[4] = {
+		 0.6002720360238832528230907598e+07,
+		 0.3277251518082914423057964422e+06,
+		 0.1749287689093076403844945335e+04,
+		 0.1000000000000000000000000000e+01
+	};
+
+	int negative = x < 0;
+	int ipart, large = 0;
+	double xsqr, xPxx, Qxx;
+	extern double floor(), ldexp();
+
+	if (x <= M_LN_MIN_D) {
+		if (x < M_LN_MIN_D) errno = ERANGE;
+		return M_MIN_D;
+	}
+	if (x >= M_LN_MAX_D) {
+		if (x < M_LN_MAX_D) errno = ERANGE;
+		return M_MAX_D;
+	}
+
+	if (negative) {
+		x = -x;
+	}
+	x /= M_LN2;
+	ipart = floor(x);
+	x -= ipart;
+	if (x > 0.5) {
+		large = 1;
+		x -= 0.5;
+	}
+	xsqr = x * x;
+	xPxx = x * POLYNOM2(xsqr, p);
+	Qxx = POLYNOM3(xsqr, q);
+	x = (Qxx + xPxx) / (Qxx - xPxx);
+	if (large) x *= M_SQRT2;
+	x = ldexp(x, ipart);
+	if (negative) return 1.0/x;
+	return x;
+}

lang/cem/libcc/math/fabs.c

@@ -0,0 +1,15 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+double
+fabs(x)
+	double x;
+{
+	return  x < 0 ? -x : x;
+}

lang/cem/libcc/math/floor.c

@@ -0,0 +1,21 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+double
+floor(x)
+	double x;
+{
+	extern double modf();
+	double val;
+
+	return modf(x, &val) < 0 ? val - 1.0 : val ;
+	/*	this also works if modf always returns a positive
+		fractional part
+	*/
+}

lang/cem/libcc/math/gamma.c

@@ -0,0 +1,137 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+static double
+smallpos_gamma(x)
+	double x;
+{
+	/*	Approximation of gamma function using
+		gamma(x) = P(x-2) / Q(x-2) for x in [2,3]
+	*/
+	/* Hart & Cheney # 5251 */
+
+	static double p[11] = {
+		-0.2983543278574342138830437659e+06,
+		-0.2384953970018198872468734423e+06,
+		-0.1170494760121780688403854445e+06,
+		-0.3949445048301571936421824091e+05,
+		-0.1046699423827521405330650531e+05,
+		-0.2188218110071816359394795998e+04,
+		-0.3805112208641734657584922631e+03,
+		-0.5283123755635845383718978382e+02,
+		-0.6128571763704498306889428212e+01,
+		-0.5028018054416812467364198750e+00,
+		-0.3343060322330595274515660112e-01
+	};
+
+	static double q[9] = {
+		-0.2983543278574342138830438524e+06,
+		-0.1123558608748644911342306408e+06,
+		 0.5332716689118142157485686311e+05,
+		 0.8571160498907043851961147763e+04,
+		-0.4734865977028211706556819770e+04,
+		 0.1960497612885585838997039621e+03,
+		 0.1257733367869888645966647426e+03,
+		-0.2053126153100672764513929067e+02,
+		 0.1000000000000000000000000000e+01
+	};
+
+	double result = 1.0;
+
+	while (x > 3) {
+		x -= 1.0;
+		result *= x;
+	}
+	while (x < 2) {
+		result /= x;
+		x += 1.0;
+	}
+
+	x -= 2.0;
+
+	return result * POLYNOM10(x, p) / POLYNOM8(x, q);
+}
+
+#define log_sqrt_2pi 0.91893853320467274178032973640561763
+
+int	signgam;
+
+static double
+bigpos_loggamma(x)
+	double x;
+{
+	/*	computes the log(gamma(x)) function for big arguments
+		using the Stirling form
+		  log(gamma(x)) = (x - 0.5)log(x) - x + log(sqrt(2*pi)) + fi(x)
+		where fi(x) = (1/x)*P(1/(x*x))/Q(1/(x*x)) for x in [12,1000]
+	*/
+	/* Hart & Cheney # 5468 */
+
+	static double p[4] = {
+		 0.12398282342474941538685913e+00,
+		 0.67082783834332134961461700e+00,
+		 0.64507302912892202513890000e+00,
+		 0.66662907040200752600000000e-01
+	};
+
+	static double q[4] = {
+		 0.14877938810969929846815600e+01,
+		 0.80995271894897557472821400e+01,
+		 0.79966911236636441947720000e+01,
+		 0.10000000000000000000000000e+01
+	};
+
+	double rsq = 1.0/(x*x);
+	extern double log();
+
+	return (x-0.5)*log(x)-x+log_sqrt_2pi+POLYNOM3(rsq, p)/(x*POLYNOM3(rsq, q));
+}
+
+static double
+neg_loggamma(x)
+	double x;
+{
+	/*	compute the log(gamma(x)) function for negative values of x,
+		using the rule:
+			-x*gamma(x)*gamma(-x) = pi/sin(z*pi)
+	*/
+	extern double sin(), log();
+	double sinpix;
+
+	x = -x;
+	sinpix = sin(M_PI * x);
+	if (sinpix == 0.0) {
+		errno = EDOM;
+		return HUGE;
+	}
+	if (sinpix < 0) sinpix = -sinpix;
+	else signgam = -1;
+	return log(M_PI/(x * smallpos_gamma(x) * sinpix));
+}
+
+double
+gamma(x)
+	double x;
+{
+	/*	Wrong name; Actually computes log(gamma(x))
+	*/
+	extern double log();
+
+	signgam = 1;
+	if (x <= 0) {
+		return neg_loggamma(x);
+	}
+	if (x > 12.0) {
+		return bigpos_loggamma(x);
+	}
+	return log(smallpos_gamma(x));
+}

lang/cem/libcc/math/hypot.c

@@ -0,0 +1,39 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+double
+hypot(x,y)
+	double x,y;
+{
+	/*	Computes sqrt(x*x+y*y), avoiding overflow */
+
+	extern double sqrt();
+
+	if (x < 0) x = -x;
+	if (y < 0) y = -y;
+	if (x > y) {
+		double t = y;
+		y = x;
+		x = t;
+	}
+	/* sqrt(x*x+y*y) = sqrt(y*y*(x*x/(y*y)+1.0)) = y*sqrt(x*x/(y*y)+1.0) */
+	x /= y;
+	return y*sqrt(x*x+1.0);
+}
+
+struct complex {
+	double r,i;
+};
+
+double
+cabs(p_compl)
+	struct complex p_compl;
+{
+	return hypot(p_compl.r, p_compl.i);
+}

lang/cem/libcc/math/j0.c

@@ -0,0 +1,203 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+
+static double
+P0(x)
+	double x;
+{
+	/*	P0(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
+	/*	Hart & Cheney # 6554 */
+
+	static double p[9] = {
+		 0.9999999999999999999999995647e+00,
+		 0.5638253933310769952531889297e+01,
+		 0.1124846237418285392887270013e+02,
+		 0.1009280644639441488899111404e+02,
+		 0.4290591487686900980651458361e+01,
+		 0.8374209971661497198619102718e+00,
+		 0.6702347074465611456598882534e-01,
+		 0.1696260729396856143084502774e-02,
+		 0.6463970103128382090713889584e-05
+	};
+
+	static double q[9] = {
+		 0.9999999999999999999999999999e+00,
+		 0.5639352566123269952531467562e+01,
+		 0.1125463057106955935416066535e+02,
+		 0.1010501892629524191262518048e+02,
+		 0.4301396985171094350444425443e+01,
+		 0.8418926086780046799127094223e+00,
+		 0.6784915305473610998681570734e-01,
+		 0.1754416614608056207958880988e-02,
+		 0.7482977995134121064747276923e-05
+	};
+
+	double zsq = 64.0/(x*x);
+
+	return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
+}
+
+static double
+Q0(x)
+	double x;
+{
+	/*	Q0(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
+	/*	Hart & Cheney # 6955 */
+	/*	Probably typerror in Hart & Cheney; it sais:
+		Q0(x) = x*P(z*z)/Q(z*z)
+	*/
+
+	static double p[9] = {
+		-0.1562499999999999999999995808e-01,
+		-0.1111285583113679178917024959e+00,
+		-0.2877685516355036842789761274e+00,
+		-0.3477683453166454475665803194e+00,
+		-0.2093031978191084473537206358e+00,
+		-0.6209520943730206312601003832e-01,
+		-0.8434508346572023650653353729e-02,
+		-0.4414848186188819989871882393e-03,
+		-0.5768946278415631134804064871e-05
+	};
+
+	static double q[10] = {
+		 0.9999999999999999999999999999e+00,
+		 0.7121383005365046745065850254e+01,
+		 0.1848194194302368046679068851e+02,
+		 0.2242327522435983712994071530e+02,
+		 0.1359286169255959339963319677e+02,
+		 0.4089489268101204780080944780e+01,
+		 0.5722140925672174525430730669e+00,
+		 0.3219814230905924725810683346e-01,
+		 0.5299687475496044642364124073e-03,
+		 0.9423249021001925212258428217e-06
+	};
+
+	double zsq = 64.0/(x*x);
+
+	return (8.0/x) * POLYNOM8(zsq, p) / POLYNOM9(zsq, q);
+}
+
+static double
+smallj0(x)
+	double x;
+{
+	/*	J0(x) = P(x*x)/Q(x*x) for x in [0,8] */
+	/*	Hart & Cheney # 5852 */
+
+	static double p[10] = {
+		 0.1641556014884554385346147435e+25,
+		-0.3943559664767296636012616471e+24,
+		 0.2172018385924539313982287997e+23,
+		-0.4814859952069817648285245941e+21,
+		 0.5345457598841972345381674607e+19,
+		-0.3301538925689637686465426220e+17,
+		 0.1187390681211042949874031474e+15,
+		-0.2479851167896144439689877514e+12,
+		 0.2803148940831953934479400118e+09,
+		-0.1336625500481224741885945416e+06
+	};
+
+	static double q[10] = {
+		 0.1641556014884554385346137617e+25,
+		 0.1603303724440893273539045602e+23,
+		 0.7913043777646405204323616203e+20,
+		 0.2613165313325153278086066185e+18,
+		 0.6429607918826017759289213100e+15,
+		 0.1237672982083407903483177730e+13,
+		 0.1893012093677918995179541438e+10,
+		 0.2263381356781110003609399116e+07,
+		 0.1974019272727281783930443513e+04,
+		 0.1000000000000000000000000000e+01
+	};
+
+	double xsq = x*x;
+
+	return POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
+}
+
+double
+j0(x)
+	double x;
+{
+	/*	Use J0(x) = sqrt(2/(pi*x))*(P0(x)*cos(X0)-Q0(x)*sin(X0))
+			where X0 = x - pi/4 for |x| > 8.
+		Use J0(-x) = J0(x).
+		Use direct approximation of smallj0 for |x| <= 8.
+	*/
+	extern double sqrt(), sin(), cos();
+
+	if (x < 0) x = -x;
+	if (x > 8.0) {
+		double X0 = x - M_PI_4;
+		return sqrt(M_2_PI/x)*(P0(x)*cos(X0) - Q0(x)*sin(X0));
+	}
+	return smallj0(x);
+}
+
+static double
+smally0_bar(x)
+	double x;
+{
+	/*	Y0(x) = Y0BAR(x)+(2/pi)*J0(x)ln(x)
+		Approximation of Y0BAR for 0 <= x <= 8:
+			Y0BAR(x) = P(x*x)/Q(x*x)
+		Hart & Cheney #6250
+	*/
+
+	static double p[14] = {
+		-0.2692670958801060448840356941e+14,
+		 0.6467231173109037044444917683e+14,
+		-0.5563036156275660297303897296e+13,
+		 0.1698403391975239335187832821e+12,
+		-0.2606282788256139370857687880e+10,
+		 0.2352841334491277505699488812e+08,
+		-0.1365184412186963659690851354e+06,
+		 0.5371538422626582142170627457e+03,
+		-0.1478903875146718839145348490e+01,
+		 0.2887840299886172125955719069e-02,
+		-0.3977426824263991024666116123e-05,
+		 0.3738169731655229006655176866e-08,
+		-0.2194460874896856106887900645e-11,
+		 0.6208996973821484304384239393e-15
+	};
+
+	static double q[6] = {
+		 0.3648393301278364629844168660e+15,
+		 0.1698390180526960997295118328e+13,
+		 0.3587111679107612117789088586e+10,
+		 0.4337760840406994515845890005e+07,
+		 0.3037977771964348276793136205e+04,
+		 0.1000000000000000000000000000e+01
+	};
+
+	double xsq = x*x;
+
+	return POLYNOM13(xsq, p) / POLYNOM5(xsq, q);
+}
+
+double
+y0(x)
+	double x;
+{
+	extern double sqrt(), sin(), cos(), log();
+
+	if (x <= 0.0) {
+		errno = EDOM;
+		return -HUGE;
+	}
+	if (x > 8.0) {
+		double X0 = x - M_PI_4;
+		return sqrt(M_2_PI/x) * (P0(x)*sin(X0)+Q0(x)*cos(X0));
+	}
+	return smally0_bar(x) + M_2_PI*j0(x)*log(x);
+}

lang/cem/libcc/math/j1.c

@@ -0,0 +1,206 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+
+static double
+P1(x)
+	double x;
+{
+	/*	P1(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
+	/*	Hart & Cheney # 6755 */
+
+	static double p[9] = {
+		 0.1000000000000000000000000489e+01,
+		 0.5581663300347182292169450071e+01,
+		 0.1100186625131173123750501118e+02,
+		 0.9727139359130463694593683431e+01,
+		 0.4060011483142278994462590992e+01,
+		 0.7742832212665311906917358099e+00,
+		 0.6021617752811098752098248630e-01,
+		 0.1482350677236405118074646993e-02,
+		 0.6094215148131061431667573909e-05
+	};
+
+	static double q[9] = {
+		 0.9999999999999999999999999999e+00,
+		 0.5579832245659682292169922224e+01,
+		 0.1099168447731617288972771040e+02,
+		 0.9707206835125961446797916892e+01,
+		 0.4042610016540342097334497865e+01,
+		 0.7671965204303836019508430169e+00,
+		 0.5893258668794493100786371406e-01,
+		 0.1393993644981256852404222530e-02,
+		 0.4585597769784750669754696825e-05
+	};
+
+	double zsq = 64.0/(x*x);
+
+	return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
+}
+
+static double
+Q1(x)
+	double x;
+{
+	/*	Q1(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
+	/*	Hart & Cheney # 7157 */
+	/*	Probably typerror in Hart & Cheney; it sais:
+		Q1(x) = x*P(z*z)/Q(z*z)
+	*/
+
+	static double p[9] = {
+		0.4687499999999999999999995275e-01,
+		0.3302394516691663879252493748e+00,
+		0.8456888491208195767613862428e+00,
+		0.1008551084218946085420665147e+01,
+		0.5973407972399900690521296181e+00,
+		0.1737697433393258207540273097e+00,
+		0.2303862814819568573893610740e-01,
+		0.1171224207976250587945594946e-02,
+		0.1486418220337492918307904804e-04
+	};
+
+	static double q[10] = {
+		0.9999999999999999999999999999e+00,
+		0.7049380763213049609070823421e+01,
+		0.1807129960468949760845562209e+02,
+		0.2159171174362827330505421695e+02,
+		0.1283239297740546866114600499e+02,
+		0.3758349275324260869598403931e+01,
+		0.5055985453754739528620657666e+00,
+		0.2665604326323907148063400439e-01,
+		0.3821140353404633025596424652e-03,
+		0.3206696590241261037875154062e-06
+	};
+
+	double zsq = 64.0/(x*x);
+
+	return (8.0/x) * POLYNOM8(zsq, p) / POLYNOM9(zsq, q);
+}
+
+static double
+smallj1(x)
+	double x;
+{
+	/*	J1(x) = x*P(x*x)/Q(x*x) for x in [0,8] */
+	/*	Hart & Cheney # 6054 */
+
+	static double p[10] = {
+		 0.1921176307760798128049021316e+25,
+		-0.2226092031387396254771375773e+24,
+		 0.7894463902082476734673226741e+22,
+		-0.1269424373753606065436561036e+21,
+		 0.1092152214043184787101134641e+19,
+		-0.5454629264396819144157448868e+16,
+		 0.1634659487571284628830445048e+14,
+		-0.2909662785381647825756152444e+11,
+		 0.2853433451054763915026471449e+08,
+		-0.1197705712815379389149134705e+05
+	};
+
+	static double q[10] = {
+		 0.3842352615521596256098041912e+25,
+		 0.3507567066272028105798868716e+23,
+		 0.1611334311633414344007062889e+21,
+		 0.4929612313959850319632645381e+18,
+		 0.1117536965288162684489793105e+16,
+		 0.1969278625584719037168592923e+13,
+		 0.2735606122949877990248154504e+10,
+		 0.2940957355049651347475558106e+07,
+		 0.2274736606126590905134610965e+04,
+		 0.1000000000000000000000000000e+01
+	};
+
+	double xsq = x*x;
+
+	return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
+}
+
+double
+j1(x)
+	double x;
+{
+	/*	Use J1(x) = sqrt(2/(pi*x))*(P1(x)*cos(X1)-Q1(x)*sin(X1))
+			where X1 = x - 3*pi/4 for |x| > 8.
+		Use J1(-x) = -J1(x).
+		Use direct approximation of smallj1 for |x| <= 8.
+	*/
+	extern double sqrt(), sin(), cos();
+	int negative = x < 0.0;
+
+	if (negative) x = -x;
+	if (x > 8.0) {
+		double X1 = x - (M_PI - M_PI_4);
+		x = sqrt(M_2_PI/x)*(P1(x)*cos(X1) - Q1(x)*sin(X1));
+	}
+	else x = smallj1(x);
+	if (negative) return -x;
+	return x;
+}
+
+static double
+smally1_bar(x)
+	double x;
+{
+	/*	Y1(x) = Y1BAR(x)+(2/pi)*(J1(x)ln(x) - 1/x)
+		Approximation of Y1BAR for 0 <= x <= 8:
+			Y1BAR(x) = x*P(x*x)/Q(x*x)
+		Hart & Cheney # 6449
+	*/
+
+	static double p[10] = {
+		-0.5862655424363443992938931700e+24,
+		 0.1570668341992328458208364904e+24,
+		-0.7351681299005467428400402479e+22,
+		 0.1390658785759080111485190942e+21,
+		-0.1339544201526785345938109179e+19,
+		 0.7290257386242270629526344379e+16,
+		-0.2340575603057015935501295099e+14,
+		 0.4411516199185230690878878903e+11,
+		-0.4542128738770213026987060358e+08,
+		 0.1988612563465350530472715888e+05
+	};
+
+	static double q[10] = {
+		 0.2990279721605116022908679994e+25,
+		 0.2780285010357803058127175655e+23,
+		 0.1302687474507355553192845146e+21,
+		 0.4071330372239164349602952937e+18,
+		 0.9446611865086570116528399283e+15,
+		 0.1707657951197456205887347694e+13,
+		 0.2440358986882941823431612517e+10,
+		 0.2708852767034077697963790196e+07,
+		 0.2174361138333330803617969305e+04,
+		 0.1000000000000000000000000000e+01
+	};
+
+	double xsq = x*x;
+
+	return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
+}
+
+double
+y1(x)
+	double x;
+{
+	extern double sqrt(), sin(), cos(), log();
+
+	if (x <= 0.0) {
+		errno = EDOM;
+		return -HUGE;
+	}
+	if (x > 8.0) {
+		double X1 = x - (M_PI - M_PI_4);
+		return sqrt(M_2_PI/x) * (P1(x)*sin(X1)+Q1(x)*cos(X1));
+	}
+	return smally1_bar(x) + M_2_PI*(j1(x)*log(x) - 1/x);
+}

lang/cem/libcc/math/jn.c

@@ -0,0 +1,121 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+double
+yn(n, x)
+	double x;
+{
+	/*	Use y0, y1, and the recurrence relation
+		y(n+1,x) = 2*n*y(n,x)/x - y(n-1, x).
+		According to Hart & Cheney, this is stable for all
+		x, n.
+		Also use: y(-n,x) = (-1)^n * y(n, x)
+	*/
+
+	int negative = 0;
+	extern double y0(), y1();
+	double yn1, yn2;
+	register int i;
+
+	if (x <= 0) {
+		errno = EDOM;
+		return -HUGE;
+	}
+
+	if (n < 0) {
+		n = -n;
+		negative = (n % 2);
+	}
+
+	if (n == 0) return y0(x);
+	if (n == 1) return y1(x);
+
+	yn2 = y0(x);
+	yn1 = y1(x);
+	for (i = 1; i < n; i++) {
+		double tmp = yn1;
+		yn1 = (i*2)*yn1/x - yn2;
+		yn2 = tmp;
+	}
+	if (negative) return -yn1;
+	return yn1;
+}
+
+double
+jn(n, x)
+	double x;
+{
+	/*	Unfortunately, according to Hart & Cheney, the recurrence
+		j(n+1,x) = 2*n*j(n,x)/x - j(n-1,x) is unstable for
+		increasing n, except when x > n.
+		However, j(n,x)/j(n-1,x) = 2/(2*n-x*x/(2*(n+1)-x*x/( .... 
+		(a continued fraction).
+		We can use this to determine KJn and KJn-1, where K is a
+		normalization constant not yet known. This enables us
+		to determine KJn-2, ...., KJ1, KJ0. Now we can use the
+		J0 or J1 approximation to determine K.
+		Use: j(-n, x) = (-1)^n * j(n, x)
+		     j(n, -x) = (-1)^n * j(n, x)
+	*/
+
+	extern double j0(), j1();
+
+	if (n < 0) {
+		n = -n;
+		x = -x;
+	}
+
+	if (n == 0) return j0(x);
+	if (n == 1) return j1(x);
+	if (x > n) {
+		/* in this case, the recurrence relation is stable for
+		   increasing n, so we use that.
+		*/
+		double jn2 = j0(x), jn1 = j1(x);
+		register int i;
+
+		for (i = 1; i < n; i++) {
+			double tmp = jn1;
+			jn1 = (2*i)*jn1/x - jn2;
+			jn2 = tmp;
+		}
+		return jn1;
+	}
+	{
+		/* we first compute j(n,x)/j(n-1,x) */
+		register int i;
+		double quotient = 0.0;
+		double xsqr = x*x;
+		double jn1, jn2;
+
+		for (i = 20;	/* ??? how many do we need ??? */
+		     i > 0; i--) {
+			quotient = xsqr/(2*(i+n) - quotient);
+		}
+		quotient = x / (2*n - quotient);
+
+		jn1 = quotient;
+		jn2 = 1.0;
+		for (i = n-1; i > 0; i--) {
+			/* recurrence relation is stable for decreasing n
+			*/
+			double tmp = jn2;
+			jn2 = (2*i)*jn2/x - jn1;
+			jn1 = tmp;
+		}
+		/* So, now we have K*Jn = quotient and K*J0 = jn2.
+		   Now it is easy; compute real j0, this gives K = jn2/j0,
+		   and this then gives Jn = quotient/K = j0 * quotient / jn2.
+		*/
+		return j0(x)*quotient/jn2;
+	}
+}

lang/cem/libcc/math/log.c

@@ -0,0 +1,56 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+log(x)
+	double x;
+{
+	/* log(x) = z*P(z*z)/Q(z*z), z = (x-1)/(x+1), x in [1/sqrt(2), sqrt(2)]
+	*/
+	/*	Hart & Cheney #2707 */
+
+	static double p[5] = {
+		 0.7504094990777122217455611007e+02,
+		-0.1345669115050430235318253537e+03,
+		 0.7413719213248602512779336470e+02,
+		-0.1277249755012330819984385000e+02,
+		 0.3327108381087686938144000000e+00
+	};
+
+	static double q[5] = {
+		 0.3752047495388561108727775374e+02,
+		-0.7979028073715004879439951583e+02,
+		 0.5616126132118257292058560360e+02,
+		-0.1450868091858082685362325000e+02,
+		 0.1000000000000000000000000000e+01
+	};
+
+	extern double frexp();
+	double z, zsqr;
+	int exponent;
+
+	if (x <= 0) {
+		errno = EDOM;
+		return 0;
+	}
+
+	x = frexp(x, &exponent);
+	while (x < M_1_SQRT2) {
+		x += x;
+		exponent--;
+	}
+	z = (x-1)/(x+1);
+	zsqr = z*z;
+	return z * POLYNOM4(zsqr, p) / POLYNOM4(zsqr, q) + exponent * M_LN2;
+}

lang/cem/libcc/math/log10.c

@@ -0,0 +1,27 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+log10(x)
+	double x;
+{
+	extern double log();
+
+	if (x <= 0) {
+		errno = EDOM;
+		return 0;
+	}
+
+	return log(x) / M_LN10;
+}

lang/cem/libcc/math/pow.c

@@ -0,0 +1,40 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+pow(x,y)
+	double x,y;
+{
+	double dummy;
+	extern double modf(), exp(), log();
+
+	if ((x == 0 && y == 0) ||
+	    (x < 0 && modf(y, &dummy) != 0)) {
+		errno = EDOM;
+		return 0;
+	}
+
+	if (x == 0) return x;
+
+	if (x < 0) {
+		double val = exp(log(-x) * y);
+		if (modf(y/2.0, &dummy) != 0) {
+			/* y was odd */
+			val = - val;
+		}
+		return val;
+	}
+
+	return exp(log(x) * y);
+}

lang/cem/libcc/math/sin.c

@@ -0,0 +1,115 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+static double
+sinus(x, quadrant)
+	double x;
+{
+	/*	sin(0.5*pi*x) = x * P(x*x)/Q(x*x) for x in [0,1] */
+	/*	Hart & Cheney # 3374 */
+
+	static double p[6] = {
+		 0.4857791909822798473837058825e+10,
+		-0.1808816670894030772075877725e+10,
+		 0.1724314784722489597789244188e+09,
+		-0.6351331748520454245913645971e+07,
+		 0.1002087631419532326179108883e+06,
+		-0.5830988897678192576148973679e+03
+	};
+
+	static double q[6] = {
+		 0.3092566379840468199410228418e+10,
+		 0.1202384907680254190870913060e+09,
+		 0.2321427631602460953669856368e+07,
+		 0.2848331644063908832127222835e+05,
+		 0.2287602116741682420054505174e+03,
+		 0.1000000000000000000000000000e+01
+	};
+
+	double xsqr;
+	int t;
+
+	if (x < 0) {
+		quadrant += 2;
+		x = -x;
+	}
+	if (M_PI_2 - x == M_PI_2) {
+		switch(quadrant) {
+		case 0:
+		case 2:
+			return 0.0;
+		case 1:
+			return 1.0;
+		case 3:
+			return -1.0;
+		}
+	}
+	if (x >= M_2PI) {
+	    if (x <= 0x7fffffff) {
+		/*	Use extended precision to calculate reduced argument.
+			Split 2pi in 2 parts a1 and a2, of which the first only
+			uses some bits of the mantissa, so that n * a1 is
+			exactly representable, where n is the integer part of
+			x/pi.
+			Here we used 12 bits of the mantissa for a1.
+			Also split x in integer part x1 and fraction part x2.
+			We then compute x-n*2pi as ((x1 - n*a1) + x2) - n*a2.
+		*/
+#define A1 6.2822265625
+#define A2 0.00095874467958647692528676655900576
+		double n = (long) (x / M_2PI);
+		double x1 = (long) x;
+		double x2 = x - x1;
+		x = x1 - n * A1;
+		x += x2;
+		x -= n * A2;
+#undef A1
+#undef A2
+	    }
+	    else {
+		extern double modf();
+		double dummy;
+
+		x = modf(x/M_2PI, &dummy) * M_2PI;
+	    }
+	}
+	x /= M_PI_2;
+	t = x;
+	x -= t;
+	quadrant = (quadrant + (int)(t % 4)) % 4;
+	if (quadrant & 01) {
+		x = 1 - x;
+	}
+	if (quadrant > 1) {
+		x = -x;
+	}
+	xsqr = x * x;
+	x = x * POLYNOM5(xsqr, p) / POLYNOM5(xsqr, q);
+	return x;
+}
+
+double
+sin(x)
+	double x;
+{
+	return sinus(x, 0);
+}
+
+double
+cos(x)
+	double x;
+{
+	if (x < 0) x = -x;
+	return sinus(x, 1);
+}

lang/cem/libcc/math/sinh.c

@@ -0,0 +1,42 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+sinh(x)
+	double x;
+{
+	int negx = x < 0;
+	extern double exp();
+
+	if (negx) {
+		x = -x;
+	}
+	if (x > M_LN_MAX_D) {
+		/* exp(x) would overflow */
+		if (x >= M_LN_MAX_D + M_LN2) {
+			/* not representable */
+			x = HUGE;
+			errno = ERANGE;
+		}
+		else	x = exp (x - M_LN2);
+	}
+	else {
+		double expx = exp(x);
+		x = 0.5 * (expx - 1.0/expx);
+	}
+	if (negx) {
+		return -x;
+	}
+	return x;
+}

lang/cem/libcc/math/sqrt.c

@@ -0,0 +1,41 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+#define NITER	5
+
+double
+sqrt(x)
+	double x;
+{
+	extern double frexp(), ldexp();
+	int exponent;
+	double val;
+
+	if (x <= 0) {
+		if (x < 0) errno = EDOM;
+		return 0;
+	}
+
+	val = frexp(x, &exponent);
+	if (exponent & 1) {
+		exponent--;
+		val *= 2;
+	}
+	val = ldexp(val + 1.0, exponent/2 - 1);
+	/* was: val = (val + 1.0)/2.0; val = ldexp(val, exponent/2); */
+	for (exponent = NITER - 1; exponent >= 0; exponent--) {
+		val = (val + x / val) / 2.0;
+	}
+	return val;
+}

lang/cem/libcc/math/tan.c

@@ -0,0 +1,126 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+tan(x)
+	double x;
+{
+	/*	First reduce range to [0, pi/4].
+		Then use approximation tan(x*pi/4) = x * P(x*x)/Q(x*x).
+		Hart & Cheney # 4288
+		Use: tan(x) = 1/tan(pi/2 - x) 
+		     tan(-x) = -tan(x)
+		     tan(x+k*pi) = tan(x)
+	*/
+
+	static double p[5] = {
+		-0.5712939549476836914932149599e+10,
+		 0.4946855977542506692946040594e+09,
+		-0.9429037070546336747758930844e+07,
+		 0.5282725819868891894772108334e+05,
+		-0.6983913274721550913090621370e+02
+	};
+
+	static double q[6] = {
+		-0.7273940551075393257142652672e+10,
+		 0.2125497341858248436051062591e+10,
+		-0.8000791217568674135274814656e+08,
+		 0.8232855955751828560307269007e+06,
+		-0.2396576810261093558391373322e+04,
+		 0.1000000000000000000000000000e+01
+	};
+
+	int negative = x < 0;
+	double tmp, tmp1, tmp2;
+	double xsq;
+	int invert = 0;
+	int ip;
+
+	if (negative) x = -x;
+
+	/*	first reduce to [0, pi)	*/
+	if (x >= M_PI) {
+	    if (x <= 0x7fffffff) {
+		/*	Use extended precision to calculate reduced argument.
+			Split pi in 2 parts a1 and a2, of which the first only
+			uses some bits of the mantissa, so that n * a1 is
+			exactly representable, where n is the integer part of
+			x/pi.
+			Here we used 12 bits of the mantissa for a1.
+			Also split x in integer part x1 and fraction part x2.
+			We then compute x-n*pi as ((x1 - n*a1) + x2) - n*a2.
+		*/
+#define A1 3.14111328125
+#define A2 0.00047937233979323846264338327950288
+		double n = (long) (x / M_PI);
+		double x1 = (long) x;
+		double x2 = x - x1;
+		x = x1 - n * A1;
+		x += x2;
+		x -= n * A2;
+#undef A1
+#undef A2
+	    }
+	    else {
+		extern double modf();
+
+		x = modf(x/M_PI, &tmp) * M_PI;
+	    }
+	}
+	/*	because the approximation uses x*pi/4, we reverse this */
+	x /= M_PI_4;
+	ip = (int) x;
+	x -= ip;
+
+	switch(ip) {
+	case 0:
+		/* [0,pi/4] */
+		break;
+	case 1:
+		/* [pi/4, pi/2]
+		   tan(x+pi/4) = 1/tan(pi/2 - (x+pi/4)) = 1/tan(pi/4 - x)
+		*/
+		invert = 1;
+		x = 1.0 - x;
+		break;
+	case 2:
+		/* [pi/2, 3pi/4]
+		   tan(x+pi/2) = tan((x+pi/2)-pi) = -tan(pi/2 - x) =
+		   -1/tan(x)
+		*/
+		negative = ! negative;
+		invert = 1;
+		break;
+	case 3:
+		/* [3pi/4, pi)
+		   tan(x+3pi/4) = tan(x-pi/4) = - tan(pi/4-x)
+		*/
+		x = 1.0 - x;
+		negative = ! negative;
+		break;
+	}
+	xsq = x * x;
+	tmp1 = x*POLYNOM4(xsq, p);
+	tmp2 = POLYNOM5(xsq, q);
+	tmp = tmp1 / tmp2;
+	if (invert) {
+		if (tmp == 0.0) {
+			errno = ERANGE;
+			tmp = HUGE;
+		}
+		else tmp = tmp2 / tmp1;
+	}
+
+	return negative ? -tmp : tmp;
+}

lang/cem/libcc/math/tanh.c

@@ -0,0 +1,27 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+/* $Header$ */
+
+#include <math.h>
+#include <errno.h>
+
+double
+tanh(x)
+	double x;
+{
+	extern double exp();
+
+	if (x <= 0.5*M_LN_MIN_D) {
+		return -1;
+	}
+	if (x >= 0.5*M_LN_MAX_D) {
+		return 1;
+	}
+	x = exp(x + x);
+	return (x - 1.0)/(x + 1.0);
+}

lang/cem/libcc/math/test.c

@@ -0,0 +1,193 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Ceriel J.H. Jacobs
+ */
+
+#include <math.h>
+#include <stdio.h>
+
+#define EPS_D	5.0e-14
+main()
+{
+	testsqrt();
+	testtrig();
+	testexplog();
+	testgamma();
+	testbessel();
+}
+
+dotest(s, x, d, v)
+	char *s;
+	double x, d, v;
+{
+	double fabs();
+
+	if (fabs((v - d) / (fabs(v) < EPS_D ? 1.0 : v)) > EPS_D) {
+		printf(s, x);
+		printf(" = %.16e, should be %.16e\n", d, v);
+	}
+}
+
+testsqrt()
+{
+#define SQRT2	M_SQRT2
+#define SQRT10	3.16227766016837933199889354443271853
+
+	double x, val;
+	extern double sqrt();
+
+	dotest("sqrt(%.1f)", 2.0, sqrt(2.0), SQRT2);
+	dotest("sqrt(%.1f)", 10.0, sqrt(10.0), SQRT10);
+
+	for (x = 0.1; x < 0.1e20; x += x) {
+		val = sqrt(x);
+		dotest("sqrt(%.1f)^2", x, val*val, x);
+	}
+}
+
+testtrig()
+{
+#define SINPI_24	0.13052619222005159154840622789548901
+#define SINPI_16	0.19509032201612826784828486847702224
+#define SINPI_12	0.25881904510252076234889883762404832
+#define SINPI_6		0.5
+#define SINPI_4		M_1_SQRT2
+#define SINPI_3		0.86602540378443864676372317075293618
+#define SINPI_2		1.0
+#define SIN0		0.0
+
+	double x;
+	extern double sin(), cos(), tan(), asin(), acos(), atan(), fabs();
+
+	dotest("sin(0)", 0.0, sin(0.0), SIN0);
+	dotest("sin(pi/24)", M_PI/24 , sin(M_PI/24), SINPI_24);
+	dotest("sin(pi/16)", M_PI/16 , sin(M_PI/16), SINPI_16);
+	dotest("sin(pi/12)", M_PI/12 , sin(M_PI/12), SINPI_12);
+	dotest("sin(pi/6)", M_PI/6 , sin(M_PI/6), SINPI_6);
+	dotest("sin(pi/4)", M_PI_4 , sin(M_PI_4), SINPI_4);
+	dotest("sin(pi/3)", M_PI/3 , sin(M_PI/3), SINPI_3);
+	dotest("sin(pi/2)", M_PI_2 , sin(M_PI_2), SINPI_2);
+	dotest("sin(pi)", 0.0, sin(M_PI), SIN0);
+	dotest("sin(3*pi/2)", 0.0, sin(M_PI+M_PI_2), -SINPI_2);
+
+	dotest("sin(-pi/24)", -M_PI/24 , sin(-M_PI/24), -SINPI_24);
+	dotest("sin(-pi/16)", -M_PI/16 , sin(-M_PI/16), -SINPI_16);
+	dotest("sin(-pi/12)", -M_PI/12 , sin(-M_PI/12), -SINPI_12);
+	dotest("sin(-pi/6)", -M_PI/6 , sin(-M_PI/6), -SINPI_6);
+	dotest("sin(-pi/4)", -M_PI_4 , sin(-M_PI_4), -SINPI_4);
+	dotest("sin(-pi/3)", -M_PI/3 , sin(-M_PI/3), -SINPI_3);
+	dotest("sin(-pi/2)", -M_PI_2 , sin(-M_PI_2), -SINPI_2);
+
+	dotest("cos(pi/2)", M_PI_2, cos(M_PI_2), SIN0);
+	dotest("cos(11pi/24)", M_PI/24 , cos(11*M_PI/24), SINPI_24);
+	dotest("cos(7pi/16)", M_PI/16 , cos(7*M_PI/16), SINPI_16);
+	dotest("cos(5pi/12)", M_PI/12 , cos(5*M_PI/12), SINPI_12);
+	dotest("cos(pi/3)", M_PI/6 , cos(M_PI/3), SINPI_6);
+	dotest("cos(pi/4)", M_PI_4 , cos(M_PI_4), SINPI_4);
+	dotest("cos(pi/6)", M_PI/3 , cos(M_PI/6), SINPI_3);
+	dotest("cos(0)", M_PI_2 , cos(0), SINPI_2);
+	dotest("cos(pi)", M_PI , cos(M_PI), -SINPI_2);
+	dotest("cos(3pi/2)", M_PI , cos(M_PI+M_PI_2), SIN0);
+
+	dotest("cos(-pi/2)", M_PI_2, cos(-M_PI_2), SIN0);
+	dotest("cos(-11pi/24)", M_PI/24 , cos(-11*M_PI/24), SINPI_24);
+	dotest("cos(-7pi/16)", M_PI/16 , cos(-7*M_PI/16), SINPI_16);
+	dotest("cos(-5pi/12)", M_PI/12 , cos(-5*M_PI/12), SINPI_12);
+	dotest("cos(-pi/3)", M_PI/6 , cos(-M_PI/3), SINPI_6);
+	dotest("cos(-pi/4)", M_PI_4 , cos(-M_PI_4), SINPI_4);
+	dotest("cos(-pi/6)", M_PI/3 , cos(-M_PI/6), SINPI_3);
+
+	for (x = -10; x <= 10; x += 0.5) {
+		dotest("sin+2*pi-sin(%.2f)", x, sin(x+M_2PI)-sin(x), 0.0);
+		dotest("cos+2*pi-cos(%.2f)", x, cos(x+M_2PI)-cos(x), 0.0);
+		dotest("tan+2*pi-tan(%.2f)", x, tan(x+M_2PI)-tan(x), 0.0);
+		dotest("tan+pi-tan(%.2f)", x, tan(x+M_PI)-tan(x), 0.0);
+	}
+
+	for (x = -1.5; x <= 1.5; x += 0.1) {
+		dotest("asin(sin(%.2f))", x, asin(sin(x)), x);
+		dotest("acos(cos(%.2f))", x, acos(cos(x)), fabs(x));
+		dotest("atan(tan(%.2f))", x, atan(tan(x)), x);
+	}
+}
+
+testexplog()
+{
+#define EXPMIN1		0.36787944117144232159552377016146087	/* exp(-1) */
+#define EXPMIN1_4	0.77880078307140486824517026697832065	/* exp(-1/4) */
+#define EXP0		1.0					/* exp(0) */
+#define EXP1_4		1.28402541668774148407342056806243646	/* exp(1/4) */
+#define EXP1		M_E					/* exp(1) */
+#define LN1		0.0					/* log(1) */
+#define LN2		M_LN2					/* log(2) */
+#define LN4		1.38629436111989061883446424291635313	/* log(4) */
+#define LNE		1.0					/* log(e) */
+#define LN10		M_LN10					/* log(10) */
+
+	extern double exp(), log();
+	double x;
+
+	dotest("exp(%.2f)", -1.0, exp(-1.0), EXPMIN1);
+	dotest("exp(%.2f)", -0.25, exp(-0.25), EXPMIN1_4);
+	dotest("exp(%.2f)", 0.0, exp(0.0), EXP0);
+	dotest("exp(%.2f)", 0.25, exp(0.25), EXP1_4);
+	dotest("exp(%.2f)", 1.0, exp(1.0), EXP1);
+
+	dotest("log(%.2f)", 1.0, log(1.0), LN1);
+	dotest("log(%.2f)", 2.0, log(2.0), LN2);
+	dotest("log(%.2f)", 4.0, log(4.0), LN4);
+	dotest("log(%.2f)", 10.0, log(10.0), LN10);
+	dotest("log(e)", M_E, log(M_E), LNE);
+
+	for (x = -30.0; x <= 30.0; x += 0.5) {
+		dotest("log(exp(%.2f))", x, log(exp(x)), x);
+	}
+}
+
+testgamma()
+{
+	double x, xfac;
+	extern double gamma(), exp();
+
+	for (x = 1.0, xfac = 1.0; x < 30.0; x += 1.0) {
+		dotest("exp(gamma(%.2f))", x, exp(gamma(x)), xfac);
+		xfac *= x;
+	}
+}
+
+testbessel()
+{
+#define J0__PI_4	0.85163191370480801270040601506092607 /* j0(pi/4) */
+#define J0__PI_2	0.47200121576823476744766838787250096 /* j0(pi/2) */
+#define J1__PI_4	0.36318783834686733179559374778892472 /* j1(pi/4) */
+#define J1__PI_2	0.56682408890587393771124496346716028 /* j1(pi/2) */
+#define J10__PI_4	0.00000000002369974904082422018721148 /* j10(p1/4) */
+#define J10__PI_2	0.00000002326614794865976450546482206 /* j10(pi/2) */
+
+	extern double j0(), j1(), jn(), yn();
+	register int n;
+	double x;
+	extern char *sprintf();
+	char buf[100];
+
+	dotest("j0(pi/4)", M_PI_4, j0(M_PI_4), J0__PI_4);
+	dotest("j0(pi/2)", M_PI_2, j0(M_PI_2), J0__PI_2);
+	dotest("j1(pi/4)", M_PI_4, j1(M_PI_4), J1__PI_4);
+	dotest("j1(pi/2)", M_PI_2, j1(M_PI_2), J1__PI_2);
+	dotest("j10(pi/4)", M_PI_4, jn(10,M_PI_4), J10__PI_4);
+	dotest("j10(pi/2)", M_PI_2, jn(10,M_PI_2), J10__PI_2);
+
+	/* Also check consistency using the Wronskian relation
+		jn(n+1,x)*yn(n, x) - jn(n,x)*yn(n+1,x) = 2/(pi*x)
+	*/
+
+	for (x = 0.1; x < 20.0; x += 0.5) {
+		double two_over_pix = M_2_PI/x;
+
+		for (n = 0; n <= 10; n++) {
+			dotest(sprintf(buf, "jn(%d,%.2f)*yn(%d,%.2f)-jn(%d,%.2f)*yn(%d,%.2f)",n+1,x,n,x,n,x,n+1,x), x, jn(n+1,x)*yn(n,x)-jn(n,x)*yn(n+1,x),M_2_PI/x);
+		}
+	}
+}