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

   1 /* integration/qelg.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 struct extrapolation_table
  21   {
  22     size_t n;
  23     double rlist2[52];
  24     size_t nres;
  25     double res3la[3];
  26   };
  27
  28 static void
  29   initialise_table (struct extrapolation_table *table);
  30
  31 static void
  32   append_table (struct extrapolation_table *table, double y);
  33
  34 static void
  35 initialise_table (struct extrapolation_table *table)
  36 {
  37   table->n = 0;
  38   table->nres = 0;
  39 }
  40 #ifdef JUNK
  41 static void
  42 initialise_table (struct extrapolation_table *table, double y)
  43 {
  44   table->n = 0;
  45   table->rlist2[0] = y;
  46   table->nres = 0;
  47 }
  48 #endif
  49 static void
  50 append_table (struct extrapolation_table *table, double y)
  51 {
  52   size_t n;
  53   n = table->n;
  54   table->rlist2[n] = y;
  55   table->n++;
  56 }
  57
  58 /* static inline void
  59    qelg (size_t * n, double epstab[],
  60    double * result, double * abserr,
  61    double res3la[], size_t * nres); */
  62
  63 static inline void
  64   qelg (struct extrapolation_table *table, double *result, double *abserr);
  65
  66 static inline void
  67 qelg (struct extrapolation_table *table, double *result, double *abserr)
  68 {
  69   double *epstab = table->rlist2;
  70   double *res3la = table->res3la;
  71   const size_t n = table->n - 1;
  72
  73   const double current = epstab[n];
  74
  75   double absolute = GSL_DBL_MAX;
  76   double relative = 5 * GSL_DBL_EPSILON * fabs (current);
  77
  78   const size_t newelm = n / 2;
  79   const size_t n_orig = n;
  80   size_t n_final = n;
  81   size_t i;
  82
  83   const size_t nres_orig = table->nres;
  84
  85   *result = current;
  86   *abserr = GSL_DBL_MAX;
  87
  88   if (n < 2)
  89     {
  90       *result = current;
  91       *abserr = GSL_MAX_DBL (absolute, relative);
  92       return;
  93     }
  94
  95   epstab[n + 2] = epstab[n];
  96   epstab[n] = GSL_DBL_MAX;
  97
  98   for (i = 0; i < newelm; i++)
  99     {
 100       double res = epstab[n - 2 * i + 2];
 101       double e0 = epstab[n - 2 * i - 2];
 102       double e1 = epstab[n - 2 * i - 1];
 103       double e2 = res;
 104
 105       double e1abs = fabs (e1);
 106       double delta2 = e2 - e1;
 107       double err2 = fabs (delta2);
 108       double tol2 = GSL_MAX_DBL (fabs (e2), e1abs) * GSL_DBL_EPSILON;
 109       double delta3 = e1 - e0;
 110       double err3 = fabs (delta3);
 111       double tol3 = GSL_MAX_DBL (e1abs, fabs (e0)) * GSL_DBL_EPSILON;
 112
 113       double e3, delta1, err1, tol1, ss;
 114
 115       if (err2 <= tol2 && err3 <= tol3)
 116         {
 117           /* If e0, e1 and e2 are equal to within machine accuracy,
 118              convergence is assumed.  */
 119
 120           *result = res;
 121           absolute = err2 + err3;
 122           relative = 5 * GSL_DBL_EPSILON * fabs (res);
 123           *abserr = GSL_MAX_DBL (absolute, relative);
 124           return;
 125         }
 126
 127       e3 = epstab[n - 2 * i];
 128       epstab[n - 2 * i] = e1;
 129       delta1 = e1 - e3;
 130       err1 = fabs (delta1);
 131       tol1 = GSL_MAX_DBL (e1abs, fabs (e3)) * GSL_DBL_EPSILON;
 132
 133       /* If two elements are very close to each other, omit a part of
 134          the table by adjusting the value of n */
 135
 136       if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3)
 137         {
 138           n_final = 2 * i;
 139           break;
 140         }
 141
 142       ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
 143
 144       /* Test to detect irregular behaviour in the table, and
 145          eventually omit a part of the table by adjusting the value of
 146          n. */
 147
 148       if (fabs (ss * e1) <= 0.0001)
 149         {
 150           n_final = 2 * i;
 151           break;
 152         }
 153
 154       /* Compute a new element and eventually adjust the value of
 155          result. */
 156
 157       res = e1 + 1 / ss;
 158       epstab[n - 2 * i] = res;
 159
 160       {
 161         const double error = err2 + fabs (res - e2) + err3;
 162
 163         if (error <= *abserr)
 164           {
 165             *abserr = error;
 166             *result = res;
 167           }
 168       }
 169     }
 170
 171   /* Shift the table */
 172
 173   {
 174     const size_t limexp = 50 - 1;
 175
 176     if (n_final == limexp)
 177       {
 178         n_final = 2 * (limexp / 2);
 179       }
 180   }
 181
 182   if (n_orig % 2 == 1)
 183     {
 184       for (i = 0; i <= newelm; i++)
 185         {
 186           epstab[1 + i * 2] = epstab[i * 2 + 3];
 187         }
 188     }
 189   else
 190     {
 191       for (i = 0; i <= newelm; i++)
 192         {
 193           epstab[i * 2] = epstab[i * 2 + 2];
 194         }
 195     }
 196
 197   if (n_orig != n_final)
 198     {
 199       for (i = 0; i <= n_final; i++)
 200         {
 201           epstab[i] = epstab[n_orig - n_final + i];
 202         }
 203     }
 204
 205   table->n = n_final + 1;
 206
 207   if (nres_orig < 3)
 208     {
 209       res3la[nres_orig] = *result;
 210       *abserr = GSL_DBL_MAX;
 211     }
 212   else
 213     {                           /* Compute error estimate */
 214       *abserr = (fabs (*result - res3la[2]) + fabs (*result - res3la[1])
 215                  + fabs (*result - res3la[0]));
 216
 217       res3la[0] = res3la[1];
 218       res3la[1] = res3la[2];
 219       res3la[2] = *result;
 220     }
 221
 222   /* In QUADPACK the variable table->nres is incremented at the top of
 223      qelg, so it increases on every call. This leads to the array
 224      res3la being accessed when its elements are still undefined, so I
 225      have moved the update to this point so that its value more
 226      useful. */
 227
 228   table->nres = nres_orig + 1;
 229
 230   *abserr = GSL_MAX_DBL (*abserr, 5 * GSL_DBL_EPSILON * fabs (*result));
 231
 232   return;
 233 }