current egcs miscompiles glibc-2.1's inline math

[Zack Weinberg]

The appended source file, when compiled on x86 by current CVS with
optimization on, exits unsuccessfully.  When compiled by 1.1.2, it
exits successfully.  The source was generated by massaging
preprocessor output of this file:

#include <math.h>

int
main(void)
{
  double x;

  x = copysign (0.0, -1.0);  /* -0.0 */
  x = asin(x);
 
  return !isfinite(x);
}

using the math.h from glibc 2.1.  asin(-0.0) is supposed to be -0.0,
but we get QNaN from the code generated by current CVS.

Here's an assembly dump of the relevant chunk:

        fldl .LC0	; .LCO contains double -1.0
        subl $8,%esp
        fstpl (%esp)
        pushl $0
        pushl $0
        call copysign
*
        fmul %st(0),%st
        fld1
        fsubp %st,%st(1)
#APP
        fsqrt
        fpatan
#NO_APP

Why doesn't this work?  Because when we hit the fpatan call, there is
only one value on the stack: 1.0.  fpatan wants two values on the
stack, so it outputs QNaN and sets the `fp-stack underflow' flag in
the FPSR (you can see this by stepping instructions in gdb).

At the position marked by a star, egcs 1.1.2 emits an additional

	fld %st(0)

which gives fpatan the second argument it wants.

The bug appears to be in regstack.  Right before regstack, we have

(set (reg:DF 8 %st(0))
        (call (mem:QI (symbol_ref:SI ("copysign")) 0)
            (const_int 16 [0x10])))

(set (reg:XF 10 %st(2))
        (float_extend:XF (reg/v:DF 8 %st(0))))

(set (reg/v:DF 8 %st(0))
        (mult:DF (reg/v:DF 8 %st(0))
            (reg/v:DF 8 %st(0))))

(set (reg:DF 9 %st(1))
        (const_double:DF (cc0) 0 [0x0] 0 [0x0] 1073709056 [0x3fff8000]))

(set (reg:DF 9 %st(1))
        (minus:DF (reg:DF 9 %st(1))
            (reg/v:DF 8 %st(0))))

(set (reg:XF 8 %st(0))
        (float_extend:XF (reg:DF 9 %st(1))))

(set (reg:XF 8 %st(0))
        (asm_operands/v ("fsqrt") ("=t") 0[ 
                (reg:XF 8 %st(0))
            ] 
            [ 
                (asm_input:XF ("0"))
            ]  ("matht.i") 26))

(set (reg:XF 9 %st(1))
        (reg:XF 10 %st(2)))

(parallel[ 
            (set (reg:XF 8 %st(0))
                (asm_operands/v ("fpatan
        ") ("=t") 0[ 
                        (reg:XF 8 %st(0))
                        (reg:XF 9 %st(1))
                    ] 
                    [ 
                        (asm_input:XF ("0"))
                        (asm_input:XF ("u"))
                    ]  ("matht.i") 18))
            (clobber (reg:QI 9 %st(1)))
        ] )

At this point the compiler still thinks that %st(0) etc are normal
floating point registers.  You can see that it preserves in %st(2) the
result of copysign, and puts it back into %st(1) in time to satisfy
the fpatan asm.

After regstack, the copy to and from %st(2) have simply been deleted.

zw

-- test case --
extern double copysign (double, double);

extern inline int
__finite (double __x)
{
  return (__extension__
          (((((union { double __d; int __i[2]; }) {__d: __x}).__i[1]
             | 0x800fffff) + 1) >> 31));
}

extern __inline long double
__atan2l (long double __y, long double __x)
{
  register long double __value;
  __asm __volatile__ ("fpatan\n\t"
		      : "=t" (__value)
		      : "0" (__x), "u" (__y)
		      : "st(1)");
  return __value;
}

extern __inline long double
__sqrtl (long double __x)
{
  register long double __result;
  __asm __volatile__ ("fsqrt" : "=t" (__result) : "0" (__x));
  return __result;
}

extern __inline double
asin (double __x)
{
  return __atan2l (__x, __sqrtl (1.0 - __x * __x));
}

int
main (void)
{
  double x;

  x = copysign (0.0, -1.0);
  x = asin (x);

  return !__finite (x);
}