[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: new module 'fmodl'
From: |
Bruno Haible |
Subject: |
Re: new module 'fmodl' |
Date: |
Sun, 26 Feb 2012 00:37:25 +0100 |
User-agent: |
KMail/4.7.4 (Linux/3.1.0-1.2-desktop; KDE/4.7.4; x86_64; ; ) |
> long double
> fmodl (long double x, long double y)
> {
> long double i = truncl (x / y);
> return fmal (- i, y, x);
> }
Oops, 1 bug in 2 lines of code, is a bad ratio. I was not thinking at the
effects of rounding errors. This should fix it:
2012-02-25 Bruno Haible <address@hidden>
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.
--- lib/fmodl.c.orig Sun Feb 26 00:31:10 2012
+++ lib/fmodl.c Sat Feb 25 23:28:43 2012
@@ -32,8 +32,48 @@
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
--- lib/remainder.c.orig Sun Feb 26 00:31:10 2012
+++ lib/remainder.c Sat Feb 25 23:31:24 2012
@@ -22,6 +22,25 @@
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;
}
--- lib/remainderf.c.orig Sun Feb 26 00:31:10 2012
+++ lib/remainderf.c Sat Feb 25 23:31:24 2012
@@ -25,7 +25,26 @@
#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
}
--- lib/remainderl.c.orig Sun Feb 26 00:31:10 2012
+++ lib/remainderl.c Sat Feb 25 23:31:24 2012
@@ -32,8 +32,27 @@
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