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; }