[Bug math/19869] New: floor and fmod not consistent with each other

[robert.dodier at gmail dot com]

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

            Bug ID: 19869
           Summary: floor and fmod not consistent with each other
           Product: glibc
           Version: 2.23
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
          Assignee: unassigned at sourceware dot org
          Reporter: robert.dodier at gmail dot com
  Target Milestone: ---

Created attachment 9130
  --> https://sourceware.org/bugzilla/attachment.cgi?id=9130&action=edit
program to find x and y s.t. floor(x/y) and fmod(x,y) are inconsistent

When x is very close to an integer multiple of y, floor and fmod can give
inconsistent results in the sense that floor(x/y)*y + fmod(x,y) is not anywhere
close to x.

For example, with x = 8e0 and y = 1.6e0, fmod(x,y) = 1.5999999999999996e+0 but
floor(x/y) = 5e0, so floor(x/y)*y + fmod(x,y) = 9.6e0 approximately.

The problem is that y is very slightly larger than x/floor(x/y), so
floor(x/y)*y is a little bigger than x. It appears that fmod(x, y) is returned
as y + (x - floor(x/y)*y) so fmod(x, y) is positive (it is a little smaller
than y). But floor(x/y) is not likewise adjusted. 

It seems like either fmod(x,y) should not be adjusted, and returned as a small
negative number (the choice of some other libraries) or zero (the choice of
still others), or else floor(x/y) should be reduced by 1, and fmod(x,y)
returned as it is now. There are probably other possibilities. I won't try to
sort out the arguments for different options.

There are many other examples; see attached program test-fmod.c to find them.

I am working with glibc 2.23. uname -a says: Linux freekbox 3.13.0-44-generic
#73-Ubuntu SMP Tue Dec 16 00:23:46 UTC 2014 i686 i686 i686 GNU/Linux

I only tried it with doubles.

PS. Here is some detailed output.

$ cat test-fmod-1.log 
gnu_get_libc_version => 2.23
8.0000000000000000e+00/1.6000000000000001e+00 = 5.0000000000000000e+00
floor(8.0000000000000000e+00/1.6000000000000001e+00) = 5.0000000000000000e+00
floor(8.0000000000000000e+00/1.6000000000000001e+00)*1.6000000000000001e+00 =
8.0000000000000000e+00
8.0000000000000000e+00 -
floor(8.0000000000000000e+00/1.6000000000000001e+00)*1.6000000000000001e+00 =
-4.4408920985006262e-16
fmod(8.0000000000000000e+00, 1.6000000000000001e+00) = 1.5999999999999996e+00
8.0000000000000000e+00 -
(floor(8.0000000000000000e+00/1.6000000000000001e+00)*1.6000000000000001e+00 +
fmod(8.0000000000000000e+00, 1.6000000000000001e+00)) = -1.6000000000000001e+00

PPS. Here's the program to generate the above output.

#include <stdio.h>
#include <math.h>
#include <gnu/libc-version.h>

int main (int argc, char** argv)
{
    double x = 8e0, y = 1.6e0;
    double xy, floor_xy, fmod_xy, r;

    printf ("gnu_get_libc_version => %s\n", gnu_get_libc_version ());

    xy = x/y;
    floor_xy = floor(xy);
    fmod_xy = fmod(x, y);

    r = x - (floor_xy*y + fmod_xy);

    printf ("%.16le/%.16le = %.16le\n", x, y, xy);
    printf ("floor(%.16le/%.16le) = %.16le\n", x, y, floor_xy);
    printf ("floor(%.16le/%.16le)*%.16le = %.16le\n", x, y, y, floor_xy*y);
    printf ("%.16le - floor(%.16le/%.16le)*%.16le = %.16le\n", x, x, y, y, x -
floor_xy*y);
    printf ("fmod(%.16le, %.16le) = %.16le\n", x, y, fmod_xy);
    printf ("%.16le - (floor(%.16le/%.16le)*%.16le + fmod(%.16le, %.16le)) =
%.16le\n", x, x, y, y, x, y, r);

    return 0;
}

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