| /* sum/levin_utrunc.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, 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. |
| */ |
| |
| /* Author: G. Jungman */ |
| |
| #include <config.h> |
| #include <gsl/gsl_math.h> |
| #include <gsl/gsl_test.h> |
| #include <gsl/gsl_errno.h> |
| #include <gsl/gsl_sum.h> |
| |
| int |
| gsl_sum_levin_utrunc_accel (const double *array, |
| const size_t array_size, |
| gsl_sum_levin_utrunc_workspace * w, |
| double *sum_accel, double *abserr_trunc) |
| { |
| return gsl_sum_levin_utrunc_minmax (array, array_size, |
| 0, array_size - 1, |
| w, sum_accel, abserr_trunc); |
| } |
| |
| |
| int |
| gsl_sum_levin_utrunc_minmax (const double *array, |
| const size_t array_size, |
| const size_t min_terms, |
| const size_t max_terms, |
| gsl_sum_levin_utrunc_workspace * w, |
| double *sum_accel, double *abserr_trunc) |
| { |
| if (array_size == 0) |
| { |
| *sum_accel = 0.0; |
| *abserr_trunc = 0.0; |
| w->sum_plain = 0.0; |
| w->terms_used = 0; |
| return GSL_SUCCESS; |
| } |
| else if (array_size == 1) |
| { |
| *sum_accel = array[0]; |
| *abserr_trunc = GSL_POSINF; |
| w->sum_plain = array[0]; |
| w->terms_used = 1; |
| return GSL_SUCCESS; |
| } |
| else |
| { |
| const double SMALL = 0.01; |
| const size_t nmax = GSL_MAX (max_terms, array_size) - 1; |
| double trunc_n = 0.0, trunc_nm1 = 0.0; |
| double actual_trunc_n = 0.0, actual_trunc_nm1 = 0.0; |
| double result_n = 0.0, result_nm1 = 0.0; |
| size_t n; |
| int better = 0; |
| int before = 0; |
| int converging = 0; |
| double least_trunc = GSL_DBL_MAX; |
| double result_least_trunc; |
| |
| /* Calculate specified minimum number of terms. No convergence |
| tests are made, and no truncation information is stored. */ |
| |
| for (n = 0; n < min_terms; n++) |
| { |
| const double t = array[n]; |
| |
| result_nm1 = result_n; |
| gsl_sum_levin_utrunc_step (t, n, w, &result_n); |
| } |
| |
| /* Assume the result after the minimum calculation is the best. */ |
| |
| result_least_trunc = result_n; |
| |
| /* Calculate up to maximum number of terms. Check truncation |
| condition. */ |
| |
| for (; n <= nmax; n++) |
| { |
| const double t = array[n]; |
| |
| result_nm1 = result_n; |
| gsl_sum_levin_utrunc_step (t, n, w, &result_n); |
| |
| /* Compute the truncation error directly */ |
| |
| actual_trunc_nm1 = actual_trunc_n; |
| actual_trunc_n = fabs (result_n - result_nm1); |
| |
| /* Average results to make a more reliable estimate of the |
| real truncation error */ |
| |
| trunc_nm1 = trunc_n; |
| trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1); |
| |
| /* Determine if we are in the convergence region. */ |
| |
| better = (trunc_n < trunc_nm1 || trunc_n < SMALL * fabs (result_n)); |
| converging = converging || (better && before); |
| before = better; |
| |
| if (converging) |
| { |
| if (trunc_n < least_trunc) |
| { |
| /* Found a low truncation point in the convergence |
| region. Save it. */ |
| |
| least_trunc = trunc_n; |
| result_least_trunc = result_n; |
| } |
| |
| if (fabs (trunc_n / result_n) < 10.0 * GSL_MACH_EPS) |
| break; |
| } |
| } |
| |
| if (converging) |
| { |
| /* Stopped in the convergence region. Return result and |
| error estimate. */ |
| |
| *sum_accel = result_least_trunc; |
| *abserr_trunc = least_trunc; |
| w->terms_used = n; |
| return GSL_SUCCESS; |
| } |
| else |
| { |
| /* Never reached the convergence region. Use the last |
| calculated values. */ |
| |
| *sum_accel = result_n; |
| *abserr_trunc = trunc_n; |
| w->terms_used = n; |
| return GSL_SUCCESS; |
| } |
| } |
| } |
| |
| int |
| gsl_sum_levin_utrunc_step (const double term, |
| const size_t n, |
| gsl_sum_levin_utrunc_workspace * w, double *sum_accel) |
| { |
| if (term == 0.0) |
| { |
| /* This is actually harmless when treated in this way. A term |
| which is exactly zero is simply ignored; the state is not |
| changed. We return GSL_EZERODIV as an indicator that this |
| occured. */ |
| |
| return GSL_EZERODIV; |
| } |
| else if (n == 0) |
| { |
| *sum_accel = term; |
| w->sum_plain = term; |
| w->q_den[0] = 1.0 / term; |
| w->q_num[0] = 1.0; |
| return GSL_SUCCESS; |
| } |
| else |
| { |
| double factor = 1.0; |
| double ratio = (double) n / (n + 1.0); |
| int j; |
| |
| w->sum_plain += term; |
| w->q_den[n] = 1.0 / (term * (n + 1.0) * (n + 1.0)); |
| w->q_num[n] = w->sum_plain * w->q_den[n]; |
| |
| for (j = n - 1; j >= 0; j--) |
| { |
| double c = factor * (j + 1) / (n + 1); |
| factor *= ratio; |
| w->q_den[j] = w->q_den[j + 1] - c * w->q_den[j]; |
| w->q_num[j] = w->q_num[j + 1] - c * w->q_num[j]; |
| } |
| |
| *sum_accel = w->q_num[0] / w->q_den[0]; |
| return GSL_SUCCESS; |
| } |
| } |