https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63661

--- Comment #15 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
Slightly reduced testcase:
static void
foo (double a, double q, double xx, double *ff, double *gx, int e, int ni)
{
  int ii, nn;
  double g1 = *gx;
  ff[ni] = xx;
  if (e == 0)
    {
      for (ii = 0; ii < ni; ii++)
    {
      nn = 100 - ii - 1;
      ff[ni - ii - 1]
        = -1.0 / ((4 * (nn + 1) * (nn + 1) - a) / q + ff[ni - ii]);
    }
    }
  else
    {
      for (ii = 0; ii < ni; ii++)
    {
      nn = 100 - ii - 1;
      ff[ni - ii - 1]
        = -1.0 / (((2 * nn + 1) * (2 * nn + 1) - a) / q + ff[ni - ii]);
    }
    }
  *gx = ff[0] - g1;
}

int
bar (int order, double q, double a, double c[])
{
  int ni, nn, i, e;
  double eps, g1, g2, x1, x2, e1, e2, de, xh, s, ratio, ff[100];
  eps = 1e-10;
  c[0] = 1.0;
  e = 0;
  if (order % 2 != 0)
    e = 1;
  if (q == 0.0)
    {
      for (i = 0; i < 100; i++)
    c[i] = 0.0;
      c[(order - 1) / 2] = 1.0;
      return 0;
    }
  if (order < 5)
    {
      nn = 0;
      s = 0.0;
      if (e == 0)
    ratio = (a - 4) / q;
      else
    ratio = (a - 1 - q) / q;
    }
  else
    {
      if (e == 0)
    {
      c[1] = (a - 4) / q;
      s = 2 * c[0] + 4 * c[1];
      for (i = 2; i < order / 2; i++)
        {
          c[i] = (a - 4 * i * i) / q * c[i - 1] - c[i - 2];
          s += 2 * (i + 1) * c[i];
        }
    }
      else
    {
      c[1] = (a - 1) / q + 1;
      s = c[0] + 3 * c[1];
      for (i = 2; i < order / 2 + 1; i++)
        {
          c[i] = (a - (2 * i - 1) * (2 * i - 1)) / q * c[i - 1] - c[i - 2];
          s += (2 * (i + 1) - 1) * c[i];
        }
    }
      nn = i - 1;
      ratio = c[nn] / c[nn - 1];
    }
  ni = 100 - nn - 1;
  if (e == 0)
    x1 = -q / (4.0 * (100 + 1.0) * (100 + 1.0));
  else
    x1 = -q / ((2.0 * 100 + 1.0) * (2.0 * 100 + 1.0));
  g1 = ratio;
  foo (a, q, x1, ff, &g1, e, ni);
  x2 = g1;
  g2 = ratio;
  while (1)
    {
      __builtin_printf ("loop %g\n", g2);
      e1 = g1 - x1;
      e2 = g2 - x2;
      de = e1 - e2;
      if (__builtin_fabs (de) < eps)
    break;
      xh = (e1 * x2 - e2 * x1) / de;
      x1 = x2;
      g1 = g2;
      x2 = xh;
      g2 = ratio;
      foo (a, q, x2, ff, &g2, e, ni);
    }
  s += 2 * (nn + 1) * c[nn];
  for (i = nn + 1; i < 100; i++)
    {
      c[i] = ff[i - nn - 1] * c[i - 1];
      s += 2 * (i + 1) * c[i];
      if (__builtin_fabs (c[i]) < 1e-20)
    for (; i < 100;)
      c[i++] = 0.0;
    }
  for (i = 0; i < 100; i++)
    c[i] /= s;
  return 0;
}

int
main ()
{
  double c[100];
  bar (1, 5.0, -5.790080598637771, c);
  return 0;
}

Reply via email to