diff options
Diffstat (limited to 'gsl-1.9/integration/qcheb.c')
| -rw-r--r-- | gsl-1.9/integration/qcheb.c | 229 |
1 files changed, 229 insertions, 0 deletions
diff --git a/gsl-1.9/integration/qcheb.c b/gsl-1.9/integration/qcheb.c new file mode 100644 index 0000000..a6564f7 --- /dev/null +++ b/gsl-1.9/integration/qcheb.c @@ -0,0 +1,229 @@ +/* integration/qcheb.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. + */ + +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_integration.h> + + +/* This function computes the 12-th order and 24-th order Chebyshev + approximations to f(x) on [a,b] */ + +void +gsl_integration_qcheb (gsl_function * f, double a, double b, double *cheb12, double *cheb24) +{ + size_t i; + double fval[25], v[12]; + + /* These are the values of cos(pi*k/24) for k=1..11 needed for the + Chebyshev expansion of f(x) */ + + const double x[11] = { 0.9914448613738104, + 0.9659258262890683, + 0.9238795325112868, + 0.8660254037844386, + 0.7933533402912352, + 0.7071067811865475, + 0.6087614290087206, + 0.5000000000000000, + 0.3826834323650898, + 0.2588190451025208, + 0.1305261922200516 }; + + const double center = 0.5 * (b + a); + const double half_length = 0.5 * (b - a); + + fval[0] = 0.5 * GSL_FN_EVAL (f, b); + fval[12] = GSL_FN_EVAL (f, center); + fval[24] = 0.5 * GSL_FN_EVAL (f, a); + + for (i = 1; i < 12; i++) + { + const size_t j = 24 - i; + const double u = half_length * x[i-1]; + fval[i] = GSL_FN_EVAL(f, center + u); + fval[j] = GSL_FN_EVAL(f, center - u); + } + + for (i = 0; i < 12; i++) + { + const size_t j = 24 - i; + v[i] = fval[i] - fval[j]; + fval[i] = fval[i] + fval[j]; + } + + { + const double alam1 = v[0] - v[8]; + const double alam2 = x[5] * (v[2] - v[6] - v[10]); + + cheb12[3] = alam1 + alam2; + cheb12[9] = alam1 - alam2; + } + + { + const double alam1 = v[1] - v[7] - v[9]; + const double alam2 = v[3] - v[5] - v[11]; + { + const double alam = x[2] * alam1 + x[8] * alam2; + + cheb24[3] = cheb12[3] + alam; + cheb24[21] = cheb12[3] - alam; + } + + { + const double alam = x[8] * alam1 - x[2] * alam2; + cheb24[9] = cheb12[9] + alam; + cheb24[15] = cheb12[9] - alam; + } + } + + { + const double part1 = x[3] * v[4]; + const double part2 = x[7] * v[8]; + const double part3 = x[5] * v[6]; + + { + const double alam1 = v[0] + part1 + part2; + const double alam2 = x[1] * v[2] + part3 + x[9] * v[10]; + + cheb12[1] = alam1 + alam2; + cheb12[11] = alam1 - alam2; + } + + { + const double alam1 = v[0] - part1 + part2; + const double alam2 = x[9] * v[2] - part3 + x[1] * v[10]; + cheb12[5] = alam1 + alam2; + cheb12[7] = alam1 - alam2; + } + } + + { + const double alam = (x[0] * v[1] + x[2] * v[3] + x[4] * v[5] + + x[6] * v[7] + x[8] * v[9] + x[10] * v[11]); + cheb24[1] = cheb12[1] + alam; + cheb24[23] = cheb12[1] - alam; + } + + { + const double alam = (x[10] * v[1] - x[8] * v[3] + x[6] * v[5] + - x[4] * v[7] + x[2] * v[9] - x[0] * v[11]); + cheb24[11] = cheb12[11] + alam; + cheb24[13] = cheb12[11] - alam; + } + + { + const double alam = (x[4] * v[1] - x[8] * v[3] - x[0] * v[5] + - x[10] * v[7] + x[2] * v[9] + x[6] * v[11]); + cheb24[5] = cheb12[5] + alam; + cheb24[19] = cheb12[5] - alam; + } + + { + const double alam = (x[6] * v[1] - x[2] * v[3] - x[10] * v[5] + + x[0] * v[7] - x[8] * v[9] - x[4] * v[11]); + cheb24[7] = cheb12[7] + alam; + cheb24[17] = cheb12[7] - alam; + } + + for (i = 0; i < 6; i++) + { + const size_t j = 12 - i; + v[i] = fval[i] - fval[j]; + fval[i] = fval[i] + fval[j]; + } + + { + const double alam1 = v[0] + x[7] * v[4]; + const double alam2 = x[3] * v[2]; + + cheb12[2] = alam1 + alam2; + cheb12[10] = alam1 - alam2; + } + + cheb12[6] = v[0] - v[4]; + + { + const double alam = x[1] * v[1] + x[5] * v[3] + x[9] * v[5]; + cheb24[2] = cheb12[2] + alam; + cheb24[22] = cheb12[2] - alam; + } + + { + const double alam = x[5] * (v[1] - v[3] - v[5]); + cheb24[6] = cheb12[6] + alam; + cheb24[18] = cheb12[6] - alam; + } + + { + const double alam = x[9] * v[1] - x[5] * v[3] + x[1] * v[5]; + cheb24[10] = cheb12[10] + alam; + cheb24[14] = cheb12[10] - alam; + } + + for (i = 0; i < 3; i++) + { + const size_t j = 6 - i; + v[i] = fval[i] - fval[j]; + fval[i] = fval[i] + fval[j]; + } + + cheb12[4] = v[0] + x[7] * v[2]; + cheb12[8] = fval[0] - x[7] * fval[2]; + + { + const double alam = x[3] * v[1]; + cheb24[4] = cheb12[4] + alam; + cheb24[20] = cheb12[4] - alam; + } + + { + const double alam = x[7] * fval[1] - fval[3]; + cheb24[8] = cheb12[8] + alam; + cheb24[16] = cheb12[8] - alam; + } + + cheb12[0] = fval[0] + fval[2]; + + { + const double alam = fval[1] + fval[3]; + cheb24[0] = cheb12[0] + alam; + cheb24[24] = cheb12[0] - alam; + } + + cheb12[12] = v[0] - v[2]; + cheb24[12] = cheb12[12]; + + for (i = 1; i < 12; i++) + { + cheb12[i] *= 1.0 / 6.0; + } + + cheb12[0] *= 1.0 / 12.0; + cheb12[12] *= 1.0 / 12.0; + + for (i = 1; i < 24; i++) + { + cheb24[i] *= 1.0 / 12.0; + } + + cheb24[0] *= 1.0 / 24.0; + cheb24[24] *= 1.0 / 24.0; +} |