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

   1 /* integration/qags.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 <stdlib.h>
  22 #include <gsl/gsl_errno.h>
  23 #include <gsl/gsl_integration.h>
  24
  25 #include "initialise.c"
  26 #include "set_initial.c"
  27 #include "qpsrt.c"
  28 #include "util.c"
  29 #include "reset.c"
  30 #include "qpsrt2.c"
  31 #include "qelg.c"
  32 #include "positivity.c"
  33
  34 static int qags (const gsl_function * f, const double a, const double
  35   b, const double epsabs, const double epsrel, const size_t limit,
  36   gsl_integration_workspace * workspace, double *result, double *abserr,
  37   gsl_integration_rule * q);
  38
  39 int
  40 gsl_integration_qags (const gsl_function *f,
  41                       double a, double b,
  42                       double epsabs, double epsrel, size_t limit,
  43                       gsl_integration_workspace * workspace,
  44                       double * result, double * abserr)
  45 {
  46   int status = qags (f, a, b, epsabs, epsrel, limit,
  47                      workspace,
  48                      result, abserr,
  49                      &gsl_integration_qk21) ;
  50   return status ;
  51 }
  52
  53 /* QAGI: evaluate an integral over an infinite range using the
  54    transformation
  55
  56    integrate(f(x),-Inf,Inf) = integrate((f((1-t)/t) + f(-(1-t)/t))/t^2,0,1)
  57
  58    */
  59
  60 static double i_transform (double t, void *params);
  61
  62 int
  63 gsl_integration_qagi (gsl_function * f,
  64                       double epsabs, double epsrel, size_t limit,
  65                       gsl_integration_workspace * workspace,
  66                       double *result, double *abserr)
  67 {
  68   int status;
  69
  70   gsl_function f_transform;
  71
  72   f_transform.function = &i_transform;
  73   f_transform.params = f;
  74
  75   status = qags (&f_transform, 0.0, 1.0,
  76                  epsabs, epsrel, limit,
  77                  workspace,
  78                  result, abserr,
  79                  &gsl_integration_qk15);
  80
  81   return status;
  82 }
  83
  84 static double
  85 i_transform (double t, void *params)
  86 {
  87   gsl_function *f = (gsl_function *) params;
  88   double x = (1 - t) / t;
  89   double y = GSL_FN_EVAL (f, x) + GSL_FN_EVAL (f, -x);
  90   return (y / t) / t;
  91 }
  92
  93
  94 /* QAGIL: Evaluate an integral over an infinite range using the
  95    transformation,
  96
  97    integrate(f(x),-Inf,b) = integrate(f(b-(1-t)/t)/t^2,0,1)
  98
  99    */
 100
 101 struct il_params { double b ; gsl_function * f ; } ;
 102
 103 static double il_transform (double t, void *params);
 104
 105 int
 106 gsl_integration_qagil (gsl_function * f,
 107                        double b,
 108                        double epsabs, double epsrel, size_t limit,
 109                        gsl_integration_workspace * workspace,
 110                        double *result, double *abserr)
 111 {
 112   int status;
 113
 114   gsl_function f_transform;
 115   struct il_params transform_params  ;
 116
 117   transform_params.b = b ;
 118   transform_params.f = f ;
 119
 120   f_transform.function = &il_transform;
 121   f_transform.params = &transform_params;
 122
 123   status = qags (&f_transform, 0.0, 1.0,
 124                  epsabs, epsrel, limit,
 125                  workspace,
 126                  result, abserr,
 127                  &gsl_integration_qk15);
 128
 129   return status;
 130 }
 131
 132 static double
 133 il_transform (double t, void *params)
 134 {
 135   struct il_params *p = (struct il_params *) params;
 136   double b = p->b;
 137   gsl_function * f = p->f;
 138   double x = b - (1 - t) / t;
 139   double y = GSL_FN_EVAL (f, x);
 140   return (y / t) / t;
 141 }
 142
 143 /* QAGIU: Evaluate an integral over an infinite range using the
 144    transformation
 145
 146    integrate(f(x),a,Inf) = integrate(f(a+(1-t)/t)/t^2,0,1)
 147
 148    */
 149
 150 struct iu_params { double a ; gsl_function * f ; } ;
 151
 152 static double iu_transform (double t, void *params);
 153
 154 int
 155 gsl_integration_qagiu (gsl_function * f,
 156                        double a,
 157                        double epsabs, double epsrel, size_t limit,
 158                        gsl_integration_workspace * workspace,
 159                        double *result, double *abserr)
 160 {
 161   int status;
 162
 163   gsl_function f_transform;
 164   struct iu_params transform_params  ;
 165
 166   transform_params.a = a ;
 167   transform_params.f = f ;
 168
 169   f_transform.function = &iu_transform;
 170   f_transform.params = &transform_params;
 171
 172   status = qags (&f_transform, 0.0, 1.0,
 173                  epsabs, epsrel, limit,
 174                  workspace,
 175                  result, abserr,
 176                  &gsl_integration_qk15);
 177
 178   return status;
 179 }
 180
 181 static double
 182 iu_transform (double t, void *params)
 183 {
 184   struct iu_params *p = (struct iu_params *) params;
 185   double a = p->a;
 186   gsl_function * f = p->f;
 187   double x = a + (1 - t) / t;
 188   double y = GSL_FN_EVAL (f, x);
 189   return (y / t) / t;
 190 }
 191
 192 /* Main integration function */
 193
 194 static int
 195 qags (const gsl_function * f,
 196       const double a, const double b,
 197       const double epsabs, const double epsrel,
 198       const size_t limit,
 199       gsl_integration_workspace * workspace,
 200       double *result, double *abserr,
 201       gsl_integration_rule * q)
 202 {
 203   double area, errsum;
 204   double res_ext, err_ext;
 205   double result0, abserr0, resabs0, resasc0;
 206   double tolerance;
 207
 208   double ertest = 0;
 209   double error_over_large_intervals = 0;
 210   double reseps = 0, abseps = 0, correc = 0;
 211   size_t ktmin = 0;
 212   int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
 213   int error_type = 0, error_type2 = 0;
 214
 215   size_t iteration = 0;
 216
 217   int positive_integrand = 0;
 218   int extrapolate = 0;
 219   int disallow_extrapolation = 0;
 220
 221   struct extrapolation_table table;
 222
 223   /* Initialize results */
 224
 225   initialise (workspace, a, b);
 226
 227   *result = 0;
 228   *abserr = 0;
 229
 230   if (limit > workspace->limit)
 231     {
 232       GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
 233     }
 234
 235   /* Test on accuracy */
 236
 237   if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
 238     {
 239       GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
 240                  GSL_EBADTOL);
 241     }
 242
 243   /* Perform the first integration */
 244
 245   q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
 246
 247   set_initial_result (workspace, result0, abserr0);
 248
 249   tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
 250
 251   if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance)
 252     {
 253       *result = result0;
 254       *abserr = abserr0;
 255
 256       GSL_ERROR ("cannot reach tolerance because of roundoff error"
 257                  "on first attempt", GSL_EROUND);
 258     }
 259   else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
 260     {
 261       *result = result0;
 262       *abserr = abserr0;
 263
 264       return GSL_SUCCESS;
 265     }
 266   else if (limit == 1)
 267     {
 268       *result = result0;
 269       *abserr = abserr0;
 270
 271       GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
 272     }
 273
 274   /* Initialization */
 275
 276   initialise_table (&table);
 277   append_table (&table, result0);
 278
 279   area = result0;
 280   errsum = abserr0;
 281
 282   res_ext = result0;
 283   err_ext = GSL_DBL_MAX;
 284
 285   positive_integrand = test_positivity (result0, resabs0);
 286
 287   iteration = 1;
 288
 289   do
 290     {
 291       size_t current_level;
 292       double a1, b1, a2, b2;
 293       double a_i, b_i, r_i, e_i;
 294       double area1 = 0, area2 = 0, area12 = 0;
 295       double error1 = 0, error2 = 0, error12 = 0;
 296       double resasc1, resasc2;
 297       double resabs1, resabs2;
 298       double last_e_i;
 299
 300       /* Bisect the subinterval with the largest error estimate */
 301
 302       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
 303
 304       current_level = workspace->level[workspace->i] + 1;
 305
 306       a1 = a_i;
 307       b1 = 0.5 * (a_i + b_i);
 308       a2 = b1;
 309       b2 = b_i;
 310
 311       iteration++;
 312
 313       q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
 314       q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
 315
 316       area12 = area1 + area2;
 317       error12 = error1 + error2;
 318       last_e_i = e_i;
 319
 320       /* Improve previous approximations to the integral and test for
 321          accuracy.
 322
 323          We write these expressions in the same way as the original
 324          QUADPACK code so that the rounding errors are the same, which
 325          makes testing easier. */
 326
 327       errsum = errsum + error12 - e_i;
 328       area = area + area12 - r_i;
 329
 330       tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
 331
 332       if (resasc1 != error1 && resasc2 != error2)
 333         {
 334           double delta = r_i - area12;
 335
 336           if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
 337             {
 338               if (!extrapolate)
 339                 {
 340                   roundoff_type1++;
 341                 }
 342               else
 343                 {
 344                   roundoff_type2++;
 345                 }
 346             }
 347           if (iteration > 10 && error12 > e_i)
 348             {
 349               roundoff_type3++;
 350             }
 351         }
 352
 353       /* Test for roundoff and eventually set error flag */
 354
 355       if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
 356         {
 357           error_type = 2;       /* round off error */
 358         }
 359
 360       if (roundoff_type2 >= 5)
 361         {
 362           error_type2 = 1;
 363         }
 364
 365       /* set error flag in the case of bad integrand behaviour at
 366          a point of the integration range */
 367
 368       if (subinterval_too_small (a1, a2, b2))
 369         {
 370           error_type = 4;
 371         }
 372
 373       /* append the newly-created intervals to the list */
 374
 375       update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
 376
 377       if (errsum <= tolerance)
 378         {
 379           goto compute_result;
 380         }
 381
 382       if (error_type)
 383         {
 384           break;
 385         }
 386
 387       if (iteration >= limit - 1)
 388         {
 389           error_type = 1;
 390           break;
 391         }
 392
 393       if (iteration == 2)       /* set up variables on first iteration */
 394         {
 395           error_over_large_intervals = errsum;
 396           ertest = tolerance;
 397           append_table (&table, area);
 398           continue;
 399         }
 400
 401       if (disallow_extrapolation)
 402         {
 403           continue;
 404         }
 405
 406       error_over_large_intervals += -last_e_i;
 407
 408       if (current_level < workspace->maximum_level)
 409         {
 410           error_over_large_intervals += error12;
 411         }
 412
 413       if (!extrapolate)
 414         {
 415           /* test whether the interval to be bisected next is the
 416              smallest interval. */
 417
 418           if (large_interval (workspace))
 419             continue;
 420
 421           extrapolate = 1;
 422           workspace->nrmax = 1;
 423         }
 424
 425       if (!error_type2 && error_over_large_intervals > ertest)
 426         {
 427           if (increase_nrmax (workspace))
 428             continue;
 429         }
 430
 431       /* Perform extrapolation */
 432
 433       append_table (&table, area);
 434
 435       qelg (&table, &reseps, &abseps);
 436
 437       ktmin++;
 438
 439       if (ktmin > 5 && err_ext < 0.001 * errsum)
 440         {
 441           error_type = 5;
 442         }
 443
 444       if (abseps < err_ext)
 445         {
 446           ktmin = 0;
 447           err_ext = abseps;
 448           res_ext = reseps;
 449           correc = error_over_large_intervals;
 450           ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
 451           if (err_ext <= ertest)
 452             break;
 453         }
 454
 455       /* Prepare bisection of the smallest interval. */
 456
 457       if (table.n == 1)
 458         {
 459           disallow_extrapolation = 1;
 460         }
 461
 462       if (error_type == 5)
 463         {
 464           break;
 465         }
 466
 467       /* work on interval with largest error */
 468
 469       reset_nrmax (workspace);
 470       extrapolate = 0;
 471       error_over_large_intervals = errsum;
 472
 473     }
 474   while (iteration < limit);
 475
 476   *result = res_ext;
 477   *abserr = err_ext;
 478
 479   if (err_ext == GSL_DBL_MAX)
 480     goto compute_result;
 481
 482   if (error_type || error_type2)
 483     {
 484       if (error_type2)
 485         {
 486           err_ext += correc;
 487         }
 488
 489       if (error_type == 0)
 490         error_type = 3;
 491
 492       if (res_ext != 0.0 && area != 0.0)
 493         {
 494           if (err_ext / fabs (res_ext) > errsum / fabs (area))
 495             goto compute_result;
 496         }
 497       else if (err_ext > errsum)
 498         {
 499           goto compute_result;
 500         }
 501       else if (area == 0.0)
 502         {
 503           goto return_error;
 504         }
 505     }
 506
 507   /*  Test on divergence. */
 508
 509   {
 510     double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
 511
 512     if (!positive_integrand && max_area < 0.01 * resabs0)
 513       goto return_error;
 514   }
 515
 516   {
 517     double ratio = res_ext / area;
 518
 519     if (ratio < 0.01 || ratio > 100.0 || errsum > fabs (area))
 520       error_type = 6;
 521   }
 522
 523   goto return_error;
 524
 525 compute_result:
 526
 527   *result = sum_results (workspace);
 528   *abserr = errsum;
 529
 530 return_error:
 531
 532   if (error_type > 2)
 533     error_type--;
 534
 535
 536
 537   if (error_type == 0)
 538     {
 539       return GSL_SUCCESS;
 540     }
 541   else if (error_type == 1)
 542     {
 543       GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
 544     }
 545   else if (error_type == 2)
 546     {
 547       GSL_ERROR ("cannot reach tolerance because of roundoff error",
 548                  GSL_EROUND);
 549     }
 550   else if (error_type == 3)
 551     {
 552       GSL_ERROR ("bad integrand behavior found in the integration interval",
 553                  GSL_ESING);
 554     }
 555   else if (error_type == 4)
 556     {
 557       GSL_ERROR ("roundoff error detected in the extrapolation table",
 558                  GSL_EROUND);
 559     }
 560   else if (error_type == 5)
 561     {
 562       GSL_ERROR ("integral is divergent, or slowly convergent",
 563                  GSL_EDIVERGE);
 564     }
 565   else
 566     {
 567       GSL_ERROR ("could not integrate function", GSL_EFAILED);
 568     }
 569
 570 }