/* sum/gsl_sum.h
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */


#ifndef __GSL_SUM_H__
#define __GSL_SUM_H__

#include <stdlib.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS          /* empty */
# define __END_DECLS            /* empty */
#endif

__BEGIN_DECLS

/*  Workspace for Levin U Transform with error estimation,
 *   
 *   size        = number of terms the workspace can handle
 *   sum_plain   = simple sum of series
 *   q_num       = backward diagonal of numerator; length = size
 *   q_den       = backward diagonal of denominator; length = size
 *   dq_num      = table of numerator derivatives; length = size**2
 *   dq_den      = table of denominator derivatives; length = size**2
 *   dsum        = derivative of sum wrt term i; length = size
 */

typedef struct
{
  size_t size;
  size_t i;                     /* position in array */
  size_t terms_used;            /* number of calls */
  double sum_plain;
  double *q_num;
  double *q_den;
  double *dq_num;
  double *dq_den;
  double *dsum;
}
gsl_sum_levin_u_workspace;

gsl_sum_levin_u_workspace *gsl_sum_levin_u_alloc (size_t n);
void gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * w);

/* Basic Levin-u acceleration method.
 *
 *   array       = array of series elements
 *   n           = size of array
 *   sum_accel   = result of summation acceleration
 *   err         = estimated error   
 *
 * See [Fessler et al., ACM TOMS 9, 346 (1983) and TOMS-602]
 */

int gsl_sum_levin_u_accel (const double *array,
                           const size_t n,
                           gsl_sum_levin_u_workspace * w,
                           double *sum_accel, double *abserr);

/* Basic Levin-u acceleration method with constraints on the terms
 * used,
 *
 *   array       = array of series elements
 *   n           = size of array
 *   min_terms   = minimum number of terms to sum
 *   max_terms   = maximum number of terms to sum
 *   sum_accel   = result of summation acceleration
 *   err         = estimated error   
 *
 * See [Fessler et al., ACM TOMS 9, 346 (1983) and TOMS-602] 
 */

int gsl_sum_levin_u_minmax (const double *array,
                            const size_t n,
                            const size_t min_terms,
                            const size_t max_terms,
                            gsl_sum_levin_u_workspace * w,
                            double *sum_accel, double *abserr);

/* Basic Levin-u step w/o reference to the array of terms.
 * We only need to specify the value of the current term
 * to execute the step. See TOMS-745.
 *
 * sum = t0 + ... + t_{n-1} + term;  term = t_{n}
 *
 *   term   = value of the series term to be added
 *   n      = position of term in series (starting from 0)
 *   sum_accel = result of summation acceleration
 *   sum_plain = simple sum of series
 */

int
gsl_sum_levin_u_step (const double term,
                      const size_t n,
                      const size_t nmax,
                      gsl_sum_levin_u_workspace * w, 
                      double *sum_accel);

/* The following functions perform the same calculation without
   estimating the errors. They require O(N) storage instead of O(N^2).
   This may be useful for summing many similar series where the size
   of the error has already been estimated reliably and is not
   expected to change.  */

typedef struct
{
  size_t size;
  size_t i;                     /* position in array */
  size_t terms_used;            /* number of calls */
  double sum_plain;
  double *q_num;
  double *q_den;
  double *dsum;
}
gsl_sum_levin_utrunc_workspace;

gsl_sum_levin_utrunc_workspace *gsl_sum_levin_utrunc_alloc (size_t n);
void gsl_sum_levin_utrunc_free (gsl_sum_levin_utrunc_workspace * w);

int gsl_sum_levin_utrunc_accel (const double *array,
                                const size_t n,
                                gsl_sum_levin_utrunc_workspace * w,
                                double *sum_accel, double *abserr_trunc);

int gsl_sum_levin_utrunc_minmax (const double *array,
                                 const size_t n,
                                 const size_t min_terms,
                                 const size_t max_terms,
                                 gsl_sum_levin_utrunc_workspace * w,
                                 double *sum_accel, double *abserr_trunc);

int gsl_sum_levin_utrunc_step (const double term,
                               const size_t n,
                               gsl_sum_levin_utrunc_workspace * w, 
                               double *sum_accel);

__END_DECLS

#endif /* __GSL_SUM_H__ */