More wrong math results w/ inline

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

Here is an example for more wrong math results:

garloff@kg1:/home/garloff/C > egcc -O2 -o bug_math bug_math.c -lm
garloff@kg1:/home/garloff/C > bug_math 42 12 288 372.94964; bug_math
5.7718313756965, 4.885604e+00
5.7718313755102, 4.885604e+00

Here are correct results:

garloff@kg1:/home/garloff/C > egcc -o bug_math bug_math.c -lm
garloff@kg1:/home/garloff/C > bug_math 42 12 288 372.94964; bug_math
0.8862269255100, 5.720524e-11
0.8862269254528, 2.109424e-15
garloff@kg1:/home/garloff/C > egcc -O2 -o bug_math bug_math.c -lm -fno-inline
garloff@kg1:/home/garloff/C > bug_math 42 12 288 372.94964; bug_math
0.8862269255100, 5.720535e-11
0.8862269254528, 2.109424e-15

egcs/gcc-28/gcc-27 all behave the same. Ughh !
(Machine: 6x86, Linux-2190, glibc-206, egcs-980315)

Inlining of math functions seem to produce this. Is it a problem of the
small ix86 FPU stack? -ffloat-store fixes the problem as well as -O3 or
-fno-inline.

Sourcecode is appended. (It is an implementation for the gamma function.)

I hope this can be fixed soon. Otherwise I have to recompile everything w/o
inlining to be sure to get the correct math results from my numerical
simulation code. I fear this would drop performance by a factor of 3.

-- 
Kurt Garloff, Dortmund 
<K.Garloff@ping.de>
PGP key on http://student.physik.uni-dortmund.de/homepages/garloff
/* bug_math.c: gcc math inlining bug */
/* calc gamma fct using stirling formula */

#include <math.h>
#include <stdlib.h>

const double pi = 3.1415926535897932384626433832795028841968;
const double E1 = 2.71828182845904523536028747135;

double cut;

// Stirling is only fine for large numbers
double stir (const double x, const double s0, const double s1, const double s2)
{
	  double r = 1.0/x;
	  //double corr = 1.0 + r/12 + r/(x*288) - 139.0*r/(x*x*51840); //Abramov/Stegun
	  double corr = 1.0 + r/s0 + r/(x*s1) - r/(x*x*s2);
	  return (pow(x/E1,x) * sqrt(x*pi*2.0) * corr);
};


double fct (const double x, const double s0, const double s1, const double s2)
{
	double x0, res;
	if (x >= cut) return stir (x-1.0, s0, s1, s2);
	x0 = x + cut - 2*floor (x/2) - 1.0;
	res = stir (x0, s0, s1, s2);
	while ((x-x0) < 0.5) { res /= x0; x0 -= 1.0; }
	return res;
};

			
int main (int argc, char* argv[])
{
	const double arg = 21.0; double s0, s1, s2; char *dum; 
	if (argc > 1) cut = strtod (argv[1], &dum); else cut = 42;
	if (argc > 2) s0 = strtod (argv[2], &dum); else s0 = 12.00000034;
	if (argc > 3) s1 = strtod (argv[3], &dum); else s1 = 287.97452028;
	if (argc > 4) s2 = strtod (argv[4], &dum); else s2 = 371.045958768;
	printf ("%15.13f, %e\n", fct(1.5,s0,s1,s2), fct(1.5,s0,s1,s2)-0.5*sqrt(pi));
};