Mercurial > hg > octave-lojdl > gnulib-hg
changeset 16418:d10d4da19813
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.
author | Bruno Haible <bruno@clisp.org> |
---|---|
date | Sun, 26 Feb 2012 00:36:01 +0100 |
parents | 861cd2de5c55 |
children | 271a263949d9 |
files | ChangeLog lib/fmodl.c lib/remainder.c lib/remainderf.c lib/remainderl.c |
diffstat | 5 files changed, 114 insertions(+), 8 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,12 @@ +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'.
--- a/lib/fmodl.c +++ b/lib/fmodl.c @@ -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
--- a/lib/remainder.c +++ b/lib/remainder.c @@ -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; }
--- a/lib/remainderf.c +++ b/lib/remainderf.c @@ -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 }
--- a/lib/remainderl.c +++ b/lib/remainderl.c @@ -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