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

   1 /* integration/qng.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 #include <config.h>
  21 #include <math.h>
  22 #include <float.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_integration.h>
  26
  27 #include "err.c"
  28 #include "qng.h"
  29
  30 int
  31 gsl_integration_qng (const gsl_function *f,
  32                      double a, double b,
  33                      double epsabs, double epsrel,
  34                      double * result, double * abserr, size_t * neval)
  35 {
  36   double fv1[5], fv2[5], fv3[5], fv4[5];
  37   double savfun[21];  /* array of function values which have been computed */
  38   double res10, res21, res43, res87;    /* 10, 21, 43 and 87 point results */
  39   double result_kronrod, err ;
  40   double resabs; /* approximation to the integral of abs(f) */
  41   double resasc; /* approximation to the integral of abs(f-i/(b-a)) */
  42
  43   const double half_length =  0.5 * (b - a);
  44   const double abs_half_length = fabs (half_length);
  45   const double center = 0.5 * (b + a);
  46   const double f_center = GSL_FN_EVAL(f, center);
  47
  48   int k ;
  49
  50   if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
  51     {
  52       * result = 0;
  53       * abserr = 0;
  54       * neval = 0;
  55       GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
  56                  GSL_EBADTOL);
  57     };
  58
  59   /* Compute the integral using the 10- and 21-point formula. */
  60
  61   res10 = 0;
  62   res21 = w21b[5] * f_center;
  63   resabs = w21b[5] * fabs (f_center);
  64
  65   for (k = 0; k < 5; k++)
  66     {
  67       const double abscissa = half_length * x1[k];
  68       const double fval1 = GSL_FN_EVAL(f, center + abscissa);
  69       const double fval2 = GSL_FN_EVAL(f, center - abscissa);
  70       const double fval = fval1 + fval2;
  71       res10 += w10[k] * fval;
  72       res21 += w21a[k] * fval;
  73       resabs += w21a[k] * (fabs (fval1) + fabs (fval2));
  74       savfun[k] = fval;
  75       fv1[k] = fval1;
  76       fv2[k] = fval2;
  77     }
  78
  79   for (k = 0; k < 5; k++)
  80     {
  81       const double abscissa = half_length * x2[k];
  82       const double fval1 = GSL_FN_EVAL(f, center + abscissa);
  83       const double fval2 = GSL_FN_EVAL(f, center - abscissa);
  84       const double fval = fval1 + fval2;
  85       res21 += w21b[k] * fval;
  86       resabs += w21b[k] * (fabs (fval1) + fabs (fval2));
  87       savfun[k + 5] = fval;
  88       fv3[k] = fval1;
  89       fv4[k] = fval2;
  90     }
  91
  92   resabs *= abs_half_length ;
  93
  94   {
  95     const double mean = 0.5 * res21;
  96
  97     resasc = w21b[5] * fabs (f_center - mean);
  98
  99     for (k = 0; k < 5; k++)
 100       {
 101         resasc +=
 102           (w21a[k] * (fabs (fv1[k] - mean) + fabs (fv2[k] - mean))
 103           + w21b[k] * (fabs (fv3[k] - mean) + fabs (fv4[k] - mean)));
 104       }
 105     resasc *= abs_half_length ;
 106   }
 107
 108   result_kronrod = res21 * half_length;
 109
 110   err = rescale_error ((res21 - res10) * half_length, resabs, resasc) ;
 111
 112   /*   test for convergence. */
 113
 114   if (err < epsabs || err < epsrel * fabs (result_kronrod))
 115     {
 116       * result = result_kronrod ;
 117       * abserr = err ;
 118       * neval = 21;
 119       return GSL_SUCCESS;
 120     }
 121
 122   /* compute the integral using the 43-point formula. */
 123
 124   res43 = w43b[11] * f_center;
 125
 126   for (k = 0; k < 10; k++)
 127     {
 128       res43 += savfun[k] * w43a[k];
 129     }
 130
 131   for (k = 0; k < 11; k++)
 132     {
 133       const double abscissa = half_length * x3[k];
 134       const double fval = (GSL_FN_EVAL(f, center + abscissa)
 135                            + GSL_FN_EVAL(f, center - abscissa));
 136       res43 += fval * w43b[k];
 137       savfun[k + 10] = fval;
 138     }
 139
 140   /*  test for convergence */
 141
 142   result_kronrod = res43 * half_length;
 143   err = rescale_error ((res43 - res21) * half_length, resabs, resasc);
 144
 145   if (err < epsabs || err < epsrel * fabs (result_kronrod))
 146     {
 147       * result = result_kronrod ;
 148       * abserr = err ;
 149       * neval = 43;
 150       return GSL_SUCCESS;
 151     }
 152
 153   /* compute the integral using the 87-point formula. */
 154
 155   res87 = w87b[22] * f_center;
 156
 157   for (k = 0; k < 21; k++)
 158     {
 159       res87 += savfun[k] * w87a[k];
 160     }
 161
 162   for (k = 0; k < 22; k++)
 163     {
 164       const double abscissa = half_length * x4[k];
 165       res87 += w87b[k] * (GSL_FN_EVAL(f, center + abscissa)
 166                           + GSL_FN_EVAL(f, center - abscissa));
 167     }
 168
 169   /*  test for convergence */
 170
 171   result_kronrod = res87 * half_length ;
 172
 173   err = rescale_error ((res87 - res43) * half_length, resabs, resasc);
 174
 175   if (err < epsabs || err < epsrel * fabs (result_kronrod))
 176     {
 177       * result = result_kronrod ;
 178       * abserr = err ;
 179       * neval = 87;
 180       return GSL_SUCCESS;
 181     }
 182
 183   /* failed to converge */
 184
 185   * result = result_kronrod ;
 186   * abserr = err ;
 187   * neval = 87;
 188
 189   GSL_ERROR("failed to reach tolerance with highest-order rule", GSL_ETOL) ;
 190 }