An x86 bug in egcs 1.0.2?

[H.J. Lu]

Can someone please tell me if it is an egcs 1.0.2 bug?

Thanks.


-- 
H.J. Lu (hjl@gnu.org)
---
/*
Return-Path: <bregor@aegis.anu.edu.au>
Date: Fri, 10 Apr 1998 14:12:36 +1000 (EST)
Sender: bregor@aegis.anu.edu.au
From: "Roger W. Brown" <bregor@anusf.anu.edu.au>
To: hjl@gnu.org
Subject: Re: egcs-1.0.2 i686-pc-linux-gnu math optimization pb
Status: RO

> 
> I have some c++ code that produces incorrect results when compiled
> with -O2 or -O3 (but not -O).  It is floating-point math.  Have cases
> like this been reported before?
> 
> This is egcs-1.0.2 i686-pc-linux-gnu glibc-2.0.6.
> 

Can you post a small test case?

-- 
H.J. Lu (hjl@gnu.org)

=======================================================================

  H.J.

       Did you see my example (posted in March egcs-bugs)

  No problem under glibc-2.1  Trouble with glibc-2.0.6 & glibc-2.0.7pre2
  The problen is not numerical precision (Jim Wilson), but in-lining asm.

  gcc -o tst -O2 -mpentiumpro tst.c     Works
  gcc -o tst -O2 -mpentium    tst.c     Fails
  gcc -o tst -O2 -m486        tst.c       "

  Regards,

          Roger Brown

==================================================================
*/
extern   double floor  (double __x);
extern   double exp  (double __x);
extern   double lgamma   (double __x);
extern __inline  double
__log2 (double __x)
{
  register double __value;
  __asm __volatile__
    ("fld1\n\t"
     "fxch\n\t"
     "fyl2x"
     : "=t" (__value) : "0" (__x));
  return __value;
}
extern __inline  double
sqrt (double __x)
{
  register double __value;
  __asm __volatile__
    ("fsqrt"
     : "=t" (__value) : "0" (__x));
  return __value;
}
extern __inline  double
pow (double __x, double __y)
{
  register double __value, __exponent;
  long __p = (long) __y;
  if (__x == 0.0 && __y > 0.0)
    return 0.0;
  if (__y == (double) __p)
    {
      double __r = 1.0;
      if (__p == 0)
	return 1.0;
      if (__p < 0)
	{
	  __p = -__p;
	  __x = 1.0 / __x;
	}
      while (1)
	{
	  if (__p & 1)
	    __r *= __x;
	  __p >>= 1;
	  if (__p == 0)
	    return __r;
	  __x *= __x;
	}
    }
  __asm __volatile__
    ("fmul	%%st(1)		# y * log2(x)\n\t"
     "fstl	%%st(1)\n\t"
     "frndint			# int(y * log2(x))\n\t"
     "fxch\n\t"
     "fsub	%%st(1)		# fract(y * log2(x))\n\t"
     "f2xm1			# 2^(fract(y * log2(x))) - 1\n\t"
     : "=t" (__value), "=u" (__exponent) :  "0" (__log2 (__x)), "1" (__y));
  __value += 1.0;
  __asm __volatile__
    ("fscale"
     : "=t" (__value) : "0" (__value), "u" (__exponent));
  return __value;
}

const double pi = 3.14159265358979323846264338328;
const double E1 = 2.71828182845904523536028747135;
double fact (double x)
{
    double corr, t;
    t = 1.0/x;
    corr = 1.0 + t*(1.0/12.0 + t*(1.0/288.0 - t*139.0/51840.0));
    t = pow(x/E1, x);
    t *= sqrt(x*pi*2.0);
    return t*corr;
};
const double cut_off = 46.0;


/* the gamma fct */
static
double gamma (double x)
{
    double res, x0;

    x0 = x-1;
    if (x0 > cut_off) return fact (x0);
    x0 += cut_off - floor(x0);
    res = fact (x0);
    while ((x-x0) < 0.5) { res /= x0; x0 -= 1.0; };
    return res;
};


int main ()
{
       printf(" Gamma(1.2): %16.9e  %16.9e\n", gamma(1.2), exp(lgamma(1.2)));
}