New, improved fmod.c

diff --git a/lang/cem/libcc.ansi/math/fmod.c b/lang/cem/libcc.ansi/math/fmod.c
index 1e0c178..7dd5091 100644 (file)
--- a/lang/cem/libcc.ansi/math/fmod.c
+++ b/lang/cem/libcc.ansi/math/fmod.c
@@ -1,33 +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
- */
+/* fmod function */
+/* Author Robert R. Hall (hall@crach.cts.com) */
+
 /* $Header$ */
 
-#include       <math.h>
-#include       <errno.h>
+#include <math.h>
+#include <errno.h>
 
 double
-fmod(double x, double y)
-{
-       double val;
-       double frac;
+(fmod)(double x, double y)
+{                              /* compute fmod(x, y) */
+       double t;
+       int n, neg;
+       int ychar, xchar;
 
-       if (y == 0) {
+       if (y == 0.0) {
                errno = EDOM;
-               return 0;
-       }
-       frac = modf( x / y, &val);
-
-       return frac * y;
+               return 0.0;
+       }
+                               /* fmod(finite, finite) */
+       if (y < 0.0) y = -y;
+       if (x < 0.0) x = -x, neg = 1;
+       else neg = 0;
+       t = frexp(y, &ychar);
+                               /* substract |y| until |x| < |y| */
 
-/*
-       val = x / y;
-       if (val > LONG_MIN && val < LONG_MAX) {
-               long i = val;
-               return x - i * y;
+       t = frexp(x, &xchar);
+       for (n = xchar - ychar; 0 <= n; --n) {
+                               /* try to substract |y|*2^n */
+           t = ldexp(y, n);
+           if (t <= x) x -= t;
        }
-*/
+       return (neg ? -x : x);
 }