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

["robert.dodier at gmail dot com" <sourceware-bugzilla at sourceware dot org>]

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

--- Comment #3 from Robert Dodier <robert.dodier at gmail dot com> ---
Thanks for the hints. With the flags -msse -mfpmath=sse -fno-builtin, I am
still seeing the unexpected result (namely that the final output, labeled "x -
(floor(x/y)*y + fmod(x,y))" below), is pretty big; I was expecting zero or a
small multiple of the floating point epsilon.

The test program is now working with specific bit patterns for its inputs.

$ ./test-fmod-1 
gnu_get_libc_version => 2.23
x = 0x4020000000000000 = 8.0000000000000000e+00, y = 0x3ff999999999999a =
1.6000000000000001e+00
x/y = 8.0000000000000000e+00/1.6000000000000001e+00 = 5.0000000000000000e+00
floor(x/y) = floor(8.0000000000000000e+00/1.6000000000000001e+00) =
5.0000000000000000e+00
floor(x/y)*y =
floor(8.0000000000000000e+00/1.6000000000000001e+00)*1.6000000000000001e+00 =
8.0000000000000000e+00
x - floor(x/y)*y = 8.0000000000000000e+00 -
floor(8.0000000000000000e+00/1.6000000000000001e+00)*1.6000000000000001e+00 =
-4.4408920985006262e-16
fmod(x,y) = fmod(8.0000000000000000e+00, 1.6000000000000001e+00) =
1.5999999999999996e+00
x - (floor(x/y)*y + fmod(x,y)) = 8.0000000000000000e+00 -
(floor(8.0000000000000000e+00/1.6000000000000001e+00)*1.6000000000000001e+00 +
fmod(8.0000000000000000e+00, 1.6000000000000001e+00)) = -1.6000000000000001e+00

Here is the current version of the test program:

$ cat test-fmod-1.c 
#include <stdio.h>
#include <math.h>
#include <gnu/libc-version.h>

unsigned char xbytes [] = { 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x20, 0x40 };
// 0x4020000000000000 = 8.0e0
unsigned char ybytes [] = { 0x9A, 0x99, 0x99, 0x99, 0x99, 0x99, 0xF9, 0x3F };
// 0x3FF999999999999A = 1.6e0

int main (int argc, char** argv)
{
    double x = *(double *)xbytes, y = *(double *)ybytes;
    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 ("x = 0x%02x%02x%02x%02x%02x%02x%02x%02x = %.16le, y =
0x%02x%02x%02x%02x%02x%02x%02x%02x = %.16le\n",
        xbytes[7], xbytes[6], xbytes[5], xbytes[4], xbytes[3], xbytes[2],
xbytes[1], xbytes[0], x,
        ybytes[7], ybytes[6], ybytes[5], ybytes[4], ybytes[3], ybytes[2],
ybytes[1], ybytes[0], y);
    printf ("x/y = %.16le/%.16le = %.16le\n", x, y, xy);
    printf ("floor(x/y) = floor(%.16le/%.16le) = %.16le\n", x, y, floor_xy);
    printf ("floor(x/y)*y = floor(%.16le/%.16le)*%.16le = %.16le\n", x, y, y,
floor_xy*y);
    printf ("x - floor(x/y)*y = %.16le - floor(%.16le/%.16le)*%.16le =
%.16le\n", x, x, y, y, x - floor_xy*y);
    printf ("fmod(x,y) = fmod(%.16le, %.16le) = %.16le\n", x, y, fmod_xy);
    printf ("x - (floor(x/y)*y + fmod(x,y)) = %.16le -
(floor(%.16le/%.16le)*%.16le + fmod(%.16le, %.16le)) = %.16le\n",
        x, x, y, y, x, y, r);

    return 0;
}

Here is how I compiled and linked it:

$ gcc -g -c -msse -mfpmath=sse -fno-builtin test-fmod-1.c
$ gcc -nostdlib -nostartfiles -static -msse -mfpmath=sse -fno-builtin -o
test-fmod-1 /usr/local/glibc-test/lib/crt1.o /usr/local/glibc-test/lib/crti.o
`gcc --print-file-name=crtbegin.o` test-fmod-1.o -Wl,--start-group
/usr/local/glibc-test/lib/libc.a /usr/local/glibc-test/lib/libm.a -lgcc
-lgcc_eh -Wl,--end-group `gcc --print-file-name=crtend.o`
/usr/local/glibc-test/lib/crtn.o

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