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

   1 /* integration/qc25c.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 fn_cauchy_params
  21 {
  22   gsl_function *function;
  23   double singularity;
  24 };
  25
  26 static double fn_cauchy (double t, void *params);
  27
  28 static void compute_moments (double cc, double *moment);
  29
  30 static void
  31 qc25c (gsl_function * f, double a, double b, double c,
  32        double *result, double *abserr, int *err_reliable);
  33
  34 static void
  35 qc25c (gsl_function * f, double a, double b, double c,
  36        double *result, double *abserr, int *err_reliable)
  37 {
  38   double cc = (2 * c - b - a) / (b - a);
  39
  40   if (fabs (cc) > 1.1)
  41     {
  42       double resabs, resasc;
  43
  44       gsl_function weighted_function;
  45       struct fn_cauchy_params fn_params;
  46
  47       fn_params.function = f;
  48       fn_params.singularity = c;
  49
  50       weighted_function.function = &fn_cauchy;
  51       weighted_function.params = &fn_params;
  52
  53       gsl_integration_qk15 (&weighted_function, a, b, result, abserr,
  54                             &resabs, &resasc);
  55
  56       if (*abserr == resasc)
  57         {
  58           *err_reliable = 0;
  59         }
  60       else
  61         {
  62           *err_reliable = 1;
  63         }
  64
  65       return;
  66     }
  67   else
  68     {
  69       double cheb12[13], cheb24[25], moment[25];
  70       double res12 = 0, res24 = 0;
  71       size_t i;
  72       gsl_integration_qcheb (f, a, b, cheb12, cheb24);
  73       compute_moments (cc, moment);
  74
  75       for (i = 0; i < 13; i++)
  76         {
  77           res12 += cheb12[i] * moment[i];
  78         }
  79
  80       for (i = 0; i < 25; i++)
  81         {
  82           res24 += cheb24[i] * moment[i];
  83         }
  84
  85       *result = res24;
  86       *abserr = fabs(res24 - res12) ;
  87       *err_reliable = 0;
  88
  89       return;
  90     }
  91 }
  92
  93 static double
  94 fn_cauchy (double x, void *params)
  95 {
  96   struct fn_cauchy_params *p = (struct fn_cauchy_params *) params;
  97   gsl_function *f = p->function;
  98   double c = p->singularity;
  99   return GSL_FN_EVAL (f, x) / (x - c);
 100 }
 101
 102 static void
 103 compute_moments (double cc, double *moment)
 104 {
 105   size_t k;
 106
 107   double a0 = log (fabs ((1.0 - cc) / (1.0 + cc)));
 108   double a1 = 2 + a0 * cc;
 109
 110   moment[0] = a0;
 111   moment[1] = a1;
 112
 113   for (k = 2; k < 25; k++)
 114     {
 115       double a2;
 116
 117       if ((k % 2) == 0)
 118         {
 119           a2 = 2.0 * cc * a1 - a0;
 120         }
 121       else
 122         {
 123           const double km1 = k - 1.0;
 124           a2 = 2.0 * cc * a1 - a0 - 4.0 / (km1 * km1 - 1.0);
 125         }
 126
 127       moment[k] = a2;
 128
 129       a0 = a1;
 130       a1 = a2;
 131     }
 132 }