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.
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_integration.h>
27 #include "initialise.c"
28 #include "set_initial.c"
34 #include "positivity.c"
39 gsl_integration_qawo (gsl_function * f,
41 const double epsabs, const double epsrel,
43 gsl_integration_workspace * workspace,
44 gsl_integration_qawo_table * wf,
45 double *result, double *abserr)
48 double res_ext, err_ext;
49 double result0, abserr0, resabs0, resasc0;
53 double error_over_large_intervals = 0;
54 double reseps = 0, abseps = 0, correc = 0;
56 int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
57 int error_type = 0, error_type2 = 0;
61 int positive_integrand = 0;
64 int disallow_extrapolation = 0;
66 struct extrapolation_table table;
68 double b = a + wf->L ;
69 double abs_omega = fabs (wf->omega) ;
71 /* Initialize results */
73 initialise (workspace, a, b);
78 if (limit > workspace->limit)
80 GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
83 /* Test on accuracy */
85 if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
87 GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
91 /* Perform the first integration */
93 qc25f (f, a, b, wf, 0, &result0, &abserr0, &resabs0, &resasc0);
95 set_initial_result (workspace, result0, abserr0);
97 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
99 if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance)
104 GSL_ERROR ("cannot reach tolerance because of roundoff error"
105 "on first attempt", GSL_EROUND);
107 else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
119 GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
124 initialise_table (&table);
126 if (0.5 * abs_omega * fabs(b - a) <= 2)
128 append_table (&table, result0);
136 err_ext = GSL_DBL_MAX;
138 positive_integrand = test_positivity (result0, resabs0);
144 size_t current_level;
145 double a1, b1, a2, b2;
146 double a_i, b_i, r_i, e_i;
147 double area1 = 0, area2 = 0, area12 = 0;
148 double error1 = 0, error2 = 0, error12 = 0;
149 double resasc1, resasc2;
150 double resabs1, resabs2;
153 /* Bisect the subinterval with the largest error estimate */
155 retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
157 current_level = workspace->level[workspace->i] + 1;
159 if (current_level >= wf->n)
161 error_type = -1 ; /* exceeded limit of table */
166 b1 = 0.5 * (a_i + b_i);
172 qc25f (f, a1, b1, wf, current_level, &area1, &error1, &resabs1, &resasc1);
173 qc25f (f, a2, b2, wf, current_level, &area2, &error2, &resabs2, &resasc2);
175 area12 = area1 + area2;
176 error12 = error1 + error2;
179 /* Improve previous approximations to the integral and test for
182 We write these expressions in the same way as the original
183 QUADPACK code so that the rounding errors are the same, which
184 makes testing easier. */
186 errsum = errsum + error12 - e_i;
187 area = area + area12 - r_i;
189 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
191 if (resasc1 != error1 && resasc2 != error2)
193 double delta = r_i - area12;
195 if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
206 if (iteration > 10 && error12 > e_i)
212 /* Test for roundoff and eventually set error flag */
214 if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
216 error_type = 2; /* round off error */
219 if (roundoff_type2 >= 5)
224 /* set error flag in the case of bad integrand behaviour at
225 a point of the integration range */
227 if (subinterval_too_small (a1, a2, b2))
232 /* append the newly-created intervals to the list */
234 update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
236 if (errsum <= tolerance)
246 if (iteration >= limit - 1)
252 /* set up variables on first iteration */
254 if (iteration == 2 && extall)
256 error_over_large_intervals = errsum;
258 append_table (&table, area);
262 if (disallow_extrapolation)
269 error_over_large_intervals += -last_e_i;
271 if (current_level < workspace->maximum_level)
273 error_over_large_intervals += error12;
280 if (large_interval(workspace))
288 workspace->nrmax = 1;
292 /* test whether the interval to be bisected next is the
293 smallest interval. */
294 size_t i = workspace->i;
295 double width = workspace->blist[i] - workspace->alist[i];
297 if (0.25 * fabs(width) * abs_omega > 2)
301 error_over_large_intervals = errsum;
307 if (!error_type2 && error_over_large_intervals > ertest)
309 if (increase_nrmax (workspace))
313 /* Perform extrapolation */
315 append_table (&table, area);
319 reset_nrmax(workspace);
321 error_over_large_intervals = errsum;
325 qelg (&table, &reseps, &abseps);
329 if (ktmin > 5 && err_ext < 0.001 * errsum)
334 if (abseps < err_ext)
339 correc = error_over_large_intervals;
340 ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
341 if (err_ext <= ertest)
345 /* Prepare bisection of the smallest interval. */
349 disallow_extrapolation = 1;
357 /* work on interval with largest error */
359 reset_nrmax (workspace);
361 error_over_large_intervals = errsum;
364 while (iteration < limit);
369 if (err_ext == GSL_DBL_MAX)
372 if (error_type || error_type2)
382 if (result != 0 && area != 0)
384 if (err_ext / fabs (res_ext) > errsum / fabs (area))
387 else if (err_ext > errsum)
391 else if (area == 0.0)
397 /* Test on divergence. */
400 double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
402 if (!positive_integrand && max_area < 0.01 * resabs0)
407 double ratio = res_ext / area;
409 if (ratio < 0.01 || ratio > 100 || errsum > fabs (area))
417 *result = sum_results (workspace);
429 else if (error_type == 1)
431 GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
433 else if (error_type == 2)
435 GSL_ERROR ("cannot reach tolerance because of roundoff error",
438 else if (error_type == 3)
440 GSL_ERROR ("bad integrand behavior found in the integration interval",
443 else if (error_type == 4)
445 GSL_ERROR ("roundoff error detected in extrapolation table", GSL_EROUND);
447 else if (error_type == 5)
449 GSL_ERROR ("integral is divergent, or slowly convergent", GSL_EDIVERGE);
451 else if (error_type == -1)
453 GSL_ERROR ("exceeded limit of trigonometric table", GSL_ETABLE);
457 GSL_ERROR ("could not integrate function", GSL_EFAILED);