series_acc.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2009 Gerard
24  * Jungman, Brian Gough
25  *
26  * This program is free software; you can redistribute it and/or modify
27  * it under the terms of the GNU General Public License as published by
28  * the Free Software Foundation; either version 3 of the License, or (at
29  * your option) any later version.
30  *
31  * This program is distributed in the hope that it will be useful, but
32  * WITHOUT ANY WARRANTY; without even the implied warranty of
33  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
34  * General Public License for more details.
35  *
36  * You should have received a copy of the GNU General Public License
37  * along with this program; if not, write to the Free Software
38  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
39  * 02110-1301, USA.
40  *
41  */
42 #ifndef O2SCL_SERIES_ACC_H
43 #define O2SCL_SERIES_ACC_H
44 
45 /** \file series_acc.h
46  \brief File defining \ref o2scl::series_acc
47 */
48 
49 #include <iostream>
50 #include <cmath>
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>
56 
57 #ifndef DOXYGEN_NO_O2NS
58 namespace o2scl {
59 #endif
60 
61  /** \brief Series acceleration by Levin u-transform (GSL)
62 
63  Given an array of terms in a sum, this attempts to evaluate the
64  entire sum with an estimate of the error.
65 
66  \future Move the workspaces to classes?
67  \future Create an example
68  */
69  class series_acc {
70 
71  public:
72 
73  /** \brief \c size is the number of terms in the series
74  */
75  series_acc(size_t size=0);
76 
77  virtual ~series_acc();
78 
79  /** \brief Return the accelerated sum of the series with
80  a simple error estimate
81 
82  The input vector \c x should be an array with \c n values
83  from <tt>x[0]</tt> to <tt>x[n-1]</tt>.
84  */
85  template<class vec_t>
86  double series_accel(size_t na, vec_t &array, double &abserr_trunc) {
87 
88  if (na!=size) set_size(na);
89 
90  double sum_accel;
91  size_t min_terms=0;
92  size_t max_terms=size-1;
93 
94  if (size == 0) {
95 
96  sum_accel=0.0;
97  abserr_trunc=0.0;
98  wt->sum_plain=0.0;
99  wt->terms_used=0;
100  return sum_accel;
101 
102  } else if (size == 1) {
103 
104  sum_accel=array[0];
105  abserr_trunc=GSL_POSINF;
106  wt->sum_plain=array[0];
107  wt->terms_used=1;
108  return sum_accel;
109 
110  } else {
111 
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;
117  size_t n;
118  int better=0;
119  int before=0;
120  int converging=0;
121  double least_trunc=GSL_DBL_MAX;
122  double result_least_trunc;
123 
124  /* Calculate specified minimum number of terms. No convergence
125  tests are made, and no truncation information is stored. */
126 
127  for (n=0; n < min_terms; n++) {
128  const double t=array[n];
129 
130  result_nm1=result_n;
131  levin_utrunc_step(t,n,result_n);
132  }
133 
134  /* Assume the result after the minimum calculation is the best. */
135 
136  result_least_trunc=result_n;
137 
138  /* Calculate up to maximum number of terms. Check truncation
139  condition. */
140 
141  for (; n <= nmax; n++) {
142  const double t=array[n];
143 
144  result_nm1=result_n;
145  levin_utrunc_step(t,n,result_n);
146 
147  /* Compute the truncation error directly */
148 
149  actual_trunc_nm1=actual_trunc_n;
150  actual_trunc_n=fabs(result_n-result_nm1);
151 
152  /* Average results to make a more reliable estimate of the
153  real truncation error */
154 
155  trunc_nm1=trunc_n;
156  trunc_n=0.5*(actual_trunc_n+actual_trunc_nm1);
157 
158  /* Determine if we are in the convergence region. */
159 
160  better=(trunc_n < trunc_nm1 || trunc_n < SMALL*fabs(result_n));
161  converging=converging || (better && before);
162  before=better;
163 
164  if (converging) {
165  if (trunc_n < least_trunc) {
166  /* Found a low truncation point in the convergence
167  region. Save it. */
168 
169  least_trunc=trunc_n;
170  result_least_trunc=result_n;
171  }
172 
173  if (fabs(trunc_n/result_n) < 10.0*GSL_MACH_EPS)
174  break;
175  }
176  }
177 
178  if (converging) {
179 
180  /* Stopped in the convergence region. Return result and
181  error estimate. */
182  sum_accel=result_least_trunc;
183  abserr_trunc=least_trunc;
184  wt->terms_used=n;
185  return sum_accel;
186 
187  } else {
188 
189  /* Never reached the convergence region. Use the last
190  calculated values. */
191  sum_accel=result_n;
192  abserr_trunc=trunc_n;
193  wt->terms_used=n;
194  return sum_accel;
195  }
196  }
197 
198 
199  return sum_accel;
200  }
201 
202  /** \brief Return the accelerated sum of the series with
203  an accurate error estimate
204 
205  The input vector \c x should be an array with \c n values
206  from <tt>x[0]</tt> to <tt>x[n-1]</tt>.
207  */
208  template<class vec_t>
209  double series_accel_err(size_t na, vec_t &array, double &abserr) {
210 
211  if (na!=size) set_size(na);
212 
213  double sum_accel;
214  size_t min_terms=0;
215  size_t max_terms=size-1;
216 
217  /* Ignore any trailing zeros in the array */
218  size_t size2=size;
219 
220  while (size2 > 0 && array[size2-1]==0) {
221  size2--;
222  }
223 
224  if (size2==0) {
225 
226  sum_accel=0.0;
227  abserr=0.0;
228  w->sum_plain=0.0;
229  w->terms_used=0;
230  return sum_accel;
231 
232  } else if (size2==1) {
233 
234  sum_accel=array[0];
235  abserr=0.0;
236  w->sum_plain=array[0];
237  w->terms_used=1;
238  return sum_accel;
239 
240  } else {
241 
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;
248  double variance=0;
249  size_t n;
250  unsigned int i;
251  int better=0;
252  int before=0;
253  int converging=0;
254  double least_trunc=GSL_DBL_MAX;
255  double least_trunc_noise=GSL_DBL_MAX;
256  double least_trunc_result;
257 
258  /* Calculate specified minimum number of terms. No convergence
259  tests are made, and no truncation information is stored. */
260 
261  for (n=0; n < min_terms; n++) {
262  const double t=array[n];
263  result_nm1=result_n;
264  levin_u_step(t,n,nmax,result_n);
265  }
266 
267  least_trunc_result=result_n;
268 
269  variance=0;
270  for (i=0; i < n; i++) {
271  double dn=w->dsum[i]*GSL_MACH_EPS*array[i];
272  variance += dn*dn;
273  }
274  noise_n=std::sqrt(variance);
275 
276  /* Calculate up to maximum number of terms. Check truncation
277  condition. */
278 
279  for (; n <= nmax; n++) {
280  const double t=array[n];
281 
282  result_nm1=result_n;
283  levin_u_step(t,n,nmax,result_n);
284 
285  /* Compute the truncation error directly */
286 
287  actual_trunc_nm1=actual_trunc_n;
288  actual_trunc_n=fabs(result_n-result_nm1);
289 
290  /* Average results to make a more reliable estimate of the
291  real truncation error */
292 
293  trunc_nm1=trunc_n;
294  trunc_n=0.5*(actual_trunc_n+actual_trunc_nm1);
295 
296  noise_nm1=noise_n;
297  variance=0;
298 
299  for (i=0; i <= n; i++) {
300  double dn=w->dsum[i]*GSL_MACH_EPS*array[i];
301  variance += dn*dn;
302  }
303 
304  noise_n=std::sqrt(variance);
305 
306  /* Determine if we are in the convergence region. */
307 
308  better=(trunc_n < trunc_nm1 ||
309  trunc_n < SMALL*fabs(result_n));
310  converging=converging || (better && before);
311  before=better;
312 
313  if (converging) {
314  if (trunc_n < least_trunc) {
315  /* Found a low truncation point in the convergence
316  region. Save it. */
317 
318  least_trunc_result=result_n;
319  least_trunc=trunc_n;
320  least_trunc_noise=noise_n;
321  }
322 
323  if (noise_n > trunc_n/3.0) {
324  break;
325  }
326 
327  if (trunc_n < 10.0*GSL_MACH_EPS*fabs(result_n)) {
328  break;
329  }
330  }
331 
332  }
333 
334  if (converging) {
335 
336  /* Stopped in the convergence region. Return result and
337  error estimate. */
338  sum_accel=least_trunc_result;
339  abserr=GSL_MAX_DBL(least_trunc,least_trunc_noise);
340  w->terms_used=n;
341  return sum_accel;
342 
343  } else {
344 
345  /* Never reached the convergence region. Use the last
346  calculated values. */
347  sum_accel=result_n;
348  abserr=GSL_MAX_DBL(trunc_n,noise_n);
349  w->terms_used=n;
350  return sum_accel;
351  }
352  }
353 
354  return sum_accel;
355  }
356 
357  /** \brief Set the number of terms */
358  void set_size(size_t new_size);
359 
360 #ifndef DOXYGEN_NO_O2NS
361 
362  protected:
363 
364  /** \brief An internal function reducing two matrix indices, i and j,
365  to index of a single array
366  */
367  size_t series_index(size_t i, size_t j, size_t nmax) {
368  return i*(nmax+1)+j;
369  }
370 
371  /** \brief Perform a step
372  */
373  int levin_u_step(const double term, const size_t n,
374  const size_t nmax, double &sum_accel);
375 
376  /** \brief Perform a step
377  */
378  int levin_utrunc_step(const double term, const size_t n,
379  double &sum_accel);
380 
381  /// The GSL workspace
382  gsl_sum_levin_u_workspace *w;
383 
384  /// The GSL workspace
385  gsl_sum_levin_utrunc_workspace *wt;
386 
387  /// The workspace size
388  size_t size;
389 
390 #endif
391 
392  };
393 
394 #ifndef DOXYGEN_NO_O2NS
395 }
396 #endif
397 
398 #endif
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
double series_accel_err(size_t na, vec_t &array, double &abserr)
Return the accelerated sum of the series with an accurate error estimate.
Definition: series_acc.h:209
gsl_sum_levin_utrunc_workspace * wt
The GSL workspace.
Definition: series_acc.h:385
double series_accel(size_t na, vec_t &array, double &abserr_trunc)
Return the accelerated sum of the series with a simple error estimate.
Definition: series_acc.h:86
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.
Definition: series_acc.h:388
Series acceleration by Levin u-transform (GSL)
Definition: series_acc.h:69
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.
Definition: series_acc.h:367
gsl_sum_levin_u_workspace * w
The GSL workspace.
Definition: series_acc.h:382

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).