Comment 1 for bug 1000498

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?