diff options
Diffstat (limited to 'gsl-1.9/integration/util.c')
-rw-r--r-- | gsl-1.9/integration/util.c | 137 |
1 files changed, 137 insertions, 0 deletions
diff --git a/gsl-1.9/integration/util.c b/gsl-1.9/integration/util.c new file mode 100644 index 0000000..8be4f64 --- /dev/null +++ b/gsl-1.9/integration/util.c @@ -0,0 +1,137 @@ +/* integration/util.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. + */ + +static inline +void update (gsl_integration_workspace * workspace, + double a1, double b1, double area1, double error1, + double a2, double b2, double area2, double error2); + +static inline void +retrieve (const gsl_integration_workspace * workspace, + double * a, double * b, double * r, double * e); + + + +static inline +void update (gsl_integration_workspace * workspace, + double a1, double b1, double area1, double error1, + double a2, double b2, double area2, double error2) +{ + double * alist = workspace->alist ; + double * blist = workspace->blist ; + double * rlist = workspace->rlist ; + double * elist = workspace->elist ; + size_t * level = workspace->level ; + + const size_t i_max = workspace->i ; + const size_t i_new = workspace->size ; + + const size_t new_level = workspace->level[i_max] + 1; + + /* append the newly-created intervals to the list */ + + if (error2 > error1) + { + alist[i_max] = a2; /* blist[maxerr] is already == b2 */ + rlist[i_max] = area2; + elist[i_max] = error2; + level[i_max] = new_level; + + alist[i_new] = a1; + blist[i_new] = b1; + rlist[i_new] = area1; + elist[i_new] = error1; + level[i_new] = new_level; + } + else + { + blist[i_max] = b1; /* alist[maxerr] is already == a1 */ + rlist[i_max] = area1; + elist[i_max] = error1; + level[i_max] = new_level; + + alist[i_new] = a2; + blist[i_new] = b2; + rlist[i_new] = area2; + elist[i_new] = error2; + level[i_new] = new_level; + } + + workspace->size++; + + if (new_level > workspace->maximum_level) + { + workspace->maximum_level = new_level; + } + + qpsrt (workspace) ; +} + +static inline void +retrieve (const gsl_integration_workspace * workspace, + double * a, double * b, double * r, double * e) +{ + const size_t i = workspace->i; + double * alist = workspace->alist; + double * blist = workspace->blist; + double * rlist = workspace->rlist; + double * elist = workspace->elist; + + *a = alist[i] ; + *b = blist[i] ; + *r = rlist[i] ; + *e = elist[i] ; +} + +static inline double +sum_results (const gsl_integration_workspace * workspace); + +static inline double +sum_results (const gsl_integration_workspace * workspace) +{ + const double * const rlist = workspace->rlist ; + const size_t n = workspace->size; + + size_t k; + double result_sum = 0; + + for (k = 0; k < n; k++) + { + result_sum += rlist[k]; + } + + return result_sum; +} + +static inline int +subinterval_too_small (double a1, double a2, double b2); + +static inline int +subinterval_too_small (double a1, double a2, double b2) +{ + const double e = GSL_DBL_EPSILON; + const double u = GSL_DBL_MIN; + + double tmp = (1 + 100 * e) * (fabs (a2) + 1000 * u); + + int status = fabs (a1) <= tmp && fabs (b2) <= tmp; + + return status; +} + |