/* integration/romberg.c

 * 

 * Copyright (C) 2018 Patrick Alken

 * 

 * 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.

 */

/* the code in this module performs Romberg integration */

#include <config.h>

#include <stdlib.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_integration.h>

#define ROMBERG_PRINT_ROW(i, r) \

  do { \

    size_t jj; \

    fprintf(stderr, "R[%zu] = ", i); \

    for (jj = 0; jj <= i; ++jj) \

      fprintf(stderr, "%.8e ", r[jj]); \

    fprintf(stderr, "\n"); \

  } while (0)

gsl_integration_romberg_workspace *

gsl_integration_romberg_alloc(const size_t n)

{

  gsl_integration_romberg_workspace *w;

  /* check inputs */

  if (n < 1)

    {

      GSL_ERROR_VAL ("workspace size n must be at least 1", GSL_EDOM, 0);

    }

  w = calloc(1, sizeof(gsl_integration_romberg_workspace));

  if (w == NULL)

    {

      GSL_ERROR_VAL ("unable to allocate workspace", GSL_ENOMEM, 0);

    }

  /* ceiling on n, since the number of points is 2^n + 1 */

  w->n = GSL_MIN(n, 30);

  w->work1 = malloc(w->n * sizeof(double));

  if (w->work1 == NULL)

    {

      gsl_integration_romberg_free(w);

      GSL_ERROR_VAL ("unable to allocate previous row", GSL_ENOMEM, 0);

    }

  w->work2 = malloc(w->n * sizeof(double));

  if (w->work2 == NULL)

    {

      gsl_integration_romberg_free(w);

      GSL_ERROR_VAL ("unable to allocate current row", GSL_ENOMEM, 0);

    }

  return w;

}

void

gsl_integration_romberg_free(gsl_integration_romberg_workspace * w)

{

  if (w->work1)

    free(w->work1);

  if (w->work2)

    free(w->work2);

  free(w);

}

int

gsl_integration_romberg(const gsl_function * f, const double a, const double b,

                        const double epsabs, const double epsrel, double * result,

                        size_t * neval, gsl_integration_romberg_workspace * w)

{

  if (epsabs < 0.0)

    {

      GSL_ERROR("epsabs must be non-negative", GSL_EDOM);

    }

  else if (epsrel < 0.0)

    {

      GSL_ERROR("epsrel must be non-negative", GSL_EDOM);

    }

  else

    {

      const size_t n = w->n;

      double *Rp = &(w->work1[0]); /* previous row */

      double *Rc = &(w->work2[0]); /* current row */

      double *Rtmp;

      double h = 0.5 * (b - a);    /* step size */

      size_t i;

      /* R(0,0) */

      Rp[0] = h * (GSL_FN_EVAL(f, a) + GSL_FN_EVAL(f, b));

      *neval = 2;

      /*ROMBERG_PRINT_ROW((size_t) 0, Rp);*/

      for (i = 1; i < n; ++i)

        {

          size_t j;

          double sum = 0.0;

          double err;

          double four_j;         /* 4^j */

          size_t two_i = 1 << i; /* 2^i */

          for (j = 1; j < two_i; j += 2)

            {

              sum += GSL_FN_EVAL(f, a + j * h);

              ++(*neval);

            }

          /* R(i,0) */

          Rc[0] = sum * h + 0.5 * Rp[0];

          four_j = 4.0;

          for (j = 1; j <= i; ++j)

            {

              Rc[j] = (four_j * Rc[j - 1] - Rp[j - 1]) / (four_j - 1.0);

              four_j *= 4.0;

            }

          /*ROMBERG_PRINT_ROW(i, Rc);*/

          /*

           * compute difference between current and previous result and

           * check for convergence

           */

          err = fabs(Rc[i] - Rp[i - 1]);

          if ((err < epsabs) || (err < epsrel * fabs(Rc[i])))

            {

              *result = Rc[i];

              return GSL_SUCCESS;

            }

          /* swap Rp and Rc */

          Rtmp = Rp;

          Rp = Rc;

          Rc = Rtmp;

          h *= 0.5;

        }

      /* did not converge - return best guess */

      *result = Rp[n - 1];

      return GSL_EMAXITER;

    }

}