3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
20 struct extrapolation_table
29 initialise_table (struct extrapolation_table *table);
32 append_table (struct extrapolation_table *table, double y);
35 initialise_table (struct extrapolation_table *table)
42 initialise_table (struct extrapolation_table *table, double y)
50 append_table (struct extrapolation_table *table, double y)
59 qelg (size_t * n, double epstab[],
60 double * result, double * abserr,
61 double res3la[], size_t * nres); */
64 qelg (struct extrapolation_table *table, double *result, double *abserr);
67 qelg (struct extrapolation_table *table, double *result, double *abserr)
69 double *epstab = table->rlist2;
70 double *res3la = table->res3la;
71 const size_t n = table->n - 1;
73 const double current = epstab[n];
75 double absolute = GSL_DBL_MAX;
76 double relative = 5 * GSL_DBL_EPSILON * fabs (current);
78 const size_t newelm = n / 2;
79 const size_t n_orig = n;
83 const size_t nres_orig = table->nres;
86 *abserr = GSL_DBL_MAX;
91 *abserr = GSL_MAX_DBL (absolute, relative);
95 epstab[n + 2] = epstab[n];
96 epstab[n] = GSL_DBL_MAX;
98 for (i = 0; i < newelm; i++)
100 double res = epstab[n - 2 * i + 2];
101 double e0 = epstab[n - 2 * i - 2];
102 double e1 = epstab[n - 2 * i - 1];
105 double e1abs = fabs (e1);
106 double delta2 = e2 - e1;
107 double err2 = fabs (delta2);
108 double tol2 = GSL_MAX_DBL (fabs (e2), e1abs) * GSL_DBL_EPSILON;
109 double delta3 = e1 - e0;
110 double err3 = fabs (delta3);
111 double tol3 = GSL_MAX_DBL (e1abs, fabs (e0)) * GSL_DBL_EPSILON;
113 double e3, delta1, err1, tol1, ss;
115 if (err2 <= tol2 && err3 <= tol3)
117 /* If e0, e1 and e2 are equal to within machine accuracy,
118 convergence is assumed. */
121 absolute = err2 + err3;
122 relative = 5 * GSL_DBL_EPSILON * fabs (res);
123 *abserr = GSL_MAX_DBL (absolute, relative);
127 e3 = epstab[n - 2 * i];
128 epstab[n - 2 * i] = e1;
130 err1 = fabs (delta1);
131 tol1 = GSL_MAX_DBL (e1abs, fabs (e3)) * GSL_DBL_EPSILON;
133 /* If two elements are very close to each other, omit a part of
134 the table by adjusting the value of n */
136 if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3)
142 ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
144 /* Test to detect irregular behaviour in the table, and
145 eventually omit a part of the table by adjusting the value of
148 if (fabs (ss * e1) <= 0.0001)
154 /* Compute a new element and eventually adjust the value of
158 epstab[n - 2 * i] = res;
161 const double error = err2 + fabs (res - e2) + err3;
163 if (error <= *abserr)
171 /* Shift the table */
174 const size_t limexp = 50 - 1;
176 if (n_final == limexp)
178 n_final = 2 * (limexp / 2);
184 for (i = 0; i <= newelm; i++)
186 epstab[1 + i * 2] = epstab[i * 2 + 3];
191 for (i = 0; i <= newelm; i++)
193 epstab[i * 2] = epstab[i * 2 + 2];
197 if (n_orig != n_final)
199 for (i = 0; i <= n_final; i++)
201 epstab[i] = epstab[n_orig - n_final + i];
205 table->n = n_final + 1;
209 res3la[nres_orig] = *result;
210 *abserr = GSL_DBL_MAX;
213 { /* Compute error estimate */
214 *abserr = (fabs (*result - res3la[2]) + fabs (*result - res3la[1])
215 + fabs (*result - res3la[0]));
217 res3la[0] = res3la[1];
218 res3la[1] = res3la[2];
222 /* In QUADPACK the variable table->nres is incremented at the top of
223 qelg, so it increases on every call. This leads to the array
224 res3la being accessed when its elements are still undefined, so I
225 have moved the update to this point so that its value more
228 table->nres = nres_orig + 1;
230 *abserr = GSL_MAX_DBL (*abserr, 5 * GSL_DBL_EPSILON * fabs (*result));