[Bug math/14048] New: fmod() incorrectly returns NaN for (some?) denormalized inputs

["jkummerow at chromium dot org" <sourceware-bugzilla at sourceware dot org>]

http://sourceware.org/bugzilla/show_bug.cgi?id=14048

             Bug #: 14048
           Summary: fmod() incorrectly returns NaN for (some?)
                    denormalized inputs
           Product: glibc
           Version: 2.15
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: unassigned@sourceware.org
        ReportedBy: jkummerow@chromium.org
    Classification: Unclassified


I've come across a case where fmod() does not return the expected result.
Reduced repro:
------
#include <math.h>
#include <stdio.h>
#include <iostream>

int main() {
  double x = 2.225073858507201e-308;
  double y = 5e-324;
  double z = fmod(x, y);
  // printf("result: %g\n", z);  // see [1] below.
  std::cout << "result: " << z << std::endl;
  return 0;
}
------
Expected result: 0
Actual result: -nan

I can repro the bug on both a Gentoo (gcc-4.5.3, kernel 3.3.4) and an Ubuntu
Precise (gcc-4.6.3, kernel 3.2.5) system, which both have glibc-2.15, and both
are x86_64. It works correctly on Ubuntu Lucid (glibc-2.11, gcc-4.4.2, kernel
2.6.38.8).

Further, it works correctly when compiling with either of -O1, -O2, -O3. It
also works correctly when removing the "std::cout << ..." and "#include
<iostream>" lines, and using the "printf(..." line (marked [1]) instead.

I've looked at the generated machine code. In the buggy case, glibc's fmod is
called directly:
main():
  400744:    55                       push   rbp
  400745:    48 89 e5                 mov    rbp,rsp
  400748:    48 83 ec 20              sub    rsp,0x20
  40074c:    48 b8 ff ff ff ff ff     movabs rax,0xfffffffffffff
  400753:    ff 0f 00 
  400756:    48 89 45 f8              mov    QWORD PTR [rbp-0x8],rax
  40075a:    b8 01 00 00 00           mov    eax,0x1
  40075f:    48 89 45 f0              mov    QWORD PTR [rbp-0x10],rax
  400763:    f2 0f 10 4d f0           movsd  xmm1,QWORD PTR [rbp-0x10]
  400768:    f2 0f 10 45 f8           movsd  xmm0,QWORD PTR [rbp-0x8]
  40076d:    e8 de fe ff ff           call   400650 <fmod@plt>
  400772:    f2 0f 11 45 e8           movsd  QWORD PTR [rbp-0x18],xmm0
  400777:    f2 0f 10 45 e8           movsd  xmm0,QWORD PTR [rbp-0x18]
  40077c:    bf dc 08 40 00           mov    edi,0x4008dc
  400781:    b8 01 00 00 00           mov    eax,0x1
  400786:    e8 75 fe ff ff           call   400600 <printf@plt>
  40078b:    b8 00 00 00 00           mov    eax,0x0
  400790:    c9                       leave  
  400791:    c3                       ret


When I do any of the changes that make it work (e.g. remove the iostream
include), g++ decides to inline an FPU-based implementation of fmod (which
seems to work as expected) and only calls out to glibc as a fallback:
main():
  400604:    55                       push   rbp
  400605:    48 89 e5                 mov    rbp,rsp
  400608:    48 83 ec 40              sub    rsp,0x40
  40060c:    48 b8 ff ff ff ff ff     movabs rax,0xfffffffffffff
  400613:    ff 0f 00 
  400616:    48 89 45 f8              mov    QWORD PTR [rbp-0x8],rax
  40061a:    b8 01 00 00 00           mov    eax,0x1
  40061f:    48 89 45 f0              mov    QWORD PTR [rbp-0x10],rax
  400623:    dd 45 f8                 fld    QWORD PTR [rbp-0x8]
  400626:    dd 45 f0                 fld    QWORD PTR [rbp-0x10]
  400629:    d9 c0                    fld    st(0)
  40062b:    d9 c2                    fld    st(2)
  40062d:    d9 f8                    fprem  
  40062f:    df e0                    fnstsw ax
  400631:    f6 c4 04                 test   ah,0x4
  400634:    75 f7                    jne    40062d <main+0x29>
  400636:    dd d9                    fstp   st(1)
  400638:    dd 5d d8                 fstp   QWORD PTR [rbp-0x28]
  40063b:    f2 0f 10 45 d8           movsd  xmm0,QWORD PTR [rbp-0x28]
  400640:    66 0f 2e c0              ucomisd xmm0,xmm0
  400644:    7a 06                    jp     40064c <main+0x48>
  400646:    66 0f 2e c0              ucomisd xmm0,xmm0
  40064a:    74 17                    je     400663 <main+0x5f>
  40064c:    dd 5d c8                 fstp   QWORD PTR [rbp-0x38]
  40064f:    f2 0f 10 4d c8           movsd  xmm1,QWORD PTR [rbp-0x38]
  400654:    dd 5d c8                 fstp   QWORD PTR [rbp-0x38]
  400657:    f2 0f 10 45 c8           movsd  xmm0,QWORD PTR [rbp-0x38]
  40065c:    e8 af fe ff ff           call   400510 <fmod@plt>
  400661:    eb 04                    jmp    400667 <main+0x63>
  400663:    dd d8                    fstp   st(0)
  400665:    dd d8                    fstp   st(0)
  400667:    f2 0f 11 45 e8           movsd  QWORD PTR [rbp-0x18],xmm0
  40066c:    f2 0f 10 45 e8           movsd  xmm0,QWORD PTR [rbp-0x18]
  400671:    bf 7c 07 40 00           mov    edi,0x40077c
  400676:    b8 01 00 00 00           mov    eax,0x1
  40067b:    e8 70 fe ff ff           call   4004f0 <printf@plt>
  400680:    b8 00 00 00 00           mov    eax,0x0
  400685:    c9                       leave  
  400686:    c3                       ret


So, it looks to me like glibc's fmod() doesn't handle this case (max
denormalized double modulo min denormalized double) correctly. Am I barking up
the wrong tree?

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.