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

   1 /* integration/qc25s.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_qaws_params
  21 {
  22   gsl_function *function;
  23   double a;
  24   double b;
  25   gsl_integration_qaws_table *table;
  26 };
  27
  28 static double fn_qaws (double t, void *params);
  29 static double fn_qaws_L (double x, void *params);
  30 static double fn_qaws_R (double x, void *params);
  31
  32 static void
  33 compute_result (const double * r, const double * cheb12, const double * cheb24,
  34                 double * result12, double * result24);
  35
  36
  37 static void
  38 qc25s (gsl_function * f, double a, double b, double a1, double b1,
  39        gsl_integration_qaws_table * t,
  40        double *result, double *abserr, int *err_reliable);
  41
  42 static void
  43 qc25s (gsl_function * f, double a, double b, double a1, double b1,
  44        gsl_integration_qaws_table * t,
  45        double *result, double *abserr, int *err_reliable)
  46 {
  47   gsl_function weighted_function;
  48   struct fn_qaws_params fn_params;
  49
  50   fn_params.function = f;
  51   fn_params.a = a;
  52   fn_params.b = b;
  53   fn_params.table = t;
  54
  55   weighted_function.params = &fn_params;
  56
  57   if (a1 == a && (t->alpha != 0.0 || t->mu != 0))
  58     {
  59       double cheb12[13], cheb24[25];
  60
  61       double factor = pow(0.5 * (b1 - a1), t->alpha + 1.0);
  62
  63       weighted_function.function = &fn_qaws_R;
  64
  65       gsl_integration_qcheb (&weighted_function, a1, b1, cheb12, cheb24);
  66
  67       if (t->mu == 0)
  68         {
  69           double res12 = 0, res24 = 0;
  70           double u = factor;
  71
  72           compute_result (t->ri, cheb12, cheb24, &res12, &res24);
  73
  74           *result = u * res24;
  75           *abserr = fabs(u * (res24 - res12));
  76         }
  77       else
  78         {
  79           double res12a = 0, res24a = 0;
  80           double res12b = 0, res24b = 0;
  81
  82           double u = factor * log(b1 - a1);
  83           double v = factor;
  84
  85           compute_result (t->ri, cheb12, cheb24, &res12a, &res24a);
  86           compute_result (t->rg, cheb12, cheb24, &res12b, &res24b);
  87
  88           *result = u * res24a + v * res24b;
  89           *abserr = fabs(u * (res24a - res12a)) + fabs(v * (res24b - res12b));
  90         }
  91
  92       *err_reliable = 0;
  93
  94       return;
  95     }
  96   else if (b1 == b && (t->beta != 0.0 || t->nu != 0))
  97     {
  98       double cheb12[13], cheb24[25];
  99       double factor = pow(0.5 * (b1 - a1), t->beta + 1.0);
 100
 101       weighted_function.function = &fn_qaws_L;
 102
 103       gsl_integration_qcheb (&weighted_function, a1, b1, cheb12, cheb24);
 104
 105       if (t->nu == 0)
 106         {
 107           double res12 = 0, res24 = 0;
 108           double u = factor;
 109
 110           compute_result (t->rj, cheb12, cheb24, &res12, &res24);
 111
 112           *result = u * res24;
 113           *abserr = fabs(u * (res24 - res12));
 114         }
 115       else
 116         {
 117           double res12a = 0, res24a = 0;
 118           double res12b = 0, res24b = 0;
 119
 120           double u = factor * log(b1 - a1);
 121           double v = factor;
 122
 123           compute_result (t->rj, cheb12, cheb24, &res12a, &res24a);
 124           compute_result (t->rh, cheb12, cheb24, &res12b, &res24b);
 125
 126           *result = u * res24a + v * res24b;
 127           *abserr = fabs(u * (res24a - res12a)) + fabs(v * (res24b - res12b));
 128         }
 129
 130       *err_reliable = 0;
 131
 132       return;
 133     }
 134   else
 135     {
 136       double resabs, resasc;
 137
 138       weighted_function.function = &fn_qaws;
 139
 140       gsl_integration_qk15 (&weighted_function, a1, b1, result, abserr,
 141                             &resabs, &resasc);
 142
 143       if (*abserr == resasc)
 144         {
 145           *err_reliable = 0;
 146         }
 147       else
 148         {
 149           *err_reliable = 1;
 150         }
 151
 152       return;
 153     }
 154
 155 }
 156
 157 static double
 158 fn_qaws (double x, void *params)
 159 {
 160   struct fn_qaws_params *p = (struct fn_qaws_params *) params;
 161   gsl_function *f = p->function;
 162   gsl_integration_qaws_table *t = p->table;
 163
 164   double factor = 1.0;
 165
 166   if (t->alpha != 0.0)
 167     factor *= pow(x - p->a, t->alpha);
 168
 169   if (t->beta != 0.0)
 170     factor *= pow(p->b - x, t->beta);
 171
 172   if (t->mu == 1)
 173     factor *= log(x - p->a);
 174
 175   if (t->nu == 1)
 176     factor *= log(p->b - x);
 177
 178   return factor * GSL_FN_EVAL (f, x);
 179 }
 180
 181 static double
 182 fn_qaws_L (double x, void *params)
 183 {
 184   struct fn_qaws_params *p = (struct fn_qaws_params *) params;
 185   gsl_function *f = p->function;
 186   gsl_integration_qaws_table *t = p->table;
 187
 188   double factor = 1.0;
 189
 190   if (t->alpha != 0.0)
 191     factor *= pow(x - p->a, t->alpha);
 192
 193   if (t->mu == 1)
 194     factor *= log(x - p->a);
 195
 196   return factor * GSL_FN_EVAL (f, x);
 197 }
 198
 199 static double
 200 fn_qaws_R (double x, void *params)
 201 {
 202   struct fn_qaws_params *p = (struct fn_qaws_params *) params;
 203   gsl_function *f = p->function;
 204   gsl_integration_qaws_table *t = p->table;
 205
 206   double factor = 1.0;
 207
 208   if (t->beta != 0.0)
 209     factor *= pow(p->b - x, t->beta);
 210
 211   if (t->nu == 1)
 212     factor *= log(p->b - x);
 213
 214   return factor * GSL_FN_EVAL (f, x);
 215 }
 216
 217
 218 static void
 219 compute_result (const double * r, const double * cheb12, const double * cheb24,
 220                 double * result12, double * result24)
 221 {
 222   size_t i;
 223   double res12 = 0;
 224   double res24 = 0;
 225
 226   for (i = 0; i < 13; i++)
 227     {
 228       res12 += r[i] * cheb12[i];
 229     }
 230
 231   for (i = 0; i < 25; i++)
 232     {
 233       res24 += r[i] * cheb24[i];
 234     }
 235
 236   *result12 = res12;
 237   *result24 = res24;
 238 }