| /* integration/qelg.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or (at |
| * your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, but |
| * WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| * General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, write to the Free Software |
| * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
| */ |
| |
| struct extrapolation_table |
| { |
| size_t n; |
| double rlist2[52]; |
| size_t nres; |
| double res3la[3]; |
| }; |
| |
| static void |
| initialise_table (struct extrapolation_table *table); |
| |
| static void |
| append_table (struct extrapolation_table *table, double y); |
| |
| static void |
| initialise_table (struct extrapolation_table *table) |
| { |
| table->n = 0; |
| table->nres = 0; |
| } |
| #ifdef JUNK |
| static void |
| initialise_table (struct extrapolation_table *table, double y) |
| { |
| table->n = 0; |
| table->rlist2[0] = y; |
| table->nres = 0; |
| } |
| #endif |
| static void |
| append_table (struct extrapolation_table *table, double y) |
| { |
| size_t n; |
| n = table->n; |
| table->rlist2[n] = y; |
| table->n++; |
| } |
| |
| /* static inline void |
| qelg (size_t * n, double epstab[], |
| double * result, double * abserr, |
| double res3la[], size_t * nres); */ |
| |
| static inline void |
| qelg (struct extrapolation_table *table, double *result, double *abserr); |
| |
| static inline void |
| qelg (struct extrapolation_table *table, double *result, double *abserr) |
| { |
| double *epstab = table->rlist2; |
| double *res3la = table->res3la; |
| const size_t n = table->n - 1; |
| |
| const double current = epstab[n]; |
| |
| double absolute = GSL_DBL_MAX; |
| double relative = 5 * GSL_DBL_EPSILON * fabs (current); |
| |
| const size_t newelm = n / 2; |
| const size_t n_orig = n; |
| size_t n_final = n; |
| size_t i; |
| |
| const size_t nres_orig = table->nres; |
| |
| *result = current; |
| *abserr = GSL_DBL_MAX; |
| |
| if (n < 2) |
| { |
| *result = current; |
| *abserr = GSL_MAX_DBL (absolute, relative); |
| return; |
| } |
| |
| epstab[n + 2] = epstab[n]; |
| epstab[n] = GSL_DBL_MAX; |
| |
| for (i = 0; i < newelm; i++) |
| { |
| double res = epstab[n - 2 * i + 2]; |
| double e0 = epstab[n - 2 * i - 2]; |
| double e1 = epstab[n - 2 * i - 1]; |
| double e2 = res; |
| |
| double e1abs = fabs (e1); |
| double delta2 = e2 - e1; |
| double err2 = fabs (delta2); |
| double tol2 = GSL_MAX_DBL (fabs (e2), e1abs) * GSL_DBL_EPSILON; |
| double delta3 = e1 - e0; |
| double err3 = fabs (delta3); |
| double tol3 = GSL_MAX_DBL (e1abs, fabs (e0)) * GSL_DBL_EPSILON; |
| |
| double e3, delta1, err1, tol1, ss; |
| |
| if (err2 <= tol2 && err3 <= tol3) |
| { |
| /* If e0, e1 and e2 are equal to within machine accuracy, |
| convergence is assumed. */ |
| |
| *result = res; |
| absolute = err2 + err3; |
| relative = 5 * GSL_DBL_EPSILON * fabs (res); |
| *abserr = GSL_MAX_DBL (absolute, relative); |
| return; |
| } |
| |
| e3 = epstab[n - 2 * i]; |
| epstab[n - 2 * i] = e1; |
| delta1 = e1 - e3; |
| err1 = fabs (delta1); |
| tol1 = GSL_MAX_DBL (e1abs, fabs (e3)) * GSL_DBL_EPSILON; |
| |
| /* If two elements are very close to each other, omit a part of |
| the table by adjusting the value of n */ |
| |
| if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3) |
| { |
| n_final = 2 * i; |
| break; |
| } |
| |
| ss = (1 / delta1 + 1 / delta2) - 1 / delta3; |
| |
| /* Test to detect irregular behaviour in the table, and |
| eventually omit a part of the table by adjusting the value of |
| n. */ |
| |
| if (fabs (ss * e1) <= 0.0001) |
| { |
| n_final = 2 * i; |
| break; |
| } |
| |
| /* Compute a new element and eventually adjust the value of |
| result. */ |
| |
| res = e1 + 1 / ss; |
| epstab[n - 2 * i] = res; |
| |
| { |
| const double error = err2 + fabs (res - e2) + err3; |
| |
| if (error <= *abserr) |
| { |
| *abserr = error; |
| *result = res; |
| } |
| } |
| } |
| |
| /* Shift the table */ |
| |
| { |
| const size_t limexp = 50 - 1; |
| |
| if (n_final == limexp) |
| { |
| n_final = 2 * (limexp / 2); |
| } |
| } |
| |
| if (n_orig % 2 == 1) |
| { |
| for (i = 0; i <= newelm; i++) |
| { |
| epstab[1 + i * 2] = epstab[i * 2 + 3]; |
| } |
| } |
| else |
| { |
| for (i = 0; i <= newelm; i++) |
| { |
| epstab[i * 2] = epstab[i * 2 + 2]; |
| } |
| } |
| |
| if (n_orig != n_final) |
| { |
| for (i = 0; i <= n_final; i++) |
| { |
| epstab[i] = epstab[n_orig - n_final + i]; |
| } |
| } |
| |
| table->n = n_final + 1; |
| |
| if (nres_orig < 3) |
| { |
| res3la[nres_orig] = *result; |
| *abserr = GSL_DBL_MAX; |
| } |
| else |
| { /* Compute error estimate */ |
| *abserr = (fabs (*result - res3la[2]) + fabs (*result - res3la[1]) |
| + fabs (*result - res3la[0])); |
| |
| res3la[0] = res3la[1]; |
| res3la[1] = res3la[2]; |
| res3la[2] = *result; |
| } |
| |
| /* In QUADPACK the variable table->nres is incremented at the top of |
| qelg, so it increases on every call. This leads to the array |
| res3la being accessed when its elements are still undefined, so I |
| have moved the update to this point so that its value more |
| useful. */ |
| |
| table->nres = nres_orig + 1; |
| |
| *abserr = GSL_MAX_DBL (*abserr, 5 * GSL_DBL_EPSILON * fabs (*result)); |
| |
| return; |
| } |