+2012-02-25 Bruno Haible <bruno@clisp.org>
+
+ fmodl, remainder*: Avoid wrong results due to rounding errors.
+ * lib/fmodl.c (fmodl): Correct the result if it is not within the
+ expected bounds.
+ * lib/remainderf.c (remainderf): Likewise.
+ * lib/remainder.c (remainder): Likewise.
+ * lib/remainderl.c (remainderl): Likewise.
+
2012-02-25 Bruno Haible <bruno@clisp.org>
Tests for module 'remainderl'.
long double
fmodl (long double x, long double y)
{
- long double i = truncl (x / y);
- return fmal (- i, y, x);
+ long double q = - truncl (x / y);
+ long double r = fmal (q, y, x); /* = x + q * y, computed in one step */
+ /* Correct possible rounding errors in the quotient x / y. */
+ if (y >= 0)
+ {
+ if (x >= 0)
+ {
+ /* Expect 0 <= r < y. */
+ if (r < 0)
+ q += 1, r = fmal (q, y, x);
+ else if (r >= y)
+ q -= 1, r = fmal (q, y, x);
+ }
+ else
+ {
+ /* Expect - y < r <= 0. */
+ if (r > 0)
+ q -= 1, r = fmal (q, y, x);
+ else if (r <= - y)
+ q += 1, r = fmal (q, y, x);
+ }
+ }
+ else
+ {
+ if (x >= 0)
+ {
+ /* Expect 0 <= r < - y. */
+ if (r < 0)
+ q -= 1, r = fmal (q, y, x);
+ else if (r >= - y)
+ q += 1, r = fmal (q, y, x);
+ }
+ else
+ {
+ /* Expect y < r <= 0. */
+ if (r > 0)
+ q += 1, r = fmal (q, y, x);
+ else if (r <= y)
+ q -= 1, r = fmal (q, y, x);
+ }
+ }
+ return r;
}
#endif
double
remainder (double x, double y)
{
- double i = round (x / y);
- return fma (- i, y, x);
+ 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)
+ {
+ /* 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);
+ }
+ 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);
+ }
+ return r;
}
#if HAVE_REMAINDER
return (float) remainder ((double) x, (double) y);
#else
- float i = roundf (x / y);
- return fmaf (- i, y, x);
+ float q = - roundf (x / y);
+ float r = fmaf (q, y, x); /* = x + q * y, computed in one step */
+ /* Correct possible rounding errors in the quotient x / y. */
+ float half_y = 0.5L * y;
+ if (y >= 0)
+ {
+ /* Expect -y/2 <= r <= y/2. */
+ if (r > half_y)
+ q -= 1, r = fmaf (q, y, x);
+ else if (r < - half_y)
+ q += 1, r = fmaf (q, y, x);
+ }
+ else
+ {
+ /* Expect y/2 <= r <= -y/2. */
+ if (r > - half_y)
+ q += 1, r = fmaf (q, y, x);
+ else if (r < half_y)
+ q -= 1, r = fmaf (q, y, x);
+ }
+ return r;
#endif
}
long double
remainderl (long double x, long double y)
{
- long double i = roundl (x / y);
- return fmal (- i, y, x);
+ long double q = - roundl (x / y);
+ long double r = fmal (q, y, x); /* = x + q * y, computed in one step */
+ /* Correct possible rounding errors in the quotient x / y. */
+ long double half_y = 0.5L * y;
+ if (y >= 0)
+ {
+ /* Expect -y/2 <= r <= y/2. */
+ if (r > half_y)
+ q -= 1, r = fmal (q, y, x);
+ else if (r < - half_y)
+ q += 1, r = fmal (q, y, x);
+ }
+ else
+ {
+ /* Expect y/2 <= r <= -y/2. */
+ if (r > - half_y)
+ q += 1, r = fmal (q, y, x);
+ else if (r < half_y)
+ q -= 1, r = fmal (q, y, x);
+ }
+ return r;
}
#endif