[Bug fortran/33197] Fortran 2008: math functions and other small changes

[fxcoudert at gcc dot gnu dot org]

------- Comment #11 from fxcoudert at gcc dot gnu dot org  2008-03-04 14:26 
-------
On ppc longdouble, convergence requires argument >= 8.6 and, again, less than
100 terms. Why does convergence radius depend on the format of the
floating-point type, I wonder... For the record, here's the small C program I
use for testing:

#include <math.h>
#include <float.h>
#include <stdio.h>

#define PI 3.1415926535897932384626433832795029L

extern long double asympt (long double);

long double asympt (long double x)
{
  long double res, f, y, z;
  int n;

  y = 1 / (2 * x * x);
  z = 1;
  res = 1;
  f = 1;
  for (n = 1; n < 100; n++)
  {
    f *= 2 * n - 1;
    z *= y;
    res += (n % 2 ? -1 : 1) * f * z;
    if (f*z <= DBL_MIN)
      break;
    if (f*z <= (nextafterl (res, __builtin_infl()) - res))
      break;
  }
  if (n >= 100)
  {
    printf ("@@@@ failed\n");
    return -1;
  }

  printf ("#### %d %Lg\n", n, f*z);
  return res / (x * sqrtl (PI));
}

int main(void)
{
  int i;

  for (i = 20; i < 10000; i++)
    printf ("%Lg %Lg\n", i * 0.1L, asympt(i * 0.1L));
  return 0;
}


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33197