The implementation of floor() and ceil() contained in sysdeps/i386/fpu/bits/mathinline.h can cause incorrect code to be generated. This problem occurs when the optimizer decides to move surrounding floating- point instructions into the region between the two "fldcw" instructions: fldcw __cwtmp frndint fldcw __cw Whenever this happens, these instructions are executed with the rounding mode forced to "round down" (floor) or "round down" (ceil), instead of the proper rounding mode in effect for the general computation. I have short C and C++ programs that demonstrate the bug. The programs get different answers when compiled with -O2 and without. Compiling with -O2 -S and examining the generated assembly code shows TWO floating-point operations moved into the space between the first fldcw and the frndint instruction. The floating-point constants used in the program are such that the differences result directly from the different rounding mode for these two instructions, not from any 53-bit versus 64-bit mantissa issues. I can send you these example programs if you decide that really need to see gcc doing this. Here is a patch that corrects the problem by generating all three of these instructions in a single "asm" directive so that the optimizer won't move any other instructions into the middle of this sequence. In my humble opinion, the code correctness obtained via the patch outweighs any possible advantages that the original implementation might have had -- for example, opening up the 3-instruction sequence for instruction scheduling, etc. Getting the wrong answer at higher speed still gets you the wrong answer. ;^> *** glibc-2.5/sysdeps/i386/fpu/bits/mathinline.h 2004-09-07 18:23:42.000000000 -0400 --- glibc-2.5.fixed/sysdeps/i386/fpu/bits/mathinline.h 2006-11-02 16:04:43.000000000 -0500 *************** *** 526,559 **** register long double __exm1 = __expm1l (-__fabsl (__x + __x)); \ return __exm1 / (__exm1 + 2.0) * __sgn1l (-__x)) #endif __inline_mathcodeNP (floor, __x, \ register long double __value; \ __volatile unsigned short int __cw; \ __volatile unsigned short int __cwtmp; \ __asm __volatile ("fnstcw %0" : "=m" (__cw)); \ __cwtmp = (__cw & 0xf3ff) | 0x0400; /* rounding down */ \ ! __asm __volatile ("fldcw %0" : : "m" (__cwtmp)); \ ! __asm __volatile ("frndint" : "=t" (__value) : "0" (__x)); \ ! __asm __volatile ("fldcw %0" : : "m" (__cw)); \ return __value) __inline_mathcodeNP (ceil, __x, \ register long double __value; \ __volatile unsigned short int __cw; \ __volatile unsigned short int __cwtmp; \ __asm __volatile ("fnstcw %0" : "=m" (__cw)); \ __cwtmp = (__cw & 0xf3ff) | 0x0800; /* rounding up */ \ ! __asm __volatile ("fldcw %0" : : "m" (__cwtmp)); \ ! __asm __volatile ("frndint" : "=t" (__value) : "0" (__x)); \ ! __asm __volatile ("fldcw %0" : : "m" (__cw)); \ return __value) #ifdef __FAST_MATH__ # define __ldexp_code \ register long double __value; \ __asm __volatile__ \ ("fscale" \ : "=t" (__value) : "0" (__x), "u" ((long double) __y)); \ return __value --- 526,559 ---- register long double __exm1 = __expm1l (-__fabsl (__x + __x)); \ return __exm1 / (__exm1 + 2.0) * __sgn1l (-__x)) #endif __inline_mathcodeNP (floor, __x, \ register long double __value; \ __volatile unsigned short int __cw; \ __volatile unsigned short int __cwtmp; \ __asm __volatile ("fnstcw %0" : "=m" (__cw)); \ __cwtmp = (__cw & 0xf3ff) | 0x0400; /* rounding down */ \ ! __asm __volatile ("fldcw %2\n\tfrndint\n\tfldcw %3" \ ! : "=t" (__value) \ ! : "0" (__x), "m" (__cwtmp), "m" (__cw)); \ return __value) __inline_mathcodeNP (ceil, __x, \ register long double __value; \ __volatile unsigned short int __cw; \ __volatile unsigned short int __cwtmp; \ __asm __volatile ("fnstcw %0" : "=m" (__cw)); \ __cwtmp = (__cw & 0xf3ff) | 0x0800; /* rounding up */ \ ! __asm __volatile ("fldcw %2\n\tfrndint\n\tfldcw %3" \ ! : "=t" (__value) \ ! : "0" (__x), "m" (__cwtmp), "m" (__cw)); \ return __value) #ifdef __FAST_MATH__ # define __ldexp_code \ register long double __value; \ __asm __volatile__ \ ("fscale" \ : "=t" (__value) : "0" (__x), "u" ((long double) __y)); \ return __value
Created attachment 1398 [details] Proposed patch to fix this bug. Sorry, here is the patch once again as an attachment (in case the original cut/pasted one got word-wrapped, or something).
The patch is not sufficient. If the compiler can pull something in between the fldcw instructions then code which changes the CW might be pulled between the fnstcw and first fldcw. I've checked in a complete patch. And next time, don't use this unreadable context diff format. Alsways use -u.
Subject: Bug 3451 CVSROOT: /cvs/glibc Module name: libc Branch: glibc-2_5-branch Changes by: jakub@sourceware.org 2007-01-12 15:31:04 Modified files: . : ChangeLog sysdeps/i386/fpu/bits: mathinline.h Log message: [BZ #3451] * sysdeps/i386/fpu/bits/mathinline.h (floor): Make rounding mode change atomic. (ceil): Likewise. Patches: http://sourceware.org/cgi-bin/cvsweb.cgi/libc/ChangeLog.diff?cvsroot=glibc&only_with_tag=glibc-2_5-branch&r1=1.10362.2.11&r2=1.10362.2.12 http://sourceware.org/cgi-bin/cvsweb.cgi/libc/sysdeps/i386/fpu/bits/mathinline.h.diff?cvsroot=glibc&only_with_tag=glibc-2_5-branch&r1=1.58&r2=1.58.8.1