htswanalysis/MACS/lib/gsl/gsl-1.11/sum/gsl_sum.h

   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__ */