/* integration/util.c

 * 

 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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.

 */

static inline

void update (gsl_integration_workspace * workspace,

                 double a1, double b1, double area1, double error1,

                 double a2, double b2, double area2, double error2);

static inline void

retrieve (const gsl_integration_workspace * workspace, 

          double * a, double * b, double * r, double * e);

static inline

void update (gsl_integration_workspace * workspace,

             double a1, double b1, double area1, double error1,

             double a2, double b2, double area2, double error2)

{

  double * alist = workspace->alist ;

  double * blist = workspace->blist ;

  double * rlist = workspace->rlist ;

  double * elist = workspace->elist ;

  size_t * level = workspace->level ;

  const size_t i_max = workspace->i ;

  const size_t i_new = workspace->size ;

  const size_t new_level = workspace->level[i_max] + 1;

  /* append the newly-created intervals to the list */

  if (error2 > error1)

    {

      alist[i_max] = a2;        /* blist[maxerr] is already == b2 */

      rlist[i_max] = area2;

      elist[i_max] = error2;

      level[i_max] = new_level;

      alist[i_new] = a1;

      blist[i_new] = b1;

      rlist[i_new] = area1;

      elist[i_new] = error1;

      level[i_new] = new_level;

    }

  else

    {

      blist[i_max] = b1;        /* alist[maxerr] is already == a1 */

      rlist[i_max] = area1;

      elist[i_max] = error1;

      level[i_max] = new_level;

      alist[i_new] = a2;

      blist[i_new] = b2;

      rlist[i_new] = area2;

      elist[i_new] = error2;

      level[i_new] = new_level;

    }

  workspace->size++;

  if (new_level > workspace->maximum_level)

    {

      workspace->maximum_level = new_level;

    }

  qpsrt (workspace) ;

}

static inline void

retrieve (const gsl_integration_workspace * workspace, 

          double * a, double * b, double * r, double * e)

{

  const size_t i = workspace->i;

  double * alist = workspace->alist;

  double * blist = workspace->blist;

  double * rlist = workspace->rlist;

  double * elist = workspace->elist;

  *a = alist[i] ;

  *b = blist[i] ;

  *r = rlist[i] ;

  *e = elist[i] ;

}

static inline double

sum_results (const gsl_integration_workspace * workspace);

static inline double

sum_results (const gsl_integration_workspace * workspace)

{

  const double * const rlist = workspace->rlist ;

  const size_t n = workspace->size;

  size_t k;

  double result_sum = 0;

  for (k = 0; k < n; k++)

    {

      result_sum += rlist[k];

    }

  return result_sum;

}

static inline int

subinterval_too_small (double a1, double a2, double b2);

static inline int

subinterval_too_small (double a1, double a2, double b2)

{

  const double e = GSL_DBL_EPSILON;

  const double u = GSL_DBL_MIN;

  double tmp = (1 + 100 * e) * (fabs (a2) + 1000 * u);

  int status = fabs (a1) <= tmp && fabs (b2) <= tmp;

  return status;

}