optimization/5017: i386 floating point optimization

[ghh at mpq dot mpg dot de]

>Number:         5017
>Category:       optimization
>Synopsis:       i386 floating point optimization
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    unassigned
>State:          open
>Class:          wrong-code
>Submitter-Id:   net
>Arrival-Date:   Wed Dec 05 04:46:00 PST 2001
>Closed-Date:
>Last-Modified:
>Originator:     Gerhard Heinzel
>Release:        3.0.2
>Organization:
Albert Einstein Institut Hannover, Germany
>Environment:
System: Linux lp27 2.4.4-4GB #1 Wed May 16 00:37:55 GMT 2001 i686 unknown
Architecture: i686

	
host: i686-pc-linux-gnu
build: i686-pc-linux-gnu
target: i686-pc-linux-gnu
configured with: ../configure 
>Description:
The included small test program fails when compiled with
gcc -O3 -o mybes mybes.c -lm
It works with -O0.
indication of failure is the following output.
The first line starting with "5.00" is garbage.
There should not be any "nan"

4.00  18      4.0944813915730677499e-11          11.301921952136330496
4.00  19      4.3099804121821765789e-12          11.301921952136330496
4.00  20      4.3099804121821765788e-13          11.301921952136330496
4.00  21      4.1047432496973110274e-14          11.301921952136330496
4.00  22      3.7315847724521009341e-15          11.301921952136330496
4 11.301921952136330773 6.62406158447265625
5.00   1                            nan                              1
5.00   2                            nan                            nan
5.00   3                            nan                            nan
5.00   4                            nan                            nan
5.00   5                            nan                            nan

>How-To-Repeat:

#include <stdio.h>
double
mybesi0 (double x)
{
  int k = 0;
  long double sum = 0;
  long double add = 1;
  long double x2 = x / 2.;

  if (x == 0)
    return 1;
  if (x < 0)
    x = -x;
  if (x > 50)
    {
      fprintf (stderr, "argument to mybesi0 too large:%g\n", x);
      exit (1);
    }
  add = 1;
  for (;;)
    {
      sum += add * add;
      add /= (++k);
      add *= x2;
      printf("%4.2f %3d %30.20Lg %30.20Lg\n",x,k,add,sum);
      if (add / sum < 1e-15)
	break;
    }
  return sum;
}

float
bessi0 (float x)
{
  float ax, ans;
  double y;

  if ((ax = fabs (x)) < 3.75)
    {
      y = x / 3.75;
      y *= y;
      ans = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
		  + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
    }
  else
    {
      y = 3.75 / ax;
      ans = (exp (ax) / sqrt (ax)) * (0.39894228 + y * (0.1328592e-1
		   + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
	       + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1
						  + y * 0.392377e-2))))))));
    }
  return ans;
}

void
main (void)
{
  double x;
  for (x = 0; x <= 50; x += 1)
   printf ("%.20g %.20g %.20g\n", x, mybesi0 (x), bessi0(x)); 
}


>Fix:
>Release-Note:
>Audit-Trail:
>Unformatted: