Initial revision

lang/cem/libcc.ansi/math/asin.c

@@ -0,0 +1,47 @@
+/*
+ * (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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+
+static double
+asin_acos(double x, int cosfl)
+{
+	int negative = x < 0;
+
+	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(double x)
+{
+	return asin_acos(x, 0);
+}
+
+double
+acos(double x)
+{
+	return asin_acos(x, 1);
+}

lang/cem/libcc.ansi/math/atan.c

@@ -0,0 +1,104 @@
+/*
+ * (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	<errno.h>
+#include	<float.h>
+#include	<math.h>
+#include	"localmath.h"
+
+double
+atan(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 },
+		{ DBL_MAX,
+		  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.ansi/math/atan2.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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+double
+atan2(double y, double x)
+{
+	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.ansi/math/ceil.c

@@ -0,0 +1,20 @@
+/*
+ * (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>
+
+double
+ceil(double x)
+{
+	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.ansi/math/cosh.c

@@ -0,0 +1,33 @@
+/*
+ * (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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+double
+cosh(double x)
+{
+	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_VAL;
+			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.ansi/math/exp.c

@@ -0,0 +1,65 @@
+/*
+ * (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	<errno.h>
+#include	<float.h>
+#include	<math.h>
+#include	"localmath.h"
+
+
+double
+exp(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;
+
+	if (x <= M_LN_MIN_D) {
+		if (x < M_LN_MIN_D) errno = ERANGE;
+		return DBL_MIN;
+	}
+	if (x >= M_LN_MAX_D) {
+		if (x > M_LN_MAX_D) errno = ERANGE;
+		return DBL_MAX;
+	}
+
+	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.ansi/math/fabs.c

@@ -0,0 +1,13 @@
+/*
+ * (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(double x)
+{
+	return  x < 0 ? -x : x;
+}

lang/cem/libcc.ansi/math/floor.c

@@ -0,0 +1,20 @@
+/*
+ * (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>
+
+double
+floor(double x)
+{
+	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.ansi/math/frexp.e

@@ -0,0 +1,20 @@
+#
+/*
+ * (c) copyright 1987 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ */
+/* $Header$ */
+
+ mes 2,EM_WSIZE,EM_PSIZE
+#ifndef NOFLOAT
+ exp $frexp
+ pro $frexp,0
+ lal 0
+ loi EM_DSIZE
+ fef EM_DSIZE
+ lal EM_DSIZE
+ loi EM_PSIZE
+ sti EM_WSIZE
+ ret EM_DSIZE
+ end
+#endif

lang/cem/libcc.ansi/math/ldexp.c

@@ -0,0 +1,36 @@
+/*
+ * (c) copyright 1987 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ */
+/* $Header$ */
+
+#include	<math.h>
+
+double
+ldexp(double fl, int exp)
+{
+	int sign = 1;
+	int currexp;
+
+	if (fl<0) {
+		fl = -fl;
+		sign = -1;
+	}
+	fl = frexp(fl,&currexp);
+	exp += currexp;
+	if (exp > 0) {
+		while (exp>30) {
+			fl *= (double) (1L << 30);
+			exp -= 30;
+		}
+		fl *= (double) (1L << exp);
+	}
+	else	{
+		while (exp<-30) {
+			fl /= (double) (1L << 30);
+			exp += 30;
+		}
+		fl /= (double) (1L << -exp);
+	}
+	return sign * fl;
+}

lang/cem/libcc.ansi/math/localmath.h

@@ -0,0 +1,42 @@
+/*
+ * localmath.h - This header is used by the mathematical library.
+ */
+/* $Header$ */
+
+/* some constants (Hart & Cheney) */
+#define	M_PI		3.14159265358979323846264338327950288
+#define	M_2PI		6.28318530717958647692528676655900576
+#define	M_3PI_4		2.35619449019234492884698253745962716
+#define	M_PI_2		1.57079632679489661923132169163975144
+#define	M_3PI_8		1.17809724509617246442349126872981358
+#define	M_PI_4		0.78539816339744830961566084581987572
+#define	M_PI_8		0.39269908169872415480783042290993786
+#define	M_1_PI		0.31830988618379067153776752674502872
+#define	M_2_PI		0.63661977236758134307553505349005744
+#define	M_4_PI		1.27323954473516268615107010698011488
+#define	M_E		2.71828182845904523536028747135266250
+#define	M_LOG2E		1.44269504088896340735992468100189213
+#define	M_LOG10E	0.43429448190325182765112891891660508
+#define	M_LN2		0.69314718055994530941723212145817657
+#define	M_LN10		2.30258509299404568401799145468436421
+#define	M_SQRT2		1.41421356237309504880168872420969808
+#define	M_1_SQRT2	0.70710678118654752440084436210484904
+#define	M_EULER		0.57721566490153286060651209008240243
+
+/* macros for constructing polynomials */
+#define	POLYNOM1(x, a)	((a)[1]*(x)+(a)[0])
+#define	POLYNOM2(x, a)	(POLYNOM1((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM3(x, a)	(POLYNOM2((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM4(x, a)	(POLYNOM3((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM5(x, a)	(POLYNOM4((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM6(x, a)	(POLYNOM5((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM7(x, a)	(POLYNOM6((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM8(x, a)	(POLYNOM7((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM9(x, a)	(POLYNOM8((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM10(x, a)	(POLYNOM9((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM11(x, a)	(POLYNOM10((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM12(x, a)	(POLYNOM11((x),(a)+1)*(x)+(a)[0])
+#define	POLYNOM13(x, a)	(POLYNOM12((x),(a)+1)*(x)+(a)[0])
+
+#define	M_LN_MAX_D	(M_LN2 * DBL_MAX_EXP)
+#define	M_LN_MIN_D	(M_LN2 * (DBL_MAX_EXP - 1))

lang/cem/libcc.ansi/math/log.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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+
+double
+log(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
+	};
+
+	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.ansi/math/log10.c

@@ -0,0 +1,22 @@
+/*
+ * (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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+double
+log10(double x)
+{
+	if (x <= 0) {
+		errno = EDOM;
+		return 0;
+	}
+
+	return log(x) / M_LN10;
+}

lang/cem/libcc.ansi/math/modf.e

@@ -0,0 +1,24 @@
+#
+/*
+ * (c) copyright 1987 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ */
+/* $Header$ */
+
+ mes 2,EM_WSIZE,EM_PSIZE
+#ifndef NOFLOAT
+ exp $modf
+ pro $modf,0
+ lal 0
+ loi EM_DSIZE
+ loc 1
+ loc EM_WSIZE
+ loc EM_DSIZE
+ cif
+ fif EM_DSIZE
+ lal EM_DSIZE
+ loi EM_PSIZE
+ sti EM_DSIZE
+ ret EM_DSIZE
+ end
+#endif

lang/cem/libcc.ansi/math/pow.c

@@ -0,0 +1,35 @@
+/*
+ * (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	<errno.h>
+#include	<math.h>
+
+double
+pow(double x, double y)
+{
+	double dummy;
+
+	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.ansi/math/sin.c

@@ -0,0 +1,109 @@
+/*
+ * (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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+static double
+sinus(double x, int quadrant)
+{
+	/*	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 {
+		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(double x)
+{
+	return sinus(x, 0);
+}
+
+double
+cos(double x)
+{
+	if (x < 0) x = -x;
+	return sinus(x, 1);
+}

lang/cem/libcc.ansi/math/sinh.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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+double
+sinh(double x)
+{
+	int negx = x < 0;
+
+	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_VAL;
+			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.ansi/math/sqrt.c

@@ -0,0 +1,36 @@
+/*
+ * (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	<errno.h>
+#include	<math.h>
+
+#define NITER	5
+
+double
+sqrt(double x)
+{
+	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.ansi/math/tan.c

@@ -0,0 +1,122 @@
+/*
+ * (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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+
+double
+tan(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 {
+		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_VAL;
+		}
+		else tmp = tmp2 / tmp1;
+	}
+
+	return negative ? -tmp : tmp;
+}

lang/cem/libcc.ansi/math/tanh.c

@@ -0,0 +1,24 @@
+/*
+ * (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	<errno.h>
+#include	<math.h>
+#include	"localmath.h"
+
+double
+tanh(double x)
+{
+	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);
+}