42 #ifndef O2SCL_SERIES_ACC_H 43 #define O2SCL_SERIES_ACC_H 51 #include <gsl/gsl_sum.h> 52 #include <gsl/gsl_minmax.h> 53 #include <gsl/gsl_nan.h> 54 #include <gsl/gsl_machine.h> 55 #include <o2scl/err_hnd.h> 57 #ifndef DOXYGEN_NO_O2NS 86 double series_accel(
size_t na, vec_t &array,
double &abserr_trunc) {
92 size_t max_terms=
size-1;
102 }
else if (
size == 1) {
105 abserr_trunc=GSL_POSINF;
106 wt->sum_plain=array[0];
112 const double SMALL=0.01;
113 const size_t nmax=GSL_MAX(max_terms,
size)-1;
114 double trunc_n=0.0, trunc_nm1=0.0;
115 double actual_trunc_n=0.0, actual_trunc_nm1=0.0;
116 double result_n=0.0, result_nm1=0.0;
121 double least_trunc=GSL_DBL_MAX;
122 double result_least_trunc;
127 for (n=0; n < min_terms; n++) {
128 const double t=array[n];
136 result_least_trunc=result_n;
141 for (; n <= nmax; n++) {
142 const double t=array[n];
149 actual_trunc_nm1=actual_trunc_n;
150 actual_trunc_n=fabs(result_n-result_nm1);
156 trunc_n=0.5*(actual_trunc_n+actual_trunc_nm1);
160 better=(trunc_n < trunc_nm1 || trunc_n < SMALL*fabs(result_n));
161 converging=converging || (better && before);
165 if (trunc_n < least_trunc) {
170 result_least_trunc=result_n;
173 if (fabs(trunc_n/result_n) < 10.0*GSL_MACH_EPS)
182 sum_accel=result_least_trunc;
183 abserr_trunc=least_trunc;
192 abserr_trunc=trunc_n;
208 template<
class vec_t>
215 size_t max_terms=
size-1;
220 while (size2 > 0 && array[size2-1]==0) {
232 }
else if (size2==1) {
236 w->sum_plain=array[0];
242 const double SMALL=0.01;
243 const size_t nmax=GSL_MAX(max_terms,
size)-1;
244 double noise_n=0.0,noise_nm1=0.0;
245 double trunc_n=0.0, trunc_nm1=0.0;
246 double actual_trunc_n=0.0, actual_trunc_nm1=0.0;
247 double result_n=0.0, result_nm1=0.0;
254 double least_trunc=GSL_DBL_MAX;
255 double least_trunc_noise=GSL_DBL_MAX;
256 double least_trunc_result;
261 for (n=0; n < min_terms; n++) {
262 const double t=array[n];
267 least_trunc_result=result_n;
270 for (i=0; i < n; i++) {
271 double dn=
w->dsum[i]*GSL_MACH_EPS*array[i];
274 noise_n=std::sqrt(variance);
279 for (; n <= nmax; n++) {
280 const double t=array[n];
287 actual_trunc_nm1=actual_trunc_n;
288 actual_trunc_n=fabs(result_n-result_nm1);
294 trunc_n=0.5*(actual_trunc_n+actual_trunc_nm1);
299 for (i=0; i <= n; i++) {
300 double dn=
w->dsum[i]*GSL_MACH_EPS*array[i];
304 noise_n=std::sqrt(variance);
308 better=(trunc_n < trunc_nm1 ||
309 trunc_n < SMALL*fabs(result_n));
310 converging=converging || (better && before);
314 if (trunc_n < least_trunc) {
318 least_trunc_result=result_n;
320 least_trunc_noise=noise_n;
323 if (noise_n > trunc_n/3.0) {
327 if (trunc_n < 10.0*GSL_MACH_EPS*fabs(result_n)) {
338 sum_accel=least_trunc_result;
339 abserr=GSL_MAX_DBL(least_trunc,least_trunc_noise);
348 abserr=GSL_MAX_DBL(trunc_n,noise_n);
360 #ifndef DOXYGEN_NO_O2NS 374 const size_t nmax,
double &sum_accel);
382 gsl_sum_levin_u_workspace *
w;
385 gsl_sum_levin_utrunc_workspace *
wt;
394 #ifndef DOXYGEN_NO_O2NS The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
double series_accel_err(size_t na, vec_t &array, double &abserr)
Return the accelerated sum of the series with an accurate error estimate.
gsl_sum_levin_utrunc_workspace * wt
The GSL workspace.
double series_accel(size_t na, vec_t &array, double &abserr_trunc)
Return the accelerated sum of the series with a simple error estimate.
int levin_u_step(const double term, const size_t n, const size_t nmax, double &sum_accel)
Perform a step.
size_t size
The workspace size.
Series acceleration by Levin u-transform (GSL)
void set_size(size_t new_size)
Set the number of terms.
series_acc(size_t size=0)
size is the number of terms in the series
int levin_utrunc_step(const double term, const size_t n, double &sum_accel)
Perform a step.
size_t series_index(size_t i, size_t j, size_t nmax)
An internal function reducing two matrix indices, i and j, to index of a single array.
gsl_sum_levin_u_workspace * w
The GSL workspace.