Improved the pow() function to give more exact results

diff --git a/lang/cem/libcc.ansi/math/pow.c b/lang/cem/libcc.ansi/math/pow.c
index 879ce85..639dc07 100644 (file)
--- 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);
 }