From 0a643bb9d0d222124ec23a3ba7505f22ce094f4f Mon Sep 17 00:00:00 2001
From: ceriel <none@none>
Date: Tue, 5 Dec 1995 12:29:36 +0000
Subject: [PATCH] Improved the pow() function to give more exact results

---
 lang/cem/libcc.ansi/math/pow.c | 96 +++++++++++++++++++++++++---------
 1 file changed, 71 insertions(+), 25 deletions(-)

diff --git a/lang/cem/libcc.ansi/math/pow.c b/lang/cem/libcc.ansi/math/pow.c
index 879ce85c4..639dc07e5 100644
--- a/lang/cem/libcc.ansi/math/pow.c
+++ b/lang/cem/libcc.ansi/math/pow.c
@@ -9,46 +9,92 @@
 #include	<math.h>
 #include	<float.h>
 #include	<errno.h>
-#include	"localmath.h"
+#include	<limits.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)) {
+	double	y_intpart, y_fractpart, fp;
+	int	negexp, negx;
+	int	ex, newexp;
+	unsigned long yi;
+
+	if (x == 1.0) return x;
+
+	if (x == 0 && y <= 0) {
 		errno = EDOM;
 		return 0;
 	}
 
-	if (x == 0) return x;
+	if (y == 0) return 1.0;
 
-	if (x < 0) {
-		if (modf(y/2.0, &dummy) != 0) {
-			/* y was odd */
-			result_neg = 1;
+	if (y < 0) {
+		y = -y;
+		negexp = 1;
+	}
+	else	negexp = 0;
+
+	y_fractpart = modf(y, &y_intpart);
+
+	if (y_fractpart != 0) {
+		if (x < 0) {
+			errno = EDOM;
+			return 0;
 		}
-		x = -x;
 	}
-	x = log(x);
 
+	negx = 0;
 	if (x < 0) {
 		x = -x;
-		y = -y;
+		negx = 1;
 	}
-	/* Beware of overflow in the multiplication */
-	if (x > 1.0 && y > DBL_MAX/x) {
-		errno = ERANGE;
-		return result_neg ? -HUGE_VAL : HUGE_VAL;
+
+	if (y_intpart > ULONG_MAX) {
+		if (negx && modf(y_intpart/2.0, &y_fractpart) == 0) {
+			negx = 0;
+		}
+
+		x = log(x);
+
+		/* Beware of overflow in the multiplication */
+		if (x > 1.0 && y > DBL_MAX/x) {
+			errno = ERANGE;
+			return HUGE_VAL;
+		}
+		if (negexp) y = -y;
+
+		if (negx) return -exp(x*y);
+		return exp(x * y);
 	}
 
-	x = exp(x * y);
-	return result_neg ? -x : x;
+	if (y_fractpart != 0) {
+		fp = exp(y_fractpart * log(x));
+	}
+	else	fp = 1.0;
+	yi = y_intpart;
+	if (! (yi & 1)) negx = 0;
+	x = frexp(x, &ex);
+	newexp = 0;
+	for (;;) {
+		if (yi & 1) {
+			fp *= x;
+			newexp += ex;
+		}
+		yi >>= 1;
+		if (yi == 0) break;
+		x *= x;
+		ex <<= 1;
+		if (x < 0.5) {
+			x += x;
+			ex -= 1;
+		}
+	}
+	if (negexp) {
+		fp = 1.0/fp;
+		newexp = -newexp;
+	}
+	if (negx) {
+		return -ldexp(fp, newexp);
+	}
+	return ldexp(fp, newexp);
 }
-- 
2.34.1