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

   1 /* sum/levin_u.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 #include <config.h>
  21 #include <gsl/gsl_math.h>
  22 #include <gsl/gsl_test.h>
  23 #include <gsl/gsl_errno.h>
  24 #include <gsl/gsl_sum.h>
  25
  26 int
  27 gsl_sum_levin_u_accel (const double *array, const size_t array_size,
  28                        gsl_sum_levin_u_workspace * w,
  29                        double *sum_accel, double *abserr)
  30 {
  31   return gsl_sum_levin_u_minmax (array, array_size,
  32                                  0, array_size - 1, w, sum_accel, abserr);
  33 }
  34
  35 int
  36 gsl_sum_levin_u_minmax (const double *array, const size_t array_size,
  37                         const size_t min_terms, const size_t max_terms,
  38                         gsl_sum_levin_u_workspace * w,
  39                         double *sum_accel, double *abserr)
  40 {
  41   if (array_size == 0)
  42     {
  43       *sum_accel = 0.0;
  44       *abserr = 0.0;
  45       w->sum_plain = 0.0;
  46       w->terms_used = 0;
  47       return GSL_SUCCESS;
  48     }
  49   else if (array_size == 1)
  50     {
  51       *sum_accel = array[0];
  52       *abserr = 0.0;
  53       w->sum_plain = array[0];
  54       w->terms_used = 1;
  55       return GSL_SUCCESS;
  56     }
  57   else
  58     {
  59       const double SMALL = 0.01;
  60       const size_t nmax = GSL_MAX (max_terms, array_size) - 1;
  61       double noise_n = 0.0, noise_nm1 = 0.0;
  62       double trunc_n = 0.0, trunc_nm1 = 0.0;
  63       double actual_trunc_n = 0.0, actual_trunc_nm1 = 0.0;
  64       double result_n = 0.0, result_nm1 = 0.0;
  65       double variance = 0;
  66       size_t n;
  67       unsigned int i;
  68       int better = 0;
  69       int before = 0;
  70       int converging = 0;
  71       double least_trunc = GSL_DBL_MAX;
  72       double least_trunc_noise = GSL_DBL_MAX;
  73       double least_trunc_result;
  74
  75       /* Calculate specified minimum number of terms.  No convergence
  76          tests are made, and no truncation information is stored.  */
  77
  78       for (n = 0; n < min_terms; n++)
  79         {
  80           const double t = array[n];
  81           result_nm1 = result_n;
  82           gsl_sum_levin_u_step (t, n, nmax, w, &result_n);
  83         }
  84
  85       least_trunc_result = result_n;
  86
  87       variance = 0;
  88       for (i = 0; i < n; i++)
  89         {
  90           double dn = w->dsum[i] * GSL_MACH_EPS * array[i];
  91           variance += dn * dn;
  92         }
  93       noise_n = sqrt (variance);
  94
  95       /* Calculate up to maximum number of terms.  Check truncation
  96          condition.  */
  97
  98       for (; n <= nmax; n++)
  99         {
 100           const double t = array[n];
 101
 102           result_nm1 = result_n;
 103           gsl_sum_levin_u_step (t, n, nmax, w, &result_n);
 104
 105           /* Compute the truncation error directly */
 106
 107           actual_trunc_nm1 = actual_trunc_n;
 108           actual_trunc_n = fabs (result_n - result_nm1);
 109
 110           /* Average results to make a more reliable estimate of the
 111              real truncation error */
 112
 113           trunc_nm1 = trunc_n;
 114           trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1);
 115
 116           noise_nm1 = noise_n;
 117           variance = 0;
 118
 119           for (i = 0; i <= n; i++)
 120             {
 121               double dn = w->dsum[i] * GSL_MACH_EPS * array[i];
 122               variance += dn * dn;
 123             }
 124
 125           noise_n = sqrt (variance);
 126
 127           /* Determine if we are in the convergence region.  */
 128
 129           better = (trunc_n < trunc_nm1 || trunc_n < SMALL * fabs (result_n));
 130           converging = converging || (better && before);
 131           before = better;
 132
 133           if (converging)
 134             {
 135               if (trunc_n < least_trunc)
 136                 {
 137                   /* Found a low truncation point in the convergence
 138                      region. Save it. */
 139
 140                   least_trunc_result = result_n;
 141                   least_trunc = trunc_n;
 142                   least_trunc_noise = noise_n;
 143                 }
 144
 145               if (noise_n > trunc_n / 3.0)
 146                 break;
 147
 148               if (trunc_n < 10.0 * GSL_MACH_EPS * fabs (result_n))
 149                 break;
 150             }
 151
 152         }
 153
 154       if (converging)
 155         {
 156           /* Stopped in the convergence region.  Return result and
 157              error estimate.  */
 158
 159           *sum_accel = least_trunc_result;
 160           *abserr = GSL_MAX_DBL (least_trunc, least_trunc_noise);
 161           w->terms_used = n;
 162           return GSL_SUCCESS;
 163         }
 164       else
 165         {
 166           /* Never reached the convergence region.  Use the last
 167              calculated values.  */
 168
 169           *sum_accel = result_n;
 170           *abserr = GSL_MAX_DBL (trunc_n, noise_n);
 171           w->terms_used = n;
 172           return GSL_SUCCESS;
 173         }
 174     }
 175 }
 176
 177
 178 int
 179 gsl_sum_levin_u_step (const double term, const size_t n, const size_t nmax,
 180                       gsl_sum_levin_u_workspace * w, double *sum_accel)
 181 {
 182
 183 #define I(i,j) ((i)*(nmax+1) + (j))
 184
 185   if (n == 0)
 186     {
 187       *sum_accel = term;
 188       w->sum_plain = term;
 189
 190       w->q_den[0] = 1.0 / term;
 191       w->q_num[0] = 1.0;
 192
 193       w->dq_den[I (0, 0)] = -1.0 / (term * term);
 194       w->dq_num[I (0, 0)] = 0.0;
 195
 196       w->dsum[0] = 1.0;
 197
 198       return GSL_SUCCESS;
 199     }
 200   else
 201     {
 202       double result;
 203       double factor = 1.0;
 204       double ratio = (double) n / (n + 1.0);
 205       unsigned int i;
 206       int j;
 207
 208       w->sum_plain += term;
 209
 210       w->q_den[n] = 1.0 / (term * (n + 1.0) * (n + 1.0));
 211       w->q_num[n] = w->sum_plain * w->q_den[n];
 212
 213       for (i = 0; i < n; i++)
 214         {
 215           w->dq_den[I (i, n)] = 0;
 216           w->dq_num[I (i, n)] = w->q_den[n];
 217         }
 218
 219       w->dq_den[I (n, n)] = -w->q_den[n] / term;
 220       w->dq_num[I (n, n)] =
 221         w->q_den[n] + w->sum_plain * (w->dq_den[I (n, n)]);
 222
 223       for (j = n - 1; j >= 0; j--)
 224         {
 225           double c = factor * (j + 1) / (n + 1);
 226           factor *= ratio;
 227           w->q_den[j] = w->q_den[j + 1] - c * w->q_den[j];
 228           w->q_num[j] = w->q_num[j + 1] - c * w->q_num[j];
 229
 230           for (i = 0; i < n; i++)
 231             {
 232               w->dq_den[I (i, j)] =
 233                 w->dq_den[I (i, j + 1)] - c * w->dq_den[I (i, j)];
 234               w->dq_num[I (i, j)] =
 235                 w->dq_num[I (i, j + 1)] - c * w->dq_num[I (i, j)];
 236             }
 237
 238           w->dq_den[I (n, j)] = w->dq_den[I (n, j + 1)];
 239           w->dq_num[I (n, j)] = w->dq_num[I (n, j + 1)];
 240         }
 241
 242       result = w->q_num[0] / w->q_den[0];
 243
 244       *sum_accel = result;
 245
 246       for (i = 0; i <= n; i++)
 247         {
 248           w->dsum[i] =
 249             (w->dq_num[I (i, 0)] -
 250              result * w->dq_den[I (i, 0)]) / w->q_den[0];
 251         }
 252
 253       return GSL_SUCCESS;
 254     }
 255 }