htswanalysis/MACS/lib/gsl/gsl-1.11/integration/util.c

   1 /* integration/util.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 static inline
  21 void update (gsl_integration_workspace * workspace,
  22                  double a1, double b1, double area1, double error1,
  23                  double a2, double b2, double area2, double error2);
  24
  25 static inline void
  26 retrieve (const gsl_integration_workspace * workspace,
  27           double * a, double * b, double * r, double * e);
  28
  29
  30
  31 static inline
  32 void update (gsl_integration_workspace * workspace,
  33              double a1, double b1, double area1, double error1,
  34              double a2, double b2, double area2, double error2)
  35 {
  36   double * alist = workspace->alist ;
  37   double * blist = workspace->blist ;
  38   double * rlist = workspace->rlist ;
  39   double * elist = workspace->elist ;
  40   size_t * level = workspace->level ;
  41
  42   const size_t i_max = workspace->i ;
  43   const size_t i_new = workspace->size ;
  44
  45   const size_t new_level = workspace->level[i_max] + 1;
  46
  47   /* append the newly-created intervals to the list */
  48
  49   if (error2 > error1)
  50     {
  51       alist[i_max] = a2;        /* blist[maxerr] is already == b2 */
  52       rlist[i_max] = area2;
  53       elist[i_max] = error2;
  54       level[i_max] = new_level;
  55
  56       alist[i_new] = a1;
  57       blist[i_new] = b1;
  58       rlist[i_new] = area1;
  59       elist[i_new] = error1;
  60       level[i_new] = new_level;
  61     }
  62   else
  63     {
  64       blist[i_max] = b1;        /* alist[maxerr] is already == a1 */
  65       rlist[i_max] = area1;
  66       elist[i_max] = error1;
  67       level[i_max] = new_level;
  68
  69       alist[i_new] = a2;
  70       blist[i_new] = b2;
  71       rlist[i_new] = area2;
  72       elist[i_new] = error2;
  73       level[i_new] = new_level;
  74     }
  75
  76   workspace->size++;
  77
  78   if (new_level > workspace->maximum_level)
  79     {
  80       workspace->maximum_level = new_level;
  81     }
  82
  83   qpsrt (workspace) ;
  84 }
  85
  86 static inline void
  87 retrieve (const gsl_integration_workspace * workspace,
  88           double * a, double * b, double * r, double * e)
  89 {
  90   const size_t i = workspace->i;
  91   double * alist = workspace->alist;
  92   double * blist = workspace->blist;
  93   double * rlist = workspace->rlist;
  94   double * elist = workspace->elist;
  95
  96   *a = alist[i] ;
  97   *b = blist[i] ;
  98   *r = rlist[i] ;
  99   *e = elist[i] ;
 100 }
 101
 102 static inline double
 103 sum_results (const gsl_integration_workspace * workspace);
 104
 105 static inline double
 106 sum_results (const gsl_integration_workspace * workspace)
 107 {
 108   const double * const rlist = workspace->rlist ;
 109   const size_t n = workspace->size;
 110
 111   size_t k;
 112   double result_sum = 0;
 113
 114   for (k = 0; k < n; k++)
 115     {
 116       result_sum += rlist[k];
 117     }
 118
 119   return result_sum;
 120 }
 121
 122 static inline int
 123 subinterval_too_small (double a1, double a2, double b2);
 124
 125 static inline int
 126 subinterval_too_small (double a1, double a2, double b2)
 127 {
 128   const double e = GSL_DBL_EPSILON;
 129   const double u = GSL_DBL_MIN;
 130
 131   double tmp = (1 + 100 * e) * (fabs (a2) + 1000 * u);
 132
 133   int status = fabs (a1) <= tmp && fabs (b2) <= tmp;
 134
 135   return status;
 136 }
 137