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;
}