Wrong math results w/ Optim. [-O3 -mpentium]

[Kurt Garloff <garloff at kg1 dot ping dot de>]

Hi,

here's a bug of the g++ compilers:
First: egcs-2.91.14 (compiled w/ haifa, threads) glibc-2.0.6, Lnx.2190

garloff@kg1:/home/garloff/C > eg++ -O3 bug_short.cc; a.out
gamma (0.333333) = 2.67894
chi_s (5,0.666667) = 0.00523951 (k=0.186641)
0.00523951

Here are correct results:

garloff@kg1:/home/garloff/C > eg++ -O2 bug_short.cc; a.out
gamma (0.333333) = 2.67894
chi_s (5,0.666667) = 0.00523951 (k=0.186641)
0.00523951
garloff@kg1:/home/garloff/C > eg++ -O3 bug_short.cc -DNOBUG; a.out
chi_s (0.666667,5)
gamma (2.5) = 1.32934
chi_s (0.666667,5) = 0.0518665 (k=0.132981)
0.0518665
garloff@kg1:/home/garloff/C > eg++ -O3 bug_short.cc -m486; a.out
gamma (2.5) = 1.32934
chi_s (0.666667,5) = 0.0518665 (k=0.132981)
0.0518665

g++-2.7.2.3 suceeds, g++-2.7.2.3 -O fails, -DNOBUG doesn't help.
g++-2.8.0 suceeds, g++-2.8-0 -O fails, ditto
egcs-2.90.23 -O3 -m486 suceeds, egcs-2.90.23 -O3 -DNOBUG fails

The arguments get overwritten in some manner or the gamma() is wrong.

How can it be that they all fail ?
Should I recompile egcs, glibc or buy a new CPU (6x86-P200+) ?

-- 
Kurt Garloff, Dortmund 
<K.Garloff@ping.de>
PGP key on http://student.physik.uni-dortmund.de/homepages/garloff

/* bug_short.cc */

#include <math.h>
#include <iostream>

const double pi = 3.14159265358979323846264338328;
const double E1 = 2.71828182845904523536028747135;


// Stirling is only fine for large numbers
double fact (const double x)
{
	double corr = 1.0 + 1.0/(x*12) + 1.0/(x*x*288) - 139.0/(x*x*x*51850);
	return (pow(x/E1,x) * sqrt(x*pi*2.0) * corr);
};

// the gamma fct
double gamma (const double x)
{
	double res; double x0 = x-1;
	if (x0 > 32.0) return fact (x0);
	x0 += 32.0 - floor(x0); // [32.0;33.0[
	res = fact (x0);
        while ((x-x0) < 0.5) { res /= x0; x0 -= 1.0; };
        cout << "gamma (" << x << ") = " << res << endl; return res;
};

double chi_s (const double x, const double n)
{
#ifdef NOBUG
	cout << "chi_s (" << x << "," << n << ")\n";
#endif
	double k = 1.0 / (pow(2.0,n/2.0) * gamma(n/2.0));
	double res = k * pow(x,(n-2.0)/2.0) * exp (-x/2.0);
	cout << "chi_s (" << x << "," << n << ") = " << res 
		<< " (k=" << k << ")" << endl; return res;
};

int main ()
{
	cout << chi_s(2.0/3,5) << endl;
};