1 /* sum/gsl_sum.h 2 * 3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough 4 * 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. 9 * 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. 14 * 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. 18 */ 19 20 /* Author: G. Jungman */ 21 22 23 #ifndef __GSL_SUM_H__ 24 #define __GSL_SUM_H__ 25 26 #include <stdlib.h> 27 28 #undef __BEGIN_DECLS 29 #undef __END_DECLS 30 #ifdef __cplusplus 31 # define __BEGIN_DECLS extern "C" { 32 # define __END_DECLS } 33 #else 34 # define __BEGIN_DECLS /* empty */ 35 # define __END_DECLS /* empty */ 36 #endif 37 38 __BEGIN_DECLS 39 40 /* Workspace for Levin U Transform with error estimation, 41 * 42 * size = number of terms the workspace can handle 43 * sum_plain = simple sum of series 44 * q_num = backward diagonal of numerator; length = size 45 * q_den = backward diagonal of denominator; length = size 46 * dq_num = table of numerator derivatives; length = size**2 47 * dq_den = table of denominator derivatives; length = size**2 48 * dsum = derivative of sum wrt term i; length = size 49 */ 50 51 typedef struct 52 { 53 size_t size; 54 size_t i; /* position in array */ 55 size_t terms_used; /* number of calls */ 56 double sum_plain; 57 double *q_num; 58 double *q_den; 59 double *dq_num; 60 double *dq_den; 61 double *dsum; 62 } 63 gsl_sum_levin_u_workspace; 64 65 gsl_sum_levin_u_workspace *gsl_sum_levin_u_alloc (size_t n); 66 void gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * w); 67 68 /* Basic Levin-u acceleration method. 69 * 70 * array = array of series elements 71 * n = size of array 72 * sum_accel = result of summation acceleration 73 * err = estimated error 74 * 75 * See [Fessler et al., ACM TOMS 9, 346 (1983) and TOMS-602] 76 */ 77 78 int gsl_sum_levin_u_accel (const double *array, 79 const size_t n, 80 gsl_sum_levin_u_workspace * w, 81 double *sum_accel, double *abserr); 82 83 /* Basic Levin-u acceleration method with constraints on the terms 84 * used, 85 * 86 * array = array of series elements 87 * n = size of array 88 * min_terms = minimum number of terms to sum 89 * max_terms = maximum number of terms to sum 90 * sum_accel = result of summation acceleration 91 * err = estimated error 92 * 93 * See [Fessler et al., ACM TOMS 9, 346 (1983) and TOMS-602] 94 */ 95 96 int gsl_sum_levin_u_minmax (const double *array, 97 const size_t n, 98 const size_t min_terms, 99 const size_t max_terms, 100 gsl_sum_levin_u_workspace * w, 101 double *sum_accel, double *abserr); 102 103 /* Basic Levin-u step w/o reference to the array of terms. 104 * We only need to specify the value of the current term 105 * to execute the step. See TOMS-745. 106 * 107 * sum = t0 + ... + t_{n-1} + term; term = t_{n} 108 * 109 * term = value of the series term to be added 110 * n = position of term in series (starting from 0) 111 * sum_accel = result of summation acceleration 112 * sum_plain = simple sum of series 113 */ 114 115 int 116 gsl_sum_levin_u_step (const double term, 117 const size_t n, 118 const size_t nmax, 119 gsl_sum_levin_u_workspace * w, 120 double *sum_accel); 121 122 /* The following functions perform the same calculation without 123 estimating the errors. They require O(N) storage instead of O(N^2). 124 This may be useful for summing many similar series where the size 125 of the error has already been estimated reliably and is not 126 expected to change. */ 127 128 typedef struct 129 { 130 size_t size; 131 size_t i; /* position in array */ 132 size_t terms_used; /* number of calls */ 133 double sum_plain; 134 double *q_num; 135 double *q_den; 136 double *dsum; 137 } 138 gsl_sum_levin_utrunc_workspace; 139 140 gsl_sum_levin_utrunc_workspace *gsl_sum_levin_utrunc_alloc (size_t n); 141 void gsl_sum_levin_utrunc_free (gsl_sum_levin_utrunc_workspace * w); 142 143 int gsl_sum_levin_utrunc_accel (const double *array, 144 const size_t n, 145 gsl_sum_levin_utrunc_workspace * w, 146 double *sum_accel, double *abserr_trunc); 147 148 int gsl_sum_levin_utrunc_minmax (const double *array, 149 const size_t n, 150 const size_t min_terms, 151 const size_t max_terms, 152 gsl_sum_levin_utrunc_workspace * w, 153 double *sum_accel, double *abserr_trunc); 154 155 int gsl_sum_levin_utrunc_step (const double term, 156 const size_t n, 157 gsl_sum_levin_utrunc_workspace * w, 158 double *sum_accel); 159 160 __END_DECLS 161 162 #endif /* __GSL_SUM_H__ */ 163