htswanalysis/MACS/lib/gsl/gsl-1.11/sum/levin_utrunc.c

   1 /* sum/levin_utrunc.c
   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 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_test.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sum.h>
  27
  28 int
  29 gsl_sum_levin_utrunc_accel (const double *array,
  30                             const size_t array_size,
  31                             gsl_sum_levin_utrunc_workspace * w,
  32                             double *sum_accel, double *abserr_trunc)
  33 {
  34   return gsl_sum_levin_utrunc_minmax (array, array_size,
  35                                       0, array_size - 1,
  36                                       w, sum_accel, abserr_trunc);
  37 }
  38
  39
  40 int
  41 gsl_sum_levin_utrunc_minmax (const double *array,
  42                              const size_t array_size,
  43                              const size_t min_terms,
  44                              const size_t max_terms,
  45                              gsl_sum_levin_utrunc_workspace * w,
  46                              double *sum_accel, double *abserr_trunc)
  47 {
  48   if (array_size == 0)
  49     {
  50       *sum_accel = 0.0;
  51       *abserr_trunc = 0.0;
  52       w->sum_plain = 0.0;
  53       w->terms_used = 0;
  54       return GSL_SUCCESS;
  55     }
  56   else if (array_size == 1)
  57     {
  58       *sum_accel = array[0];
  59       *abserr_trunc = GSL_POSINF;
  60       w->sum_plain = array[0];
  61       w->terms_used = 1;
  62       return GSL_SUCCESS;
  63     }
  64   else
  65     {
  66       const double SMALL = 0.01;
  67       const size_t nmax = GSL_MAX (max_terms, array_size) - 1;
  68       double trunc_n = 0.0, trunc_nm1 = 0.0;
  69       double actual_trunc_n = 0.0, actual_trunc_nm1 = 0.0;
  70       double result_n = 0.0, result_nm1 = 0.0;
  71       size_t n;
  72       int better = 0;
  73       int before = 0;
  74       int converging = 0;
  75       double least_trunc = GSL_DBL_MAX;
  76       double result_least_trunc;
  77
  78       /* Calculate specified minimum number of terms. No convergence
  79          tests are made, and no truncation information is stored. */
  80
  81       for (n = 0; n < min_terms; n++)
  82         {
  83           const double t = array[n];
  84
  85           result_nm1 = result_n;
  86           gsl_sum_levin_utrunc_step (t, n, w, &result_n);
  87         }
  88
  89       /* Assume the result after the minimum calculation is the best. */
  90
  91       result_least_trunc = result_n;
  92
  93       /* Calculate up to maximum number of terms. Check truncation
  94          condition. */
  95
  96       for (; n <= nmax; n++)
  97         {
  98           const double t = array[n];
  99
 100           result_nm1 = result_n;
 101           gsl_sum_levin_utrunc_step (t, n, w, &result_n);
 102
 103           /* Compute the truncation error directly */
 104
 105           actual_trunc_nm1 = actual_trunc_n;
 106           actual_trunc_n = fabs (result_n - result_nm1);
 107
 108           /* Average results to make a more reliable estimate of the
 109              real truncation error */
 110
 111           trunc_nm1 = trunc_n;
 112           trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1);
 113
 114           /* Determine if we are in the convergence region. */
 115
 116           better = (trunc_n < trunc_nm1 || trunc_n < SMALL * fabs (result_n));
 117           converging = converging || (better && before);
 118           before = better;
 119
 120           if (converging)
 121             {
 122               if (trunc_n < least_trunc)
 123                 {
 124                   /* Found a low truncation point in the convergence
 125                      region. Save it. */
 126
 127                   least_trunc = trunc_n;
 128                   result_least_trunc = result_n;
 129                 }
 130
 131               if (fabs (trunc_n / result_n) < 10.0 * GSL_MACH_EPS)
 132                 break;
 133             }
 134         }
 135
 136       if (converging)
 137         {
 138           /* Stopped in the convergence region. Return result and
 139              error estimate. */
 140
 141           *sum_accel = result_least_trunc;
 142           *abserr_trunc = least_trunc;
 143           w->terms_used = n;
 144           return GSL_SUCCESS;
 145         }
 146       else
 147         {
 148           /* Never reached the convergence region. Use the last
 149              calculated values. */
 150
 151           *sum_accel = result_n;
 152           *abserr_trunc = trunc_n;
 153           w->terms_used = n;
 154           return GSL_SUCCESS;
 155         }
 156     }
 157 }
 158
 159 int
 160 gsl_sum_levin_utrunc_step (const double term,
 161                            const size_t n,
 162                            gsl_sum_levin_utrunc_workspace * w, double *sum_accel)
 163 {
 164   if (term == 0.0)
 165     {
 166       /* This is actually harmless when treated in this way. A term
 167          which is exactly zero is simply ignored; the state is not
 168          changed. We return GSL_EZERODIV as an indicator that this
 169          occured. */
 170
 171       return GSL_EZERODIV;
 172     }
 173   else if (n == 0)
 174     {
 175       *sum_accel = term;
 176       w->sum_plain = term;
 177       w->q_den[0] = 1.0 / term;
 178       w->q_num[0] = 1.0;
 179       return GSL_SUCCESS;
 180     }
 181   else
 182     {
 183       double factor = 1.0;
 184       double ratio = (double) n / (n + 1.0);
 185       int j;
 186
 187       w->sum_plain += term;
 188       w->q_den[n] = 1.0 / (term * (n + 1.0) * (n + 1.0));
 189       w->q_num[n] = w->sum_plain * w->q_den[n];
 190
 191       for (j = n - 1; j >= 0; j--)
 192         {
 193           double c = factor * (j + 1) / (n + 1);
 194           factor *= ratio;
 195           w->q_den[j] = w->q_den[j + 1] - c * w->q_den[j];
 196           w->q_num[j] = w->q_num[j + 1] - c * w->q_num[j];
 197         }
 198
 199       *sum_accel = w->q_num[0] / w->q_den[0];
 200       return GSL_SUCCESS;
 201     }
 202 }