X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Fremainder.c;h=951040ceffca31bf801f109650aabe2df4c4cb93;hb=4779b635ef35c7b0bc4044fcb5bc746d06f158c4;hp=16f09e5db0ce8d058a657316207570c2be6eeee8;hpb=0e1c6ff93f27c939ba9e0df945b16ef98eaaeef1;p=gnulib.git

diff --git a/lib/remainder.c b/lib/remainder.c
index 16f09e5db..951040cef 100644
--- a/lib/remainder.c
+++ b/lib/remainder.c
@@ -1,5 +1,5 @@
 /* Remainder.
-   Copyright (C) 2012 Free Software Foundation, Inc.
+   Copyright (C) 2012-2013 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -14,33 +14,93 @@
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
 
-#include <config.h>
+#if ! (defined USE_LONG_DOUBLE || defined USE_FLOAT)
+# include <config.h>
+#endif
 
 /* Specification.  */
 #include <math.h>
 
-double
-remainder (double x, double y)
+#ifdef USE_LONG_DOUBLE
+# define REMAINDER remainderl
+# define DOUBLE long double
+# define L_(literal) literal##L
+# define FABS fabsl
+# define FMOD fmodl
+# define ISNAN isnanl
+#elif ! defined USE_FLOAT
+# define REMAINDER remainder
+# define DOUBLE double
+# define L_(literal) literal
+# define FABS fabs
+# define FMOD fmod
+# define ISNAN isnand
+#else /* defined USE_FLOAT */
+# define REMAINDER remainderf
+# define DOUBLE float
+# define L_(literal) literal##f
+# define FABS fabsf
+# define FMOD fmodf
+# define ISNAN isnanf
+#endif
+
+#undef NAN
+#if defined _MSC_VER
+static DOUBLE zero;
+# define NAN (zero / zero)
+#else
+# define NAN (L_(0.0) / L_(0.0))
+#endif
+
+DOUBLE
+REMAINDER (DOUBLE x, DOUBLE y)
 {
-  double q = - round (x / y);
-  double r = fma (q, y, x); /* = x + q * y, computed in one step */
-  /* Correct possible rounding errors in the quotient x / y.  */
-  double half_y = 0.5L * y;
-  if (y >= 0)
+  if (isfinite (x) && isfinite (y) && y != L_(0.0))
     {
-      /* Expect -y/2 <= r <= y/2.  */
-      if (r > half_y)
-        q -= 1, r = fma (q, y, x);
-      else if (r < - half_y)
-        q += 1, r = fma (q, y, x);
+      if (x == L_(0.0))
+        /* Return x, regardless of the sign of y.  */
+        return x;
+
+      {
+        int negate = ((!signbit (x)) ^ (!signbit (y)));
+        DOUBLE r;
+
+        /* Take the absolute value of x and y.  */
+        x = FABS (x);
+        y = FABS (y);
+
+        /* Trivial case that requires no computation.  */
+        if (x <= L_(0.5) * y)
+          return (negate ? - x : x);
+
+        /* With a fixed y, the function x -> remainder(x,y) has a period 2*y.
+           Therefore we can reduce the argument x modulo 2*y.  And it's no
+           problem if 2*y overflows, since fmod(x,Inf) = x.  */
+        x = FMOD (x, L_(2.0) * y);
+
+        /* Consider the 3 cases:
+             0 <= x <= 0.5 * y
+             0.5 * y < x < 1.5 * y
+             1.5 * y <= x <= 2.0 * y  */
+        if (x <= L_(0.5) * y)
+          r = x;
+        else
+          {
+            r = x - y;
+            if (r > L_(0.5) * y)
+              r = x - L_(2.0) * y;
+          }
+        return (negate ? - r : r);
+      }
     }
   else
     {
-      /* Expect y/2 <= r <= -y/2.  */
-      if (r > - half_y)
-        q += 1, r = fma (q, y, x);
-      else if (r < half_y)
-        q -= 1, r = fma (q, y, x);
+      if (ISNAN (x) || ISNAN (y))
+        return x + y; /* NaN */
+      else if (isinf (y))
+        return x;
+      else
+        /* x infinite or y zero */
+        return NAN;
     }
-  return r;
 }