[Bug math/23395] fmod repeated subtraction algorithm introduces round off errors

["adhemerval.zanella at linaro dot org" <sourceware-bugzilla at sourceware dot org>]

https://sourceware.org/bugzilla/show_bug.cgi?id=23395

--- Comment #5 from Adhemerval Zanella <adhemerval.zanella at linaro dot org> ---
(In reply to Irv Lustig from comment #4)
> The reason that I think it is not exact is because the result of 
> fmod(x, x/100.0) should be close to 0.
> 
> For example, fmod(1.0, 0.01) in infinite precision should be zero, since
> 0.01 divides evenly into 1.0.  You can argue that it is a precision problem,
> but then the results should get closer to zero as you increase the
> precision. Yet, glibc fmod(1.0, 0.01) is returning  9.99999999999998e-03,
> and so is mpfr, which I also think is buggy.

I think you are still missing the point here.  First, fmod is defined to
rounded toward zero to an integer, so you can't use round (which rounds away
from zero). Also, since you are using floating-point operations you can't rely
on system floating rounding mode, you need to actually set/restore it in your
function:

---
double rndfmod(double x, double y)
{
  int r = fegetround ();
  fesetround (FE_TOWARDZERO);
  double result = x - trunc(x/y)*y;
  fesetround (r);
  return result;
}
---

With this code you get:

fmod(1.0000e+01, 1.0000e-01) = 9.99999999999995e-02 
rndfmod 9.99999999999996e-02 
mpfr 9.99999999999995e-02

And then you get precision issues by using a non-default rounding mode: the
operations requires to calculate 'x - trunc(x/y)*y' without intermediary
results to get the correct exact answer (similar to fma/sqrt operation). When
you use multiple FP operations you eventually get precision issues.

Also, you need to take care of subnormal numbers (since C standard says fmod
should be exact as well) and its handling on different architectures (for
instance some architecture flush them to zero).

Also stating MPFR is wrong for this operation is quite a bold statement (it is
used on constant folding for GCC for instance and other multiprecision floating
point operations). 

> 
> The Stackoverflow reference states:
> 
> The IEEE Standard 754 defines the remainder operation x REM y as the
> mathematical operation x - (round(x/y)*y).
> 
> The function in my latest attachment new_fmod() does that computation, and
> it is almost always zero.

-- 
You are receiving this mail because:
You are on the CC list for the bug.