diff options
author | Bruno Haible <bruno@clisp.org> | 2012-02-26 00:36:01 +0100 |
---|---|---|
committer | Bruno Haible <bruno@clisp.org> | 2012-02-26 00:36:01 +0100 |
commit | 0e1c6ff93f27c939ba9e0df945b16ef98eaaeef1 (patch) | |
tree | ef99348fd84d3144205514baa0e384464cfc9e8a /lib/remainderl.c | |
parent | 87fdcafc0dd23b672f8d9496625b6c04dce1e324 (diff) | |
download | gnulib-0e1c6ff93f27c939ba9e0df945b16ef98eaaeef1.tar.gz |
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.
Diffstat (limited to 'lib/remainderl.c')
-rw-r--r-- | lib/remainderl.c | 23 |
1 files changed, 21 insertions, 2 deletions
diff --git a/lib/remainderl.c b/lib/remainderl.c index b43592a644..3476f948e2 100644 --- a/lib/remainderl.c +++ b/lib/remainderl.c @@ -32,8 +32,27 @@ remainderl (long double x, long double y) 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 |