[Bug math/3451] New: dangerous inlining of floor and ceil on i386.

["David dot Warme at Group-W-Inc dot com" <sourceware-bugzilla@sourceware.org>]

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

-- 
           Summary: dangerous inlining of floor and ceil on i386.
           Product: glibc
           Version: 2.4
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: aj at suse dot de
        ReportedBy: David dot Warme at Group-W-Inc dot com
                CC: glibc-bugs at sources dot redhat dot com
GCC target triplet: i386-*-*


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

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.