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

   1 /* integration/qk.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 <float.h>
  22 #include <math.h>
  23 #include <gsl/gsl_integration.h>
  24 #include "err.c"
  25
  26 void
  27 gsl_integration_qk (const int n,
  28                     const double xgk[], const double wg[], const double wgk[],
  29                     double fv1[], double fv2[],
  30                     const gsl_function * f, double a, double b,
  31                     double *result, double *abserr,
  32                     double *resabs, double *resasc)
  33 {
  34
  35   const double center = 0.5 * (a + b);
  36   const double half_length = 0.5 * (b - a);
  37   const double abs_half_length = fabs (half_length);
  38   const double f_center = GSL_FN_EVAL (f, center);
  39
  40   double result_gauss = 0;
  41   double result_kronrod = f_center * wgk[n - 1];
  42
  43   double result_abs = fabs (result_kronrod);
  44   double result_asc = 0;
  45   double mean = 0, err = 0;
  46
  47   int j;
  48
  49   if (n % 2 == 0)
  50     {
  51       result_gauss = f_center * wg[n / 2 - 1];
  52     }
  53
  54   for (j = 0; j < (n - 1) / 2; j++)
  55     {
  56       const int jtw = j * 2 + 1;        /* j=1,2,3 jtw=2,4,6 */
  57       const double abscissa = half_length * xgk[jtw];
  58       const double fval1 = GSL_FN_EVAL (f, center - abscissa);
  59       const double fval2 = GSL_FN_EVAL (f, center + abscissa);
  60       const double fsum = fval1 + fval2;
  61       fv1[jtw] = fval1;
  62       fv2[jtw] = fval2;
  63       result_gauss += wg[j] * fsum;
  64       result_kronrod += wgk[jtw] * fsum;
  65       result_abs += wgk[jtw] * (fabs (fval1) + fabs (fval2));
  66     }
  67
  68   for (j = 0; j < n / 2; j++)
  69     {
  70       int jtwm1 = j * 2;
  71       const double abscissa = half_length * xgk[jtwm1];
  72       const double fval1 = GSL_FN_EVAL (f, center - abscissa);
  73       const double fval2 = GSL_FN_EVAL (f, center + abscissa);
  74       fv1[jtwm1] = fval1;
  75       fv2[jtwm1] = fval2;
  76       result_kronrod += wgk[jtwm1] * (fval1 + fval2);
  77       result_abs += wgk[jtwm1] * (fabs (fval1) + fabs (fval2));
  78     };
  79
  80   mean = result_kronrod * 0.5;
  81
  82   result_asc = wgk[n - 1] * fabs (f_center - mean);
  83
  84   for (j = 0; j < n - 1; j++)
  85     {
  86       result_asc += wgk[j] * (fabs (fv1[j] - mean) + fabs (fv2[j] - mean));
  87     }
  88
  89   /* scale by the width of the integration region */
  90
  91   err = (result_kronrod - result_gauss) * half_length;
  92
  93   result_kronrod *= half_length;
  94   result_abs *= abs_half_length;
  95   result_asc *= abs_half_length;
  96
  97   *result = result_kronrod;
  98   *resabs = result_abs;
  99   *resasc = result_asc;
 100   *abserr = rescale_error (err, result_abs, result_asc);
 101
 102 }