From d43142d811334b5aec092d3d20f9ea37fcc4a974 Mon Sep 17 00:00:00 2001
From: eck <none@none>
Date: Mon, 18 Dec 1989 15:44:36 +0000
Subject: [PATCH] changed from Hart & Cheney to Cody & Waite

---
 lang/cem/libcc.ansi/math/.distr      |  22 +++++
 lang/cem/libcc.ansi/math/LIST        |   7 +-
 lang/cem/libcc.ansi/math/Makefile    |  31 ++++++
 lang/cem/libcc.ansi/math/asin.c      |  50 ++++++++--
 lang/cem/libcc.ansi/math/atan.c      | 117 ++++++++--------------
 lang/cem/libcc.ansi/math/atan2.c     |   2 +-
 lang/cem/libcc.ansi/math/exp.c       |  69 ++++++-------
 lang/cem/libcc.ansi/math/fmod.c      |  34 +++++++
 lang/cem/libcc.ansi/math/frexp.e     |  14 +--
 lang/cem/libcc.ansi/math/localmath.h |   2 +-
 lang/cem/libcc.ansi/math/log.c       |  60 ++++++------
 lang/cem/libcc.ansi/math/log10.c     |  10 +-
 lang/cem/libcc.ansi/math/modf.e      |  18 ++--
 lang/cem/libcc.ansi/math/pow.c       |  31 +++++-
 lang/cem/libcc.ansi/math/sin.c       | 119 ++++++++++-------------
 lang/cem/libcc.ansi/math/sinh.c      |  75 ++++++++++----
 lang/cem/libcc.ansi/math/sqrt.c      |   2 +-
 lang/cem/libcc.ansi/math/tan.c       | 140 ++++++++-------------------
 lang/cem/libcc.ansi/math/tanh.c      |  40 ++++++--
 19 files changed, 476 insertions(+), 367 deletions(-)
 create mode 100644 lang/cem/libcc.ansi/math/.distr
 create mode 100644 lang/cem/libcc.ansi/math/Makefile
 create mode 100644 lang/cem/libcc.ansi/math/fmod.c

diff --git a/lang/cem/libcc.ansi/math/.distr b/lang/cem/libcc.ansi/math/.distr
new file mode 100644
index 000000000..54981b057
--- /dev/null
+++ b/lang/cem/libcc.ansi/math/.distr
@@ -0,0 +1,22 @@
+LIST
+Makefile
+asin.c
+atan.c
+atan2.c
+ceil.c
+exp.c
+fabs.c
+floor.c
+fmod.c
+frexp.e
+ldexp.c
+localmath.h
+log.c
+log10.c
+modf.e
+pow.c
+sin.c
+sinh.c
+sqrt.c
+tan.c
+tanh.c
diff --git a/lang/cem/libcc.ansi/math/LIST b/lang/cem/libcc.ansi/math/LIST
index d0f012a18..835e11b14 100644
--- a/lang/cem/libcc.ansi/math/LIST
+++ b/lang/cem/libcc.ansi/math/LIST
@@ -1,10 +1,8 @@
-mlib
 localmath.h
 asin.c
 atan2.c
 atan.c
 ceil.c
-cosh.c
 fabs.c
 pow.c
 log10.c
@@ -12,8 +10,11 @@ log.c
 sin.c
 sinh.c
 sqrt.c
-ldexp.c
 tan.c
 tanh.c
 exp.c
+ldexp.c
+fmod.c
 floor.c
+frexp.e
+modf.e
diff --git a/lang/cem/libcc.ansi/math/Makefile b/lang/cem/libcc.ansi/math/Makefile
new file mode 100644
index 000000000..5da5edad8
--- /dev/null
+++ b/lang/cem/libcc.ansi/math/Makefile
@@ -0,0 +1,31 @@
+CFLAGS=-L -LIB
+
+.SUFFIXES: .o .e .c
+
+.e.o:
+	$(CC) $(CFLAGS) -c -o $@ $*.e
+
+clean:
+	rm -rf asin.o atan2.o atan.o ceil.o fabs.o pow.o log10.o \
+		log.o sin.o sinh.o sqrt.o tan.o tanh.o exp.o ldexp.o \
+		fmod.o floor.o frexp.o modf.o OLIST
+
+asin.o:
+atan2.o:
+atan.o:
+ceil.o:
+fabs.o:
+pow.o:
+log10.o:
+log.o:
+sin.o:
+sinh.o:
+sqrt.o:
+tan.o:
+tanh.o:
+exp.o:
+ldexp.o:
+fmod.o:
+floor.o:
+frexp.o:
+modf.o:
diff --git a/lang/cem/libcc.ansi/math/asin.c b/lang/cem/libcc.ansi/math/asin.c
index f5e5a15e5..c316b9b1f 100644
--- a/lang/cem/libcc.ansi/math/asin.c
+++ b/lang/cem/libcc.ansi/math/asin.c
@@ -6,31 +6,61 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<math.h>
+#include	<errno.h>
 #include	"localmath.h"
 
-
 static double
 asin_acos(double x, int cosfl)
 {
 	int negative = x < 0;
+	int     i;
+	double  g;
+	static double p[] = {
+		-0.27368494524164255994e+2,
+		 0.57208227877891731407e+2,
+		-0.39688862997540877339e+2,
+		 0.10152522233806463645e+2,
+		-0.69674573447350646411e+0
+	};
+	static double q[] = {
+		-0.16421096714498560795e+3,
+		 0.41714430248260412556e+3,
+		-0.38186303361750149284e+3,
+		 0.15095270841030604719e+3,
+		-0.23823859153670238830e+2,
+		 1.0
+	};
 
 	if (negative) {
 		x = -x;
 	}
-	if (x > 1) {
-		errno = EDOM;
-		return 0;
+	if (x > 0.5) {
+		i = 1;
+		if (x > 1) {
+			errno = EDOM;
+			return 0;
+		}
+		g = 0.5 - 0.5 * x;
+		x = - sqrt(g);
+		x += x;
 	}
-	if (x == 1) {
-		x = M_PI_2;
+	else {
+		/* ??? avoid underflow ??? */
+		i = 0;
+		g = x * x;
 	}
-	else x = atan(x/sqrt(1-x*x));
-	if (negative) x = -x;
+	x += x * g * POLYNOM4(g, p) / POLYNOM5(g, q);
 	if (cosfl) {
-		return M_PI_2 - x;
+		if (! negative) x = -x;
+	}
+	if ((cosfl == 0) == (i == 1)) {
+		x = (x + M_PI_4) + M_PI_4;
+	}
+	else if (cosfl && negative && i == 1) {
+		x = (x + M_PI_2) + M_PI_2;
 	}
+	if (! cosfl && negative) x = -x;
 	return x;
 }
 
diff --git a/lang/cem/libcc.ansi/math/atan.c b/lang/cem/libcc.ansi/math/atan.c
index 80de677f7..23cb63609 100644
--- a/lang/cem/libcc.ansi/math/atan.c
+++ b/lang/cem/libcc.ansi/math/atan.c
@@ -6,7 +6,6 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<float.h>
 #include	<math.h>
 #include	"localmath.h"
@@ -14,91 +13,55 @@
 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)
+	/*      Algorithm and coefficients from:
+			"Software manual for the elementary functions"
+			by W.J. Cody and W. Waite, Prentice-Hall, 1980
 	*/
-	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 p[] = {
+		-0.13688768894191926929e+2,
+		-0.20505855195861651981e+2,
+		-0.84946240351320683534e+1,
+		-0.83758299368150059274e+0
 	};
-
-	static double q[6] = {
-		0.7698297257888171026986294911e+03,
-		0.1813892701754635858982709369e+04,
-		0.1484049607102276827437401170e+04,
-		0.4904645326203706217748848797e+03,
-		0.5593479839280348664778328000e+02,
-		0.1000000000000000000000000000e+01
+	static double q[] = {
+		 0.41066306682575781263e+2,
+		 0.86157349597130242515e+2,
+		 0.59578436142597344465e+2,
+		 0.15024001160028576121e+2,
+		 1.0
+	};
+	static double a[] = {
+		0.0,
+		0.52359877559829887307710723554658381,  /* pi/6 */
+		M_PI_2,
+		1.04719755119659774615421446109316763   /* pi/3 */
 	};
 
-	int negative = x < 0.0;
-	register struct precomputed *pr = prec;
+	int     neg = x < 0;
+	int     n;
+	double  g;
 
-	if (negative) {
+	if (neg) {
 		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));
+	if (x > 1.0) {
+		x = 1.0/x;
+		n = 2;
 	}
-	else {
-		double xsq = x*x;
+	else    n = 0;
 
-		x = x * POLYNOM4(xsq, p)/POLYNOM5(xsq, q);
+	if (x > 0.26794919243112270647) {       /* 2-sqtr(3) */
+		n = n + 1;
+		x = (((0.73205080756887729353*x-0.5)-0.5)+x)/
+			(1.73205080756887729353+x);
 	}
-	return negative ? -x : x;
-}
 
+	/* ??? avoid underflow ??? */
+
+	g = x * x;
+	x += x * g * POLYNOM3(g, p) / POLYNOM4(g, q);
+	if (n > 1) x = -x;
+	x += a[n];
+	return neg ? -x : x;
+}
diff --git a/lang/cem/libcc.ansi/math/atan2.c b/lang/cem/libcc.ansi/math/atan2.c
index 541e91746..0e253c746 100644
--- a/lang/cem/libcc.ansi/math/atan2.c
+++ b/lang/cem/libcc.ansi/math/atan2.c
@@ -6,8 +6,8 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<math.h>
+#include	<errno.h>
 #include	"localmath.h"
 
 double
diff --git a/lang/cem/libcc.ansi/math/exp.c b/lang/cem/libcc.ansi/math/exp.c
index d731d6540..c5db4be0a 100644
--- a/lang/cem/libcc.ansi/math/exp.c
+++ b/lang/cem/libcc.ansi/math/exp.c
@@ -6,34 +6,35 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
-#include	<float.h>
 #include	<math.h>
+#include	<float.h>
+#include	<errno.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 */
+	/*	Algorithm and coefficients from:
+			"Software manual for the elementary functions"
+			by W.J. Cody and W. Waite, Prentice-Hall, 1980
+	*/
 
-	static double p[3] = {
-		 0.2080384346694663001443843411e+07,
-		 0.3028697169744036299076048876e+05,
-		 0.6061485330061080841615584556e+02
+	static double p[] = {
+		0.25000000000000000000e+0,
+		0.75753180159422776666e-2,
+		0.31555192765684646356e-4
 	};
 
-	static double q[4] = {
-		 0.6002720360238832528230907598e+07,
-		 0.3277251518082914423057964422e+06,
-		 0.1749287689093076403844945335e+04,
-		 0.1000000000000000000000000000e+01
+	static double q[] = {
+		0.50000000000000000000e+0,
+		0.56817302698551221787e-1,
+		0.63121894374398503557e-3,
+		0.75104028399870046114e-6
 	};
-
-	int negative = x < 0;
-	int ipart, large = 0;
-	double xsqr, xPxx, Qxx;
+	double	xn, g;
+	int	n;
+	int	negative = x < 0;
 
 	if (x <= M_LN_MIN_D) {
 		if (x < M_LN_MIN_D) errno = ERANGE;
@@ -44,22 +45,24 @@ exp(double x)
 		return DBL_MAX;
 	}
 
-	if (negative) {
-		x = -x;
+	if (negative) x = -x;
+ 
+	/* ??? avoid underflow ??? */
+
+	n = x * M_LOG2E + 0.5;	/* 1/ln(2) = log2(e), 0.5 added for rounding */
+	xn = n;
+	{
+		double	x1 = (long) x;
+		double	x2 = x - x1;
+
+		g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
 	}
-	x /= M_LN2;
-	ipart = floor(x);
-	x -= ipart;
-	if (x > 0.5) {
-		large = 1;
-		x -= 0.5;
+	if (negative) {
+		g = -g;
+		n = -n;
 	}
-	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;
+	xn = g * g;
+	x = g * POLYNOM2(xn, p);
+	n += 1;
+	return (ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
 }
diff --git a/lang/cem/libcc.ansi/math/fmod.c b/lang/cem/libcc.ansi/math/fmod.c
new file mode 100644
index 000000000..8117b2b6e
--- /dev/null
+++ b/lang/cem/libcc.ansi/math/fmod.c
@@ -0,0 +1,34 @@
+/*
+ * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ *
+ * Author: Hans van Eck
+ */
+/* $Header$ */
+
+#include	<math.h>
+#include	<errno.h>
+
+double
+fmod(double x, double y)
+{
+	long	i;
+	double val;
+	double frac;
+
+	if (y == 0) {
+		errno = EDOM;
+		return 0;
+	}
+	frac = modf( x / y, &val);
+
+	return frac * y;
+
+/*
+	val = x / y;
+	if (val > LONG_MIN && val < LONG_MAX) {
+		i = val;
+		return x - i * y;
+	}
+*/
+}
diff --git a/lang/cem/libcc.ansi/math/frexp.e b/lang/cem/libcc.ansi/math/frexp.e
index 922470c60..496dfa9d9 100644
--- a/lang/cem/libcc.ansi/math/frexp.e
+++ b/lang/cem/libcc.ansi/math/frexp.e
@@ -5,16 +5,16 @@
  */
 /* $Header$ */
 
- mes 2,EM_WSIZE,EM_PSIZE
+ 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
+ loi _EM_DSIZE
+ fef _EM_DSIZE
+ lal _EM_DSIZE
+ loi _EM_PSIZE
+ sti _EM_WSIZE
+ ret _EM_DSIZE
  end
 #endif
diff --git a/lang/cem/libcc.ansi/math/localmath.h b/lang/cem/libcc.ansi/math/localmath.h
index c73ab5742..0b2aaea83 100644
--- a/lang/cem/libcc.ansi/math/localmath.h
+++ b/lang/cem/libcc.ansi/math/localmath.h
@@ -39,4 +39,4 @@
 #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))
+#define	M_LN_MIN_D	(M_LN2 * (DBL_MIN_EXP - 1))
diff --git a/lang/cem/libcc.ansi/math/log.c b/lang/cem/libcc.ansi/math/log.c
index 53e709345..48f77a743 100644
--- a/lang/cem/libcc.ansi/math/log.c
+++ b/lang/cem/libcc.ansi/math/log.c
@@ -6,48 +6,54 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<math.h>
+#include	<errno.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)]
+	/*	Algorithm and coefficients from:
+			"Software manual for the elementary functions"
+			by W.J. Cody and W. Waite, Prentice-Hall, 1980
 	*/
-	/*	Hart & Cheney #2707 */
-
-	static double p[5] = {
-		 0.7504094990777122217455611007e+02,
-		-0.1345669115050430235318253537e+03,
-		 0.7413719213248602512779336470e+02,
-		-0.1277249755012330819984385000e+02,
-		 0.3327108381087686938144000000e+00
+	static double a[] = {
+		-0.64124943423745581147e2,
+		 0.16383943563021534222e2,
+		-0.78956112887491257267e0
 	};
-
-	static double q[5] = {
-		 0.3752047495388561108727775374e+02,
-		-0.7979028073715004879439951583e+02,
-		 0.5616126132118257292058560360e+02,
-		-0.1450868091858082685362325000e+02,
-		 0.1000000000000000000000000000e+01
+	static double b[] = {
+		-0.76949932108494879777e3,
+		 0.31203222091924532844e3,
+		-0.35667977739034646171e2,
+		 1.0
 	};
 
-	double z, zsqr;
-	int exponent;
+	double	znum, zden, z, w;
+	int	exponent;
 
-	if (x <= 0) {
+	if (x < 0) {
 		errno = EDOM;
-		return 0;
+		return -HUGE_VAL;
+	}
+	else if (x == 0) {
+		errno = ERANGE;
+		return -HUGE_VAL;
 	}
 
 	x = frexp(x, &exponent);
-	while (x < M_1_SQRT2) {
-		x += x;
+	if (x > M_1_SQRT2) {
+		znum = (x - 0.5) - 0.5;
+		zden = x * 0.5 + 0.5;
+	}
+	else {
+		znum = x - 0.5;
+		zden = znum * 0.5 + 0.5;
 		exponent--;
 	}
-	z = (x-1)/(x+1);
-	zsqr = z*z;
-	return z * POLYNOM4(zsqr, p) / POLYNOM4(zsqr, q) + exponent * M_LN2;
+	z = znum/zden; w = z * z;
+	x = z + z * w * (POLYNOM2(w,a)/POLYNOM3(w,b));
+	z = exponent;
+	x += z * (-2.121944400546905827679e-4);
+	return x + z * 0.693359375;
 }
diff --git a/lang/cem/libcc.ansi/math/log10.c b/lang/cem/libcc.ansi/math/log10.c
index 09ab6b882..147472201 100644
--- a/lang/cem/libcc.ansi/math/log10.c
+++ b/lang/cem/libcc.ansi/math/log10.c
@@ -6,16 +6,20 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<math.h>
+#include	<errno.h>
 #include	"localmath.h"
 
 double
 log10(double x)
 {
-	if (x <= 0) {
+	if (x < 0) {
 		errno = EDOM;
-		return 0;
+		return -HUGE_VAL;
+	}
+	else if (x == 0) {
+		errno = ERANGE;
+		return -HUGE_VAL;
 	}
 
 	return log(x) / M_LN10;
diff --git a/lang/cem/libcc.ansi/math/modf.e b/lang/cem/libcc.ansi/math/modf.e
index 8396117ba..51acd8e13 100644
--- a/lang/cem/libcc.ansi/math/modf.e
+++ b/lang/cem/libcc.ansi/math/modf.e
@@ -5,20 +5,20 @@
  */
 /* $Header$ */
 
- mes 2,EM_WSIZE,EM_PSIZE
+ mes 2,_EM_WSIZE,_EM_PSIZE
 #ifndef NOFLOAT
  exp $modf
  pro $modf,0
  lal 0
- loi EM_DSIZE
+ loi _EM_DSIZE
  loc 1
- loc EM_WSIZE
- loc EM_DSIZE
+ loc _EM_WSIZE
+ loc _EM_DSIZE
  cif
- fif EM_DSIZE
- lal EM_DSIZE
- loi EM_PSIZE
- sti EM_DSIZE
- ret EM_DSIZE
+ fif _EM_DSIZE
+ lal _EM_DSIZE
+ loi _EM_PSIZE
+ sti _EM_DSIZE
+ ret _EM_DSIZE
  end
 #endif
diff --git a/lang/cem/libcc.ansi/math/pow.c b/lang/cem/libcc.ansi/math/pow.c
index af04e3ea0..b9b7106b5 100644
--- a/lang/cem/libcc.ansi/math/pow.c
+++ b/lang/cem/libcc.ansi/math/pow.c
@@ -6,13 +6,21 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<math.h>
+#include	<float.h>
+#include	<errno.h>
+#include	"localmath.h"
 
 double
 pow(double x, double y)
 {
+	/*	Simple version for now. The Cody and Waite book has
+		a very complicated, much more precise version, but
+		this version has machine-dependent arrays A1 and A2,
+		and I don't know yet how to solve this ???
+	*/
 	double dummy;
+	int	result_neg = 0;
 
 	if ((x == 0 && y == 0) ||
 	    (x < 0 && modf(y, &dummy) != 0)) {
@@ -23,13 +31,26 @@ pow(double x, double y)
 	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;
+			result_neg = 1;
 		}
-		return val;
+		x = -x;
+	}
+	x = log(x);
+	if (x < 0) {
+		x = -x;
+		y = -y;
+	}
+	if (y > M_LN_MAX_D/x) {
+		errno = ERANGE;
+		return 0;
+	}
+	if (y < M_LN_MIN_D/x) {
+		errno = ERANGE;
+		return 0;
 	}
 
-	return exp(log(x) * y);
+	x = exp(x * y);
+	return result_neg ? -x : x;
 }
diff --git a/lang/cem/libcc.ansi/math/sin.c b/lang/cem/libcc.ansi/math/sin.c
index de60f99a1..cd3d204ff 100644
--- a/lang/cem/libcc.ansi/math/sin.c
+++ b/lang/cem/libcc.ansi/math/sin.c
@@ -6,93 +6,76 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<math.h>
 #include	"localmath.h"
 
 static double
-sinus(double x, int quadrant)
+sinus(double x, int cos_flag)
 {
-	/*	sin(0.5*pi*x) = x * P(x*x)/Q(x*x) for x in [0,1] */
-	/*	Hart & Cheney # 3374 */
+	/*      Algorithm and coefficients from:
+			"Software manual for the elementary functions"
+			by W.J. Cody and W. Waite, Prentice-Hall, 1980
+	*/
 
-	static double p[6] = {
-		 0.4857791909822798473837058825e+10,
-		-0.1808816670894030772075877725e+10,
-		 0.1724314784722489597789244188e+09,
-		-0.6351331748520454245913645971e+07,
-		 0.1002087631419532326179108883e+06,
-		-0.5830988897678192576148973679e+03
+	static double r[] = {
+		-0.16666666666666665052e+0,
+		 0.83333333333331650314e-2,
+		-0.19841269841201840457e-3,
+		 0.27557319210152756119e-5,
+		-0.25052106798274584544e-7,
+		 0.16058936490371589114e-9,
+		-0.76429178068910467734e-12,
+		 0.27204790957888846175e-14
 	};
 
-	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;
+	double  xsqr;
+	double  y;
+	int     neg = 0;
 
 	if (x < 0) {
-		quadrant += 2;
 		x = -x;
+		neg = 1;
 	}
-	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 (cos_flag) {
+		neg = 0;
+		y = M_PI_2 + x;
 	}
-	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;
+	else    y = x;
+
+	/* ??? avoid loss of significance, if y is too large, error ??? */
+
+	y = y * M_1_PI + 0.5;
+
+	/*      Use extended precision to calculate reduced argument.
+		Here we used 12 bits of the mantissa for a1.
+		Also split x in integer part x1 and fraction part x2.
+	*/
+#define A1 3.1416015625
+#define A2 -8.908910206761537356617e-6
+	{
+		double x1, x2;
+
+		modf(y, &y);
+		if (modf(0.5*y, &x1)) neg = !neg;
+		if (cos_flag) y -= 0.5;
+		x2 = modf(x, &x1);
+		x = x1 - y * A1;
 		x += x2;
-		x -= n * A2;
+		x -= y * 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) {
+ 
+	if (x < 0) {
+		neg = !neg;
 		x = -x;
 	}
-	xsqr = x * x;
-	x = x * POLYNOM5(xsqr, p) / POLYNOM5(xsqr, q);
-	return x;
+
+	/* ??? avoid underflow ??? */
+
+	y = x * x;
+	x += x * y * POLYNOM7(y, r);
+	return neg ? -x : x;
 }
 
 double
diff --git a/lang/cem/libcc.ansi/math/sinh.c b/lang/cem/libcc.ansi/math/sinh.c
index 39633dc52..ccb010e74 100644
--- a/lang/cem/libcc.ansi/math/sinh.c
+++ b/lang/cem/libcc.ansi/math/sinh.c
@@ -6,33 +6,72 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<math.h>
+#include	<float.h>
+#include	<errno.h>
 #include	"localmath.h"
 
-double
-sinh(double x)
+static double
+sinh_cosh(double x, int cosh_flag)
 {
-	int negx = x < 0;
+	/*	Algorithm and coefficients from:
+			"Software manual for the elementary functions"
+			by W.J. Cody and W. Waite, Prentice-Hall, 1980
+	*/
+
+	static double p[] = {
+		-0.35181283430177117881e+6,
+		-0.11563521196851768270e+5,
+		-0.16375798202630751372e+3,
+		-0.78966127417357099479e+0
+	};
+	static double q[] = {
+		-0.21108770058106271242e+7,
+		 0.36162723109421836460e+5,
+		-0.27773523119650701167e+3,
+		 1.0
+	};
+	int	negative = x < 0;
+	double	y = negative ? -x : x;
 
-	if (negx) {
-		x = -x;
+	if (! cosh_flag && y <= 1.0) {
+		/* ??? check for underflow ??? */
+		y = y * y;
+		return x + x * y * POLYNOM3(y, p)/POLYNOM3(y,q);
 	}
-	if (x > M_LN_MAX_D) {
-		/* exp(x) would overflow */
-		if (x >= M_LN_MAX_D + M_LN2) {
-			/* not representable */
-			x = HUGE_VAL;
+
+	if (y >= M_LN_MAX_D) {
+		/* exp(y) would cause overflow */
+#define LNV	0.69316101074218750000e+0
+#define VD2M1	0.52820835025874852469e-4
+		double	w = y - LNV;
+		
+		if (w < M_LN_MAX_D+M_LN2-LNV) {
+			x = exp(w);
+			x += VD2M1 * x;
+		}
+		else {
 			errno = ERANGE;
+			x = HUGE_VAL;
 		}
-		else	x = exp (x - M_LN2);
 	}
 	else {
-		double expx = exp(x);
-		x = 0.5 * (expx - 1.0/expx);
-	}
-	if (negx) {
-		return -x;
+		double	z = exp(y);
+		
+		x = 0.5 * (z + (cosh_flag ? 1.0 : -1.0)/z);
 	}
-	return x;
+	return negative ? -x : x;
+}
+
+double
+sinh(double x)
+{
+	return sinh_cosh(x, 0);
+}
+
+double
+cosh(double x)
+{
+	if (x < 0) x = -x;
+	return sinh_cosh(x, 1);
 }
diff --git a/lang/cem/libcc.ansi/math/sqrt.c b/lang/cem/libcc.ansi/math/sqrt.c
index afaf2a490..f41f3aafd 100644
--- a/lang/cem/libcc.ansi/math/sqrt.c
+++ b/lang/cem/libcc.ansi/math/sqrt.c
@@ -6,8 +6,8 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
 #include	<math.h>
+#include	<errno.h>
 
 #define NITER	5
 
diff --git a/lang/cem/libcc.ansi/math/tan.c b/lang/cem/libcc.ansi/math/tan.c
index 3d6911ff1..b125fd99c 100644
--- a/lang/cem/libcc.ansi/math/tan.c
+++ b/lang/cem/libcc.ansi/math/tan.c
@@ -6,117 +6,63 @@
  */
 /* $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)
+	/*      Algorithm and coefficients from:
+			"Software manual for the elementary functions"
+			by W.J. Cody and W. Waite, Prentice-Hall, 1980
 	*/
 
-	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;
+	double  y;
+	static double   p[] = {
+		 1.0,
+		-0.13338350006421960681e+0,
+		 0.34248878235890589960e-2,
+		-0.17861707342254426711e-4
+	};
+	static double   q[] = {
+		 1.0,
+		-0.46671683339755294240e+0,
+		 0.25663832289440112864e-1,
+		-0.31181531907010027307e-3,
+		 0.49819433993786512270e-6
+	};
 
 	if (negative) x = -x;
+ 
+	/* ??? avoid loss of significance, error if x is too large ??? */
 
-	/*	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;
+	y = x * M_2_PI + 0.5;
 
-	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;
+	/*      Use extended precision to calculate reduced argument.
+		Here we used 12 bits of the mantissa for a1.
+		Also split x in integer part x1 and fraction part x2.
+	*/
+    #define A1 1.57080078125
+    #define A2 -4.454455103380768678308e-6
+	{
+		double x1, x2;
+
+		modf(y, &y);
+		if (modf(0.5*y, &x1)) invert = 1;
+		x2 = modf(x, &x1);
+		x = x1 - y * A1;
+		x += x2;
+		x -= y * A2;
+    #undef A1
+    #undef A2
 	}
 
-	return negative ? -tmp : tmp;
+	/* ??? avoid underflow ??? */
+	y = x * x;
+	x += x * y * POLYNOM2(y, p+1);
+	y = POLYNOM4(y, q);
+	if (negative) x = -x;
+	return invert ? -y/x : x/y;
 }
diff --git a/lang/cem/libcc.ansi/math/tanh.c b/lang/cem/libcc.ansi/math/tanh.c
index cd1216a38..636a9df39 100644
--- a/lang/cem/libcc.ansi/math/tanh.c
+++ b/lang/cem/libcc.ansi/math/tanh.c
@@ -6,19 +6,45 @@
  */
 /* $Header$ */
 
-#include	<errno.h>
+#include	<float.h>
 #include	<math.h>
 #include	"localmath.h"
 
 double
 tanh(double x)
 {
-	if (x <= 0.5*M_LN_MIN_D) {
-		return -1;
-	}
+	/*	Algorithm and coefficients from:
+			"Software manual for the elementary functions"
+			by W.J. Cody and W. Waite, Prentice-Hall, 1980
+	*/
+
+	static double p[] = {
+		-0.16134119023996228053e+4,
+		-0.99225929672236083313e+2,
+		-0.96437492777225469787e+0
+	};
+	static double q[] = {
+		 0.48402357071988688686e+4,
+		 0.22337720718962312926e+4,
+		 0.11274474380534949335e+3,
+		 1.0
+	};
+	int 	negative = x < 0;
+
+	if (negative) x = -x;
+
 	if (x >= 0.5*M_LN_MAX_D) {
-		return 1;
+		x = 1.0;
+	}
+#define LN3D2	0.54930614433405484570e+0	/* ln(3)/2 */
+	else if (x > LN3D2) {
+		x = 0.5 - 1.0/(exp(x+x)+1.0);
+		x += x;
+	}
+	else {
+		/* ??? avoid underflow ??? */
+		double g = x*x;
+		x += x * g * POLYNOM2(g, p)/POLYNOM3(g, q);
 	}
-	x = exp(x + x);
-	return (x - 1.0)/(x + 1.0);
+	return negative ? -x : x;
 }
-- 
2.34.1