ROOT   6.10/09 Reference Guide
RooAdaptiveGaussKronrodIntegrator1D.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * @(#)root/roofitcore:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16
17 /**
18 \file RooAdaptiveGaussKronrodIntegrator1D.cxx
19 \class RooAdaptiveGaussKronrodIntegrator1D
20 \ingroup Roofitcore
21
22 RooAdaptiveGaussKronrodIntegrator1D implements the Gauss-Kronrod integration algorithm.
23
24 An adaptive Gaussian quadrature method for numerical integration in
25 which error is estimation based on evaluation at special points
26 known as "Kronrod points." By suitably picking these points,
27 abscissas from previous iterations can be reused as part of the new
28 set of points, whereas usual Gaussian quadrature would require
29 recomputation of all abscissas at each iteration.
30
31 This class automatically handles (-inf,+inf) integrals by dividing
32 the integration in three regions (-inf,-1), (-1,1), (1,inf) and
33 calculating the 1st and 3rd term using a x -> 1/x coordinate
34 transformation
35
36 This class embeds the adaptive Gauss-Kronrod integrator from the
37 GNU Scientific Library version 1.5 and applies a chosen rule ( 10-,
38 21-, 31-, 41, 51- or 61-point). The integration domain is
39 subdivided and recursively integrated until the required precision
40 is reached.
41
42 For integrands with integrable singulaties the Wynn epsilon rule
43 can be selected to speed up the converges of these integrals
44 **/
45
46 #include "RooFit.h"
47
48 #include <assert.h>
49 #include <stdlib.h>
50 #include "TClass.h"
51 #include "Riostream.h"
53 #include "RooArgSet.h"
54 #include "RooRealVar.h"
55 #include "RooNumber.h"
56 #include "RooNumIntFactory.h"
57 #include "RooIntegratorBinding.h"
58 #include "TMath.h"
59 #include "RooMsgService.h"
60
61 using namespace std ;
62
63
65 ;
66
67
68 // --- From GSL_MATH.h -------------------------------------------
69 struct gsl_function_struct
70 {
71  double (* function) (double x, void * params);
72  void * params;
73 };
74 typedef struct gsl_function_struct gsl_function ;
75 #define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params)
76
77 //----From GSL_INTEGRATION.h ---------------------------------------
78 typedef struct
79  {
80  size_t limit;
81  size_t size;
82  size_t nrmax;
83  size_t i;
84  size_t maximum_level;
85  double *alist;
86  double *blist;
87  double *rlist;
88  double *elist;
89  size_t *order;
90  size_t *level;
91  }
92 gsl_integration_workspace;
93
94 gsl_integration_workspace *
95  gsl_integration_workspace_alloc (const size_t n);
96
97 void
98  gsl_integration_workspace_free (gsl_integration_workspace * w);
99
100 int gsl_integration_qag (const gsl_function * f,
101  double a, double b,
102  double epsabs, double epsrel, size_t limit,
103  int key,
104  gsl_integration_workspace * workspace,
105  double *result, double *abserr);
106
107 int
109  double a, double b,
110  double epsabs, double epsrel, size_t limit,
111  gsl_integration_workspace * workspace,
112  double * result, double * abserr) ;
113
114 int
116  double epsabs, double epsrel, size_t limit,
117  gsl_integration_workspace * workspace,
118  double *result, double *abserr) ;
119
120 int
122  double b,
123  double epsabs, double epsrel, size_t limit,
124  gsl_integration_workspace * workspace,
125  double *result, double *abserr) ;
126
127 int
129  double a,
130  double epsabs, double epsrel, size_t limit,
131  gsl_integration_workspace * workspace,
132  double *result, double *abserr) ;
133
134
135 //-------------------------------------------------------------------
136
137
138 ////////////////////////////////////////////////////////////////////////////////
139 /// Register this class with RooNumIntConfig as a possible choice of numeric
140 /// integrator for one-dimensional integrals over finite and infinite domains
141
143 {
144  RooRealVar maxSeg("maxSeg","maximum number of segments",100) ;
145  RooCategory method("method","Integration method for each segment") ;
146  method.defineType("WynnEpsilon",0) ;
147  method.defineType("15Points",1) ;
148  method.defineType("21Points",2) ;
149  method.defineType("31Points",3) ;
150  method.defineType("41Points",4) ;
151  method.defineType("51Points",5) ;
152  method.defineType("61Points",6) ;
153  method.setIndex(2) ;
155 }
156
157
158
159 ////////////////////////////////////////////////////////////////////////////////
160 /// coverity[UNINIT_CTOR]
161 /// Default constructor
162
164 {
165 }
166
167
168
169
170 ////////////////////////////////////////////////////////////////////////////////
171 /// Constructor taking a function binding and a configuration object
172
174  const RooNumIntConfig& config) :
176  _epsAbs(config.epsRel()),
177  _epsRel(config.epsAbs()),
178  _workspace(0)
179 {
180  // Use this form of the constructor to integrate over the function's default range.
181  const RooArgSet& confSet = config.getConfigSection(IsA()->GetName()) ;
182  _maxSeg = (Int_t) confSet.getRealValue("maxSeg",100) ;
183  _methodKey = confSet.getCatIndex("method",2) ;
184
186  _valid= initialize();
187 }
188
189
190
191 ////////////////////////////////////////////////////////////////////////////////
192 /// Constructor taking a function binding, an integration range and a configuration object
193
196  const RooNumIntConfig& config) :
198  _epsAbs(config.epsRel()),
199  _epsRel(config.epsAbs()),
200  _workspace(0),
201  _xmin(xmin),
202  _xmax(xmax)
203 {
204  // Use this form of the constructor to integrate over the function's default range.
205  const RooArgSet& confSet = config.getConfigSection(IsA()->GetName()) ;
206  _maxSeg = (Int_t) confSet.getRealValue("maxSeg",100) ;
207  _methodKey = confSet.getCatIndex("method",2) ;
208
210  _valid= initialize();
211 }
212
213
214
215 ////////////////////////////////////////////////////////////////////////////////
216 /// Virtual constructor
217
219 {
220  return new RooAdaptiveGaussKronrodIntegrator1D(function,config) ;
221 }
222
223
224
225 ////////////////////////////////////////////////////////////////////////////////
226 /// Initialize integrator allocate buffers and setup GSL workspace
227
229 {
230  // Allocate coordinate buffer size after number of function dimensions
231  _x = new Double_t[_function->getDimension()] ;
233
234  return checkLimits();
235 }
236
237
238
239 ////////////////////////////////////////////////////////////////////////////////
240 /// Destructor.
241
243 {
244  if (_workspace) {
245  gsl_integration_workspace_free ((gsl_integration_workspace*) _workspace) ;
246  }
247  if (_x) {
248  delete[] _x ;
249  }
250 }
251
252
253
254 ////////////////////////////////////////////////////////////////////////////////
255 /// Change our integration limits. Return kTRUE if the new limits are
256 /// ok, or otherwise kFALSE. Always returns kFALSE and does nothing
257 /// if this object was constructed to always use our integrand's limits.
258
260 {
261  if(_useIntegrandLimits) {
262  coutE(Integration) << "RooAdaptiveGaussKronrodIntegrator1D::setLimits: cannot override integrand's limits" << endl;
263  return kFALSE;
264  }
265
266  _xmin= *xmin;
267  _xmax= *xmax;
268  return checkLimits();
269 }
270
271
272
273 ////////////////////////////////////////////////////////////////////////////////
274 /// Check that our integration range is finite and otherwise return kFALSE.
275 /// Update the limits from the integrand if requested.
276
278 {
279  if(_useIntegrandLimits) {
280  assert(0 != integrand() && integrand()->isValid());
281  _xmin= integrand()->getMinLimit(0);
282  _xmax= integrand()->getMaxLimit(0);
283  }
284
285  // Determine domain type
288
289  if (!infLo && !infHi) {
290  _domainType = Closed ;
291  } else if (infLo && infHi) {
292  _domainType = Open ;
293  } else if (infLo && !infHi) {
294  _domainType = OpenLo ;
295  } else {
296  _domainType = OpenHi ;
297  }
298
299
300  return kTRUE ;
301 }
302
303
304
305 ////////////////////////////////////////////////////////////////////////////////
306 /// Glue function interacing to GSL code
307
309 {
311  return instance->integrand(instance->xvec(x)) ;
312 }
313
314
315
316 ////////////////////////////////////////////////////////////////////////////////
317 /// Calculate and return integral at at given parameter values
318
320 {
321  assert(isValid());
322
323  // Copy yvec to xvec if provided
324  if (yvec) {
325  UInt_t i ; for (i=0 ; i<_function->getDimension()-1 ; i++) {
326  _x[i+1] = yvec[i] ;
327  }
328  }
329
330  // Setup glue function
331  gsl_function F;
333  F.params = this ;
334
335  // Return values
336  double result, error;
337
338  // Call GSL implementation of integeator
339  switch(_domainType) {
340  case Closed:
341  if (_methodKey==0) {
342  gsl_integration_qags (&F, _xmin, _xmax, _epsAbs, _epsRel, _maxSeg, (gsl_integration_workspace*)_workspace,&result, &error);
343  } else {
344  gsl_integration_qag (&F, _xmin, _xmax, _epsAbs, _epsRel, _maxSeg, _methodKey, (gsl_integration_workspace*)_workspace,&result, &error);
345  }
346  break ;
347  case OpenLo:
348  gsl_integration_qagil (&F, _xmax, _epsAbs, _epsRel, _maxSeg, (gsl_integration_workspace*)_workspace,&result, &error);
349  break ;
350  case OpenHi:
351  gsl_integration_qagiu (&F, _xmin, _epsAbs, _epsRel, _maxSeg, (gsl_integration_workspace*)_workspace,&result, &error);
352  break ;
353  case Open:
354  gsl_integration_qagi (&F, _epsAbs, _epsRel, _maxSeg, (gsl_integration_workspace*)_workspace,&result, &error);
355  break ;
356  }
357
358  return result;
359 }
360
361
362 // ----------------------------------------------------------------------------
363 // ---------- Code below imported from GSL version 1.5 ------------------------
364 // ----------------------------------------------------------------------------
365
366 /*
367  *
368  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
369  *
370  * This program is free software; you can redistribute it and/or modify
371  * it under the terms of the GNU General Public License as published by
372  * the Free Software Foundation; either version 2 of the License, or (at
373  * your option) any later version.
374  *
375  * This program is distributed in the hope that it will be useful, but
376  * WITHOUT ANY WARRANTY; without even the implied warranty of
377  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
378  * General Public License for more details.
379  *
380  * You should have received a copy of the GNU General Public License
381  * along with this program; if not, write to the Free Software
382  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
383  */
384
385 #define GSL_SUCCESS 0
386 #define GSL_EDOM 1 /* input domain error, e.g sqrt(-1) */
387 #define GSL_ENOMEM 8 /* malloc failed */
388 #define GSL_EBADTOL 13 /* user specified an invalid tolerance */
389 #define GSL_ETOL 14 /* failed to reach the specified tolerance */
390 #define GSL_ERROR(a,b) oocoutE((TObject*)0,Integration) << "RooAdaptiveGaussKronrodIntegrator1D::integral() ERROR: " << a << endl ; return b ;
391 #define GSL_DBL_MIN 2.2250738585072014e-308
392 #define GSL_DBL_MAX 1.7976931348623157e+308
393 #define GSL_DBL_EPSILON 2.2204460492503131e-16
394
395 #define GSL_EINVAL 2
396 #define GSL_EMAXITER 3
397 #define GSL_ESING 4
398 #define GSL_EFAILED 5
399 #define GSL_EDIVERGE 6
400 #define GSL_EROUND 7
401
402 #define GSL_ERROR_VAL(reason, gsl_errno, value) return value ;
403
404 #define GSL_MAX(a,b) ((a) > (b) ? (a) : (b))
405 extern inline double
406 GSL_MAX_DBL (double a, double b)
407 {
408  return GSL_MAX (a, b);
409 }
410
411 double gsl_coerce_double (const double x);
412
413 double
414 gsl_coerce_double (const double x)
415 {
416  volatile double y;
417  y = x;
418  return y;
419 }
420 #define GSL_COERCE_DBL(x) (gsl_coerce_double(x))
421
422 // INCLUDED BELOW #include <gsl/gsl_integration.h>
423
424
425 /* Workspace for adaptive integrators */
426 // WVE MOVED TO HEAD OF FILE
427
428
429 /* Definition of an integration rule */
430
431 typedef void gsl_integration_rule (const gsl_function * f,
432  double a, double b,
433  double *result, double *abserr,
434  double *defabs, double *resabs);
435
436 void gsl_integration_qk15 (const gsl_function * f, double a, double b,
437  double *result, double *abserr,
438  double *resabs, double *resasc);
439
440 void gsl_integration_qk21 (const gsl_function * f, double a, double b,
441  double *result, double *abserr,
442  double *resabs, double *resasc);
443
444 void gsl_integration_qk31 (const gsl_function * f, double a, double b,
445  double *result, double *abserr,
446  double *resabs, double *resasc);
447
448 void gsl_integration_qk41 (const gsl_function * f, double a, double b,
449  double *result, double *abserr,
450  double *resabs, double *resasc);
451
452 void gsl_integration_qk51 (const gsl_function * f, double a, double b,
453  double *result, double *abserr,
454  double *resabs, double *resasc);
455
456 void gsl_integration_qk61 (const gsl_function * f, double a, double b,
457  double *result, double *abserr,
458  double *resabs, double *resasc);
459
460 void gsl_integration_qcheb (gsl_function * f, double a, double b,
461  double *cheb12, double *cheb24);
462
463 /* The low-level integration rules in QUADPACK are identified by small
464  integers (1-6). We'll use symbolic constants to refer to them. */
465
466 enum
467  {
468  GSL_INTEG_GAUSS15 = 1, /* 15 point Gauss-Kronrod rule */
469  GSL_INTEG_GAUSS21 = 2, /* 21 point Gauss-Kronrod rule */
470  GSL_INTEG_GAUSS31 = 3, /* 31 point Gauss-Kronrod rule */
471  GSL_INTEG_GAUSS41 = 4, /* 41 point Gauss-Kronrod rule */
472  GSL_INTEG_GAUSS51 = 5, /* 51 point Gauss-Kronrod rule */
473  GSL_INTEG_GAUSS61 = 6 /* 61 point Gauss-Kronrod rule */
474  };
475
476
477 void
478 gsl_integration_qk (const int n, const double xgk[],
479  const double wg[], const double wgk[],
480  double fv1[], double fv2[],
481  const gsl_function *f, double a, double b,
482  double * result, double * abserr,
483  double * resabs, double * resasc);
484
485
486 int gsl_integration_qag (const gsl_function * f,
487  double a, double b,
488  double epsabs, double epsrel, size_t limit,
489  int key,
490  gsl_integration_workspace * workspace,
491  double *result, double *abserr);
492
493
494 // INCLUDED BELOW #include "initialise.c"
495 static inline
496 void initialise (gsl_integration_workspace * workspace, double a, double b);
497
498 static inline
499 void initialise (gsl_integration_workspace * workspace, double a, double b)
500 {
501  workspace->size = 0;
502  workspace->nrmax = 0;
503  workspace->i = 0;
504  workspace->alist[0] = a;
505  workspace->blist[0] = b;
506  workspace->rlist[0] = 0.0;
507  workspace->elist[0] = 0.0;
508  workspace->order[0] = 0;
509  workspace->level[0] = 0;
510
511  workspace->maximum_level = 0;
512 }
513
514 // INCLUDED BELOW #include "set_initial.c"
515 static inline
516 void set_initial_result (gsl_integration_workspace * workspace,
517  double result, double error);
518
519 static inline
520 void set_initial_result (gsl_integration_workspace * workspace,
521  double result, double error)
522 {
523  workspace->size = 1;
524  workspace->rlist[0] = result;
525  workspace->elist[0] = error;
526 }
527
528 // INCLUDED BELOW #include "qpsrt.c"
529 static inline void
530 qpsrt (gsl_integration_workspace * workspace);
531
532 static inline
533 void qpsrt (gsl_integration_workspace * workspace)
534 {
535  const size_t last = workspace->size - 1;
536  const size_t limit = workspace->limit;
537
538  double * elist = workspace->elist;
539  size_t * order = workspace->order;
540
541  double errmax ;
542  double errmin ;
543  int i, k, top;
544
545  size_t i_nrmax = workspace->nrmax;
546  size_t i_maxerr = order[i_nrmax] ;
547
548  /* Check whether the list contains more than two error estimates */
549
550  if (last < 2)
551  {
552  order[0] = 0 ;
553  order[1] = 1 ;
554  workspace->i = i_maxerr ;
555  return ;
556  }
557
558  errmax = elist[i_maxerr] ;
559
560  /* This part of the routine is only executed if, due to a difficult
561  integrand, subdivision increased the error estimate. In the normal
562  case the insert procedure should start after the nrmax-th largest
563  error estimate. */
564
565  while (i_nrmax > 0 && errmax > elist[order[i_nrmax - 1]])
566  {
567  order[i_nrmax] = order[i_nrmax - 1] ;
568  i_nrmax-- ;
569  }
570
571  /* Compute the number of elements in the list to be maintained in
572  descending order. This number depends on the number of
573  subdivisions still allowed. */
574
575  if(last < (limit/2 + 2))
576  {
577  top = last ;
578  }
579  else
580  {
581  top = limit - last + 1;
582  }
583
584  /* Insert errmax by traversing the list top-down, starting
585  comparison from the element elist(order(i_nrmax+1)). */
586
587  i = i_nrmax + 1 ;
588
589  /* The order of the tests in the following line is important to
590  prevent a segmentation fault */
591
592  while (i < top && errmax < elist[order[i]])
593  {
594  order[i-1] = order[i] ;
595  i++ ;
596  }
597
598  order[i-1] = i_maxerr ;
599
600  /* Insert errmin by traversing the list bottom-up */
601
602  errmin = elist[last] ;
603
604  k = top - 1 ;
605
606  while (k > i - 2 && errmin >= elist[order[k]])
607  {
608  order[k+1] = order[k] ;
609  k-- ;
610  }
611
612  order[k+1] = last ;
613
614  /* Set i_max and e_max */
615
616  i_maxerr = order[i_nrmax] ;
617
618  workspace->i = i_maxerr ;
619  workspace->nrmax = i_nrmax ;
620 }
621
622
623
624 // INCLUDED BELOW #include "util.c"
625 static inline
626 void update (gsl_integration_workspace * workspace,
627  double a1, double b1, double area1, double error1,
628  double a2, double b2, double area2, double error2);
629
630 static inline void
631 retrieve (const gsl_integration_workspace * workspace,
632  double * a, double * b, double * r, double * e);
633
634
635
636 static inline
637 void update (gsl_integration_workspace * workspace,
638  double a1, double b1, double area1, double error1,
639  double a2, double b2, double area2, double error2)
640 {
641  double * alist = workspace->alist ;
642  double * blist = workspace->blist ;
643  double * rlist = workspace->rlist ;
644  double * elist = workspace->elist ;
645  size_t * level = workspace->level ;
646
647  const size_t i_max = workspace->i ;
648  const size_t i_new = workspace->size ;
649
650  const size_t new_level = workspace->level[i_max] + 1;
651
652  /* append the newly-created intervals to the list */
653
654  if (error2 > error1)
655  {
656  alist[i_max] = a2; /* blist[maxerr] is already == b2 */
657  rlist[i_max] = area2;
658  elist[i_max] = error2;
659  level[i_max] = new_level;
660
661  alist[i_new] = a1;
662  blist[i_new] = b1;
663  rlist[i_new] = area1;
664  elist[i_new] = error1;
665  level[i_new] = new_level;
666  }
667  else
668  {
669  blist[i_max] = b1; /* alist[maxerr] is already == a1 */
670  rlist[i_max] = area1;
671  elist[i_max] = error1;
672  level[i_max] = new_level;
673
674  alist[i_new] = a2;
675  blist[i_new] = b2;
676  rlist[i_new] = area2;
677  elist[i_new] = error2;
678  level[i_new] = new_level;
679  }
680
681  workspace->size++;
682
683  if (new_level > workspace->maximum_level)
684  {
685  workspace->maximum_level = new_level;
686  }
687
688  qpsrt (workspace) ;
689 }
690
691 static inline void
692 retrieve (const gsl_integration_workspace * workspace,
693  double * a, double * b, double * r, double * e)
694 {
695  const size_t i = workspace->i;
696  double * alist = workspace->alist;
697  double * blist = workspace->blist;
698  double * rlist = workspace->rlist;
699  double * elist = workspace->elist;
700
701  *a = alist[i] ;
702  *b = blist[i] ;
703  *r = rlist[i] ;
704  *e = elist[i] ;
705 }
706
707 static inline double
708 sum_results (const gsl_integration_workspace * workspace);
709
710 static inline double
711 sum_results (const gsl_integration_workspace * workspace)
712 {
713  const double * const rlist = workspace->rlist ;
714  const size_t n = workspace->size;
715
716  size_t k;
717  double result_sum = 0;
718
719  for (k = 0; k < n; k++)
720  {
721  result_sum += rlist[k];
722  }
723
724  return result_sum;
725 }
726
727 static inline int
728 subinterval_too_small (double a1, double a2, double b2);
729
730 static inline int
731 subinterval_too_small (double a1, double a2, double b2)
732 {
733  const double e = GSL_DBL_EPSILON;
734  const double u = GSL_DBL_MIN;
735
736  double tmp = (1 + 100 * e) * (fabs (a2) + 1000 * u);
737
738  int status = fabs (a1) <= tmp && fabs (b2) <= tmp;
739
740  return status;
741 }
742
743
744 static int
745 qag (const gsl_function *f,
746  const double a, const double b,
747  const double epsabs, const double epsrel,
748  const size_t limit,
749  gsl_integration_workspace * workspace,
750  double * result, double * abserr,
752
753 int
755  double a, double b,
756  double epsabs, double epsrel, size_t limit,
757  int key,
758  gsl_integration_workspace * workspace,
759  double * result, double * abserr)
760 {
761  int status ;
762  gsl_integration_rule * integration_rule = gsl_integration_qk15 ;
763
764  if (key < GSL_INTEG_GAUSS15)
765  {
766  key = GSL_INTEG_GAUSS15 ;
767  }
768  else if (key > GSL_INTEG_GAUSS61)
769  {
770  key = GSL_INTEG_GAUSS61 ;
771  }
772
773  switch (key)
774  {
775  case GSL_INTEG_GAUSS15:
776  integration_rule = gsl_integration_qk15 ;
777  break ;
778  case GSL_INTEG_GAUSS21:
779  integration_rule = gsl_integration_qk21 ;
780  break ;
781  case GSL_INTEG_GAUSS31:
782  integration_rule = gsl_integration_qk31 ;
783  break ;
784  case GSL_INTEG_GAUSS41:
785  integration_rule = gsl_integration_qk41 ;
786  break ;
787  case GSL_INTEG_GAUSS51:
788  integration_rule = gsl_integration_qk51 ;
789  break ;
790  case GSL_INTEG_GAUSS61:
791  integration_rule = gsl_integration_qk61 ;
792  break ;
793  }
794
795  status = qag (f, a, b, epsabs, epsrel, limit,
796  workspace,
797  result, abserr,
798  integration_rule) ;
799
800  return status ;
801 }
802
803 static int
804 qag (const gsl_function * f,
805  const double a, const double b,
806  const double epsabs, const double epsrel,
807  const size_t limit,
808  gsl_integration_workspace * workspace,
809  double *result, double *abserr,
811 {
812  double area, errsum;
813  double result0, abserr0, resabs0, resasc0;
814  double tolerance;
815  size_t iteration = 0;
816  int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
817
818  double round_off;
819
820  /* Initialize results */
821
822  initialise (workspace, a, b);
823
824  *result = 0;
825  *abserr = 0;
826
827  if (limit > workspace->limit)
828  {
829  GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
830  }
831
832  if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
833  {
834  GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
835  GSL_EBADTOL);
836  }
837
838  /* perform the first integration */
839
840  q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
841
842  set_initial_result (workspace, result0, abserr0);
843
844  /* Test on accuracy */
845
846  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
847
848  /* need IEEE rounding here to match original quadpack behavior */
849
850  round_off = GSL_COERCE_DBL (50 * GSL_DBL_EPSILON * resabs0);
851
852  if (abserr0 <= round_off && abserr0 > tolerance)
853  {
854  *result = result0;
855  *abserr = abserr0;
856
857  GSL_ERROR ("cannot reach tolerance because of roundoff error "
858  "on first attempt", GSL_EROUND);
859  }
860  else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
861  {
862  *result = result0;
863  *abserr = abserr0;
864
865  return GSL_SUCCESS;
866  }
867  else if (limit == 1)
868  {
869  *result = result0;
870  *abserr = abserr0;
871
872  GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
873  }
874
875  area = result0;
876  errsum = abserr0;
877
878  iteration = 1;
879
880  do
881  {
882  double a1, b1, a2, b2;
883  double a_i, b_i, r_i, e_i;
884  double area1 = 0, area2 = 0, area12 = 0;
885  double error1 = 0, error2 = 0, error12 = 0;
886  double resasc1, resasc2;
887  double resabs1, resabs2;
888
889  /* Bisect the subinterval with the largest error estimate */
890
891  retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
892
893  a1 = a_i;
894  b1 = 0.5 * (a_i + b_i);
895  a2 = b1;
896  b2 = b_i;
897
898  q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
899  q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
900
901  area12 = area1 + area2;
902  error12 = error1 + error2;
903
904  errsum += (error12 - e_i);
905  area += area12 - r_i;
906
907  if (resasc1 != error1 && resasc2 != error2)
908  {
909  double delta = r_i - area12;
910
911  if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
912  {
913  roundoff_type1++;
914  }
915  if (iteration >= 10 && error12 > e_i)
916  {
917  roundoff_type2++;
918  }
919  }
920
921  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
922
923  if (errsum > tolerance)
924  {
925  if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
926  {
927  error_type = 2; /* round off error */
928  }
929
930  /* set error flag in the case of bad integrand behaviour at
931  a point of the integration range */
932
933  if (subinterval_too_small (a1, a2, b2))
934  {
935  error_type = 3;
936  }
937  }
938
939  update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
940
941  retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
942
943  iteration++;
944
945  }
946  while (iteration < limit && !error_type && errsum > tolerance);
947
948  *result = sum_results (workspace);
949  *abserr = errsum;
950
951  if (errsum <= tolerance)
952  {
953  return GSL_SUCCESS;
954  }
955  else if (error_type == 2)
956  {
957  GSL_ERROR ("roundoff error prevents tolerance from being achieved",
958  GSL_EROUND);
959  }
960  else if (error_type == 3)
961  {
962  GSL_ERROR ("bad integrand behavior found in the integration interval",
963  GSL_ESING);
964  }
965  else if (iteration == limit)
966  {
967  GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
968  }
969
970  GSL_ERROR ("could not integrate function", GSL_EFAILED);
971 }
972
973
974 // INCLUDED BELOW #include "err.c"
975 static double rescale_error (double err, const double result_abs, const double result_asc) ;
976
977 static double
978 rescale_error (double err, const double result_abs, const double result_asc)
979 {
980  err = fabs(err) ;
981
982  if (result_asc != 0 && err != 0)
983  {
984  double scale = TMath::Power((200 * err / result_asc), 1.5) ;
985
986  if (scale < 1)
987  {
988  err = result_asc * scale ;
989  }
990  else
991  {
992  err = result_asc ;
993  }
994  }
995  if (result_abs > GSL_DBL_MIN / (50 * GSL_DBL_EPSILON))
996  {
997  double min_err = 50 * GSL_DBL_EPSILON * result_abs ;
998
999  if (min_err > err)
1000  {
1001  err = min_err ;
1002  }
1003  }
1004
1005  return err ;
1006 }
1007
1008
1009 void
1010 gsl_integration_qk (const int n,
1011  const double xgk[], const double wg[], const double wgk[],
1012  double fv1[], double fv2[],
1013  const gsl_function * f, double a, double b,
1014  double *result, double *abserr,
1015  double *resabs, double *resasc)
1016 {
1017
1018  const double center = 0.5 * (a + b);
1019  const double half_length = 0.5 * (b - a);
1020  const double abs_half_length = fabs (half_length);
1021  const double f_center = GSL_FN_EVAL (f, center);
1022
1023  double result_gauss = 0;
1024  double result_kronrod = f_center * wgk[n - 1];
1025
1026  double result_abs = fabs (result_kronrod);
1027  double result_asc = 0;
1028  double mean = 0, err = 0;
1029
1030  int j;
1031
1032  if (n % 2 == 0)
1033  {
1034  result_gauss = f_center * wg[n / 2 - 1];
1035  }
1036
1037  for (j = 0; j < (n - 1) / 2; j++)
1038  {
1039  const int jtw = j * 2 + 1; /* j=1,2,3 jtw=2,4,6 */
1040  const double abscissa = half_length * xgk[jtw];
1041  const double fval1 = GSL_FN_EVAL (f, center - abscissa);
1042  const double fval2 = GSL_FN_EVAL (f, center + abscissa);
1043  const double fsum = fval1 + fval2;
1044  fv1[jtw] = fval1;
1045  fv2[jtw] = fval2;
1046  result_gauss += wg[j] * fsum;
1047  result_kronrod += wgk[jtw] * fsum;
1048  result_abs += wgk[jtw] * (fabs (fval1) + fabs (fval2));
1049  }
1050
1051  for (j = 0; j < n / 2; j++)
1052  {
1053  int jtwm1 = j * 2;
1054  const double abscissa = half_length * xgk[jtwm1];
1055  const double fval1 = GSL_FN_EVAL (f, center - abscissa);
1056  const double fval2 = GSL_FN_EVAL (f, center + abscissa);
1057  fv1[jtwm1] = fval1;
1058  fv2[jtwm1] = fval2;
1059  result_kronrod += wgk[jtwm1] * (fval1 + fval2);
1060  result_abs += wgk[jtwm1] * (fabs (fval1) + fabs (fval2));
1061  };
1062
1063  mean = result_kronrod * 0.5;
1064
1065  result_asc = wgk[n - 1] * fabs (f_center - mean);
1066
1067  for (j = 0; j < n - 1; j++)
1068  {
1069  result_asc += wgk[j] * (fabs (fv1[j] - mean) + fabs (fv2[j] - mean));
1070  }
1071
1072  /* scale by the width of the integration region */
1073
1074  err = (result_kronrod - result_gauss) * half_length;
1075
1076  result_kronrod *= half_length;
1077  result_abs *= abs_half_length;
1078  result_asc *= abs_half_length;
1079
1080  *result = result_kronrod;
1081  *resabs = result_abs;
1082  *resasc = result_asc;
1083  *abserr = rescale_error (err, result_abs, result_asc);
1084
1085 }
1086
1087 /* Gauss quadrature weights and kronrod quadrature abscissae and
1088  weights as evaluated with 80 decimal digit arithmetic by
1089  L. W. Fullerton, Bell Labs, Nov. 1981. */
1090
1091 static const double xgkA[8] = /* abscissae of the 15-point kronrod rule */
1092 {
1093  0.991455371120812639206854697526329,
1094  0.949107912342758524526189684047851,
1095  0.864864423359769072789712788640926,
1096  0.741531185599394439863864773280788,
1097  0.586087235467691130294144838258730,
1098  0.405845151377397166906606412076961,
1099  0.207784955007898467600689403773245,
1100  0.000000000000000000000000000000000
1101 };
1102
1103 /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule.
1104  xgk[0], xgk[2], ... abscissae to optimally extend the 7-point gauss rule */
1105
1106 static const double wgA[4] = /* weights of the 7-point gauss rule */
1107 {
1108  0.129484966168869693270611432679082,
1109  0.279705391489276667901467771423780,
1110  0.381830050505118944950369775488975,
1111  0.417959183673469387755102040816327
1112 };
1113
1114 static const double wgkA[8] = /* weights of the 15-point kronrod rule */
1115 {
1116  0.022935322010529224963732008058970,
1117  0.063092092629978553290700663189204,
1118  0.104790010322250183839876322541518,
1119  0.140653259715525918745189590510238,
1120  0.169004726639267902826583426598550,
1121  0.190350578064785409913256402421014,
1122  0.204432940075298892414161999234649,
1123  0.209482141084727828012999174891714
1124 };
1125
1126 void
1127 gsl_integration_qk15 (const gsl_function * f, double a, double b,
1128  double *result, double *abserr,
1129  double *resabs, double *resasc)
1130 {
1131  double fv1[8], fv2[8];
1132  // coverity[UNINIT_CTOR]
1133  gsl_integration_qk (8, xgkA, wgA, wgkA, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
1134 }
1135
1136 /* Gauss quadrature weights and kronrod quadrature abscissae and
1137  weights as evaluated with 80 decimal digit arithmetic by
1138  L. W. Fullerton, Bell Labs, Nov. 1981. */
1139
1140 static const double xgkB[11] = /* abscissae of the 21-point kronrod rule */
1141 {
1142  0.995657163025808080735527280689003,
1143  0.973906528517171720077964012084452,
1144  0.930157491355708226001207180059508,
1145  0.865063366688984510732096688423493,
1146  0.780817726586416897063717578345042,
1147  0.679409568299024406234327365114874,
1148  0.562757134668604683339000099272694,
1149  0.433395394129247190799265943165784,
1150  0.294392862701460198131126603103866,
1151  0.148874338981631210884826001129720,
1152  0.000000000000000000000000000000000
1153 };
1154
1155 /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule.
1156  xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */
1157
1158 static const double wgB[5] = /* weights of the 10-point gauss rule */
1159 {
1160  0.066671344308688137593568809893332,
1161  0.149451349150580593145776339657697,
1162  0.219086362515982043995534934228163,
1163  0.269266719309996355091226921569469,
1164  0.295524224714752870173892994651338
1165 };
1166
1167 static const double wgkB[11] = /* weights of the 21-point kronrod rule */
1168 {
1169  0.011694638867371874278064396062192,
1170  0.032558162307964727478818972459390,
1171  0.054755896574351996031381300244580,
1172  0.075039674810919952767043140916190,
1173  0.093125454583697605535065465083366,
1174  0.109387158802297641899210590325805,
1175  0.123491976262065851077958109831074,
1176  0.134709217311473325928054001771707,
1177  0.142775938577060080797094273138717,
1178  0.147739104901338491374841515972068,
1179  0.149445554002916905664936468389821
1180 };
1181
1182
1183 void
1184 gsl_integration_qk21 (const gsl_function * f, double a, double b,
1185  double *result, double *abserr,
1186  double *resabs, double *resasc)
1187 {
1188  double fv1[11], fv2[11];
1189  // coverity[UNINIT_CTOR]
1190  gsl_integration_qk (11, xgkB, wgB, wgkB, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
1191 }
1192
1193 /* Gauss quadrature weights and kronrod quadrature abscissae and
1194  weights as evaluated with 80 decimal digit arithmetic by
1195  L. W. Fullerton, Bell Labs, Nov. 1981. */
1196
1197 static const double xgkC[16] = /* abscissae of the 31-point kronrod rule */
1198 {
1199  0.998002298693397060285172840152271,
1200  0.987992518020485428489565718586613,
1201  0.967739075679139134257347978784337,
1202  0.937273392400705904307758947710209,
1203  0.897264532344081900882509656454496,
1204  0.848206583410427216200648320774217,
1205  0.790418501442465932967649294817947,
1206  0.724417731360170047416186054613938,
1207  0.650996741297416970533735895313275,
1208  0.570972172608538847537226737253911,
1209  0.485081863640239680693655740232351,
1210  0.394151347077563369897207370981045,
1211  0.299180007153168812166780024266389,
1212  0.201194093997434522300628303394596,
1213  0.101142066918717499027074231447392,
1214  0.000000000000000000000000000000000
1215 };
1216
1217 /* xgk[1], xgk[3], ... abscissae of the 15-point gauss rule.
1218  xgk[0], xgk[2], ... abscissae to optimally extend the 15-point gauss rule */
1219
1220 static const double wgC[8] = /* weights of the 15-point gauss rule */
1221 {
1222  0.030753241996117268354628393577204,
1223  0.070366047488108124709267416450667,
1224  0.107159220467171935011869546685869,
1225  0.139570677926154314447804794511028,
1226  0.166269205816993933553200860481209,
1227  0.186161000015562211026800561866423,
1228  0.198431485327111576456118326443839,
1229  0.202578241925561272880620199967519
1230 };
1231
1232 static const double wgkC[16] = /* weights of the 31-point kronrod rule */
1233 {
1234  0.005377479872923348987792051430128,
1235  0.015007947329316122538374763075807,
1236  0.025460847326715320186874001019653,
1237  0.035346360791375846222037948478360,
1238  0.044589751324764876608227299373280,
1239  0.053481524690928087265343147239430,
1240  0.062009567800670640285139230960803,
1241  0.069854121318728258709520077099147,
1242  0.076849680757720378894432777482659,
1243  0.083080502823133021038289247286104,
1244  0.088564443056211770647275443693774,
1245  0.093126598170825321225486872747346,
1246  0.096642726983623678505179907627589,
1247  0.099173598721791959332393173484603,
1248  0.100769845523875595044946662617570,
1249  0.101330007014791549017374792767493
1250 };
1251
1252 void
1253 gsl_integration_qk31 (const gsl_function * f, double a, double b,
1254  double *result, double *abserr,
1255  double *resabs, double *resasc)
1256 {
1257  double fv1[16], fv2[16];
1258  // coverity[UNINIT_CTOR]
1259  gsl_integration_qk (16, xgkC, wgC, wgkC, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
1260 }
1261
1262 /* Gauss quadrature weights and kronrod quadrature abscissae and
1263  weights as evaluated with 80 decimal digit arithmetic by
1264  L. W. Fullerton, Bell Labs, Nov. 1981. */
1265
1266 static const double xgkD[21] = /* abscissae of the 41-point kronrod rule */
1267 {
1268  0.998859031588277663838315576545863,
1269  0.993128599185094924786122388471320,
1270  0.981507877450250259193342994720217,
1271  0.963971927277913791267666131197277,
1272  0.940822633831754753519982722212443,
1273  0.912234428251325905867752441203298,
1274  0.878276811252281976077442995113078,
1275  0.839116971822218823394529061701521,
1276  0.795041428837551198350638833272788,
1277  0.746331906460150792614305070355642,
1278  0.693237656334751384805490711845932,
1279  0.636053680726515025452836696226286,
1280  0.575140446819710315342946036586425,
1281  0.510867001950827098004364050955251,
1282  0.443593175238725103199992213492640,
1283  0.373706088715419560672548177024927,
1284  0.301627868114913004320555356858592,
1285  0.227785851141645078080496195368575,
1286  0.152605465240922675505220241022678,
1287  0.076526521133497333754640409398838,
1288  0.000000000000000000000000000000000
1289 };
1290
1291 /* xgk[1], xgk[3], ... abscissae of the 20-point gauss rule.
1292  xgk[0], xgk[2], ... abscissae to optimally extend the 20-point gauss rule */
1293
1294 static const double wgD[11] = /* weights of the 20-point gauss rule */
1295 {
1296  0.017614007139152118311861962351853,
1297  0.040601429800386941331039952274932,
1298  0.062672048334109063569506535187042,
1299  0.083276741576704748724758143222046,
1300  0.101930119817240435036750135480350,
1301  0.118194531961518417312377377711382,
1302  0.131688638449176626898494499748163,
1303  0.142096109318382051329298325067165,
1304  0.149172986472603746787828737001969,
1305  0.152753387130725850698084331955098
1306 };
1307
1308 static const double wgkD[21] = /* weights of the 41-point kronrod rule */
1309 {
1310  0.003073583718520531501218293246031,
1311  0.008600269855642942198661787950102,
1312  0.014626169256971252983787960308868,
1313  0.020388373461266523598010231432755,
1314  0.025882133604951158834505067096153,
1315  0.031287306777032798958543119323801,
1316  0.036600169758200798030557240707211,
1317  0.041668873327973686263788305936895,
1318  0.046434821867497674720231880926108,
1319  0.050944573923728691932707670050345,
1320  0.055195105348285994744832372419777,
1321  0.059111400880639572374967220648594,
1322  0.062653237554781168025870122174255,
1323  0.065834597133618422111563556969398,
1324  0.068648672928521619345623411885368,
1325  0.071054423553444068305790361723210,
1326  0.073030690332786667495189417658913,
1327  0.074582875400499188986581418362488,
1328  0.075704497684556674659542775376617,
1329  0.076377867672080736705502835038061,
1330  0.076600711917999656445049901530102
1331 };
1332
1333 void
1334 gsl_integration_qk41 (const gsl_function * f, double a, double b,
1335  double *result, double *abserr,
1336  double *resabs, double *resasc)
1337 {
1338  double fv1[21], fv2[21];
1339  // coverity[UNINIT]
1340  gsl_integration_qk (21, xgkD, wgD, wgkD, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
1341 }
1342
1343 /* Gauss quadrature weights and kronrod quadrature abscissae and
1344  weights as evaluated with 80 decimal digit arithmetic by
1345  L. W. Fullerton, Bell Labs, Nov. 1981. */
1346
1347 static const double xgkE[26] = /* abscissae of the 51-point kronrod rule */
1348 {
1349  0.999262104992609834193457486540341,
1350  0.995556969790498097908784946893902,
1351  0.988035794534077247637331014577406,
1352  0.976663921459517511498315386479594,
1353  0.961614986425842512418130033660167,
1354  0.942974571228974339414011169658471,
1355  0.920747115281701561746346084546331,
1356  0.894991997878275368851042006782805,
1357  0.865847065293275595448996969588340,
1358  0.833442628760834001421021108693570,
1359  0.797873797998500059410410904994307,
1360  0.759259263037357630577282865204361,
1361  0.717766406813084388186654079773298,
1362  0.673566368473468364485120633247622,
1363  0.626810099010317412788122681624518,
1364  0.577662930241222967723689841612654,
1365  0.526325284334719182599623778158010,
1366  0.473002731445714960522182115009192,
1367  0.417885382193037748851814394594572,
1368  0.361172305809387837735821730127641,
1369  0.303089538931107830167478909980339,
1370  0.243866883720988432045190362797452,
1371  0.183718939421048892015969888759528,
1372  0.122864692610710396387359818808037,
1373  0.061544483005685078886546392366797,
1374  0.000000000000000000000000000000000
1375 };
1376
1377 /* xgk[1], xgk[3], ... abscissae of the 25-point gauss rule.
1378  xgk[0], xgk[2], ... abscissae to optimally extend the 25-point gauss rule */
1379
1380 static const double wgE[13] = /* weights of the 25-point gauss rule */
1381 {
1382  0.011393798501026287947902964113235,
1383  0.026354986615032137261901815295299,
1384  0.040939156701306312655623487711646,
1385  0.054904695975835191925936891540473,
1386  0.068038333812356917207187185656708,
1387  0.080140700335001018013234959669111,
1388  0.091028261982963649811497220702892,
1389  0.100535949067050644202206890392686,
1390  0.108519624474263653116093957050117,
1391  0.114858259145711648339325545869556,
1392  0.119455763535784772228178126512901,
1393  0.122242442990310041688959518945852,
1394  0.123176053726715451203902873079050
1395 };
1396
1397 static const double wgkE[26] = /* weights of the 51-point kronrod rule */
1398 {
1399  0.001987383892330315926507851882843,
1400  0.005561932135356713758040236901066,
1401  0.009473973386174151607207710523655,
1402  0.013236229195571674813656405846976,
1403  0.016847817709128298231516667536336,
1404  0.020435371145882835456568292235939,
1405  0.024009945606953216220092489164881,
1406  0.027475317587851737802948455517811,
1407  0.030792300167387488891109020215229,
1408  0.034002130274329337836748795229551,
1409  0.037116271483415543560330625367620,
1410  0.040083825504032382074839284467076,
1411  0.042872845020170049476895792439495,
1412  0.045502913049921788909870584752660,
1413  0.047982537138836713906392255756915,
1414  0.050277679080715671963325259433440,
1415  0.052362885806407475864366712137873,
1416  0.054251129888545490144543370459876,
1417  0.055950811220412317308240686382747,
1418  0.057437116361567832853582693939506,
1419  0.058689680022394207961974175856788,
1420  0.059720340324174059979099291932562,
1421  0.060539455376045862945360267517565,
1422  0.061128509717053048305859030416293,
1423  0.061471189871425316661544131965264,
1424  0.061580818067832935078759824240066
1425 };
1426
1427 /* wgk[25] was calculated from the values of wgk[0..24] */
1428
1429 void
1430 gsl_integration_qk51 (const gsl_function * f, double a, double b,
1431  double *result, double *abserr,
1432  double *resabs, double *resasc)
1433 {
1434  double fv1[26], fv2[26];
1435  //coverity[UNINIT]
1436  gsl_integration_qk (26, xgkE, wgE, wgkE, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
1437 }
1438
1439 /* Gauss quadrature weights and kronrod quadrature abscissae and
1440  weights as evaluated with 80 decimal digit arithmetic by
1441  L. W. Fullerton, Bell Labs, Nov. 1981. */
1442
1443 static const double xgkF[31] = /* abscissae of the 61-point kronrod rule */
1444 {
1445  0.999484410050490637571325895705811,
1446  0.996893484074649540271630050918695,
1447  0.991630996870404594858628366109486,
1448  0.983668123279747209970032581605663,
1449  0.973116322501126268374693868423707,
1450  0.960021864968307512216871025581798,
1451  0.944374444748559979415831324037439,
1452  0.926200047429274325879324277080474,
1453  0.905573307699907798546522558925958,
1454  0.882560535792052681543116462530226,
1455  0.857205233546061098958658510658944,
1456  0.829565762382768397442898119732502,
1457  0.799727835821839083013668942322683,
1458  0.767777432104826194917977340974503,
1459  0.733790062453226804726171131369528,
1460  0.697850494793315796932292388026640,
1461  0.660061064126626961370053668149271,
1462  0.620526182989242861140477556431189,
1463  0.579345235826361691756024932172540,
1464  0.536624148142019899264169793311073,
1465  0.492480467861778574993693061207709,
1466  0.447033769538089176780609900322854,
1467  0.400401254830394392535476211542661,
1468  0.352704725530878113471037207089374,
1469  0.304073202273625077372677107199257,
1470  0.254636926167889846439805129817805,
1471  0.204525116682309891438957671002025,
1472  0.153869913608583546963794672743256,
1473  0.102806937966737030147096751318001,
1474  0.051471842555317695833025213166723,
1475  0.000000000000000000000000000000000
1476 };
1477
1478 /* xgk[1], xgk[3], ... abscissae of the 30-point gauss rule.
1479  xgk[0], xgk[2], ... abscissae to optimally extend the 30-point gauss rule */
1480
1481 static const double wgF[15] = /* weights of the 30-point gauss rule */
1482 {
1483  0.007968192496166605615465883474674,
1484  0.018466468311090959142302131912047,
1485  0.028784707883323369349719179611292,
1486  0.038799192569627049596801936446348,
1487  0.048402672830594052902938140422808,
1488  0.057493156217619066481721689402056,
1489  0.065974229882180495128128515115962,
1490  0.073755974737705206268243850022191,
1491  0.080755895229420215354694938460530,
1492  0.086899787201082979802387530715126,
1493  0.092122522237786128717632707087619,
1494  0.096368737174644259639468626351810,
1495  0.099593420586795267062780282103569,
1496  0.101762389748405504596428952168554,
1497  0.102852652893558840341285636705415
1498 };
1499
1500 static const double wgkF[31] = /* weights of the 61-point kronrod rule */
1501 {
1502  0.001389013698677007624551591226760,
1503  0.003890461127099884051267201844516,
1504  0.006630703915931292173319826369750,
1505  0.009273279659517763428441146892024,
1506  0.011823015253496341742232898853251,
1507  0.014369729507045804812451432443580,
1508  0.016920889189053272627572289420322,
1509  0.019414141193942381173408951050128,
1510  0.021828035821609192297167485738339,
1511  0.024191162078080601365686370725232,
1512  0.026509954882333101610601709335075,
1513  0.028754048765041292843978785354334,
1514  0.030907257562387762472884252943092,
1515  0.032981447057483726031814191016854,
1516  0.034979338028060024137499670731468,
1517  0.036882364651821229223911065617136,
1518  0.038678945624727592950348651532281,
1519  0.040374538951535959111995279752468,
1520  0.041969810215164246147147541285970,
1521  0.043452539701356069316831728117073,
1522  0.044814800133162663192355551616723,
1523  0.046059238271006988116271735559374,
1524  0.047185546569299153945261478181099,
1525  0.048185861757087129140779492298305,
1526  0.049055434555029778887528165367238,
1527  0.049795683427074206357811569379942,
1528  0.050405921402782346840893085653585,
1529  0.050881795898749606492297473049805,
1530  0.051221547849258772170656282604944,
1531  0.051426128537459025933862879215781,
1532  0.051494729429451567558340433647099
1533 };
1534
1535 void
1536 gsl_integration_qk61 (const gsl_function * f, double a, double b,
1537  double *result, double *abserr,
1538  double *resabs, double *resasc)
1539 {
1540  double fv1[31], fv2[31];
1541  //coverity[UNINIT]
1542  gsl_integration_qk (31, xgkF, wgF, wgkF, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
1543 }
1544
1545 gsl_integration_workspace*
1547 {
1548  gsl_integration_workspace* w ;
1549
1550  if (n == 0)
1551  {
1552  GSL_ERROR_VAL ("workspace length n must be positive integer",
1553  GSL_EDOM, 0);
1554  }
1555
1556  w = (gsl_integration_workspace *)
1557  malloc (sizeof (gsl_integration_workspace));
1558
1559  if (w == 0)
1560  {
1561  GSL_ERROR_VAL ("failed to allocate space for workspace struct",
1562  GSL_ENOMEM, 0);
1563  }
1564
1565  w->alist = (double *) malloc (n * sizeof (double));
1566
1567  if (w->alist == 0)
1568  {
1569  free (w); /* exception in constructor, avoid memory leak */
1570
1571  GSL_ERROR_VAL ("failed to allocate space for alist ranges",
1572  GSL_ENOMEM, 0);
1573  }
1574
1575  w->blist = (double *) malloc (n * sizeof (double));
1576
1577  if (w->blist == 0)
1578  {
1579  free (w->alist);
1580  free (w); /* exception in constructor, avoid memory leak */
1581
1582  GSL_ERROR_VAL ("failed to allocate space for blist ranges",
1583  GSL_ENOMEM, 0);
1584  }
1585
1586  w->rlist = (double *) malloc (n * sizeof (double));
1587
1588  if (w->rlist == 0)
1589  {
1590  free (w->blist);
1591  free (w->alist);
1592  free (w); /* exception in constructor, avoid memory leak */
1593
1594  GSL_ERROR_VAL ("failed to allocate space for rlist ranges",
1595  GSL_ENOMEM, 0);
1596  }
1597
1598
1599  w->elist = (double *) malloc (n * sizeof (double));
1600
1601  if (w->elist == 0)
1602  {
1603  free (w->rlist);
1604  free (w->blist);
1605  free (w->alist);
1606  free (w); /* exception in constructor, avoid memory leak */
1607
1608  GSL_ERROR_VAL ("failed to allocate space for elist ranges",
1609  GSL_ENOMEM, 0);
1610  }
1611
1612  w->order = (size_t *) malloc (n * sizeof (size_t));
1613
1614  if (w->order == 0)
1615  {
1616  free (w->elist);
1617  free (w->rlist);
1618  free (w->blist);
1619  free (w->alist);
1620  free (w); /* exception in constructor, avoid memory leak */
1621
1622  GSL_ERROR_VAL ("failed to allocate space for order ranges",
1623  GSL_ENOMEM, 0);
1624  }
1625
1626  w->level = (size_t *) malloc (n * sizeof (size_t));
1627
1628  if (w->level == 0)
1629  {
1630  free (w->order);
1631  free (w->elist);
1632  free (w->rlist);
1633  free (w->blist);
1634  free (w->alist);
1635  free (w); /* exception in constructor, avoid memory leak */
1636
1637  GSL_ERROR_VAL ("failed to allocate space for order ranges",
1638  GSL_ENOMEM, 0);
1639  }
1640
1641  w->size = 0 ;
1642  w->limit = n ;
1643  w->maximum_level = 0 ;
1644
1645  return w ;
1646 }
1647
1648 void
1649 gsl_integration_workspace_free (gsl_integration_workspace * w)
1650 {
1651  free (w->level) ;
1652  free (w->order) ;
1653  free (w->elist) ;
1654  free (w->rlist) ;
1655  free (w->blist) ;
1656  free (w->alist) ;
1657  free (w) ;
1658 }
1659
1660
1661
1662 // INCLUDED BELOW #include "reset.c"
1663 static inline void
1664 reset_nrmax (gsl_integration_workspace * workspace);
1665
1666 static inline void
1667 reset_nrmax (gsl_integration_workspace * workspace)
1668 {
1669  workspace->nrmax = 0;
1670  workspace->i = workspace->order[0] ;
1671 }
1672
1673
1674 // INCLUDED BELOW #include "qpsrt2.c"
1675 /* The smallest interval has the largest error. Before bisecting
1676  decrease the sum of the errors over the larger intervals
1677  (error_over_large_intervals) and perform extrapolation. */
1678
1679 static int
1680 increase_nrmax (gsl_integration_workspace * workspace);
1681
1682 static int
1683 increase_nrmax (gsl_integration_workspace * workspace)
1684 {
1685  int k;
1686  int id = workspace->nrmax;
1687  int jupbnd;
1688
1689  const size_t * level = workspace->level;
1690  const size_t * order = workspace->order;
1691
1692  size_t limit = workspace->limit ;
1693  size_t last = workspace->size - 1 ;
1694
1695  if (last > (1 + limit / 2))
1696  {
1697  jupbnd = limit + 1 - last;
1698  }
1699  else
1700  {
1701  jupbnd = last;
1702  }
1703
1704  for (k = id; k <= jupbnd; k++)
1705  {
1706  size_t i_max = order[workspace->nrmax];
1707
1708  workspace->i = i_max ;
1709
1710  if (level[i_max] < workspace->maximum_level)
1711  {
1712  return 1;
1713  }
1714
1715  workspace->nrmax++;
1716
1717  }
1718  return 0;
1719 }
1720
1721 static int
1722 large_interval (gsl_integration_workspace * workspace)
1723 {
1724  size_t i = workspace->i ;
1725  const size_t * level = workspace->level;
1726
1727  if (level[i] < workspace->maximum_level)
1728  {
1729  return 1 ;
1730  }
1731  else
1732  {
1733  return 0 ;
1734  }
1735 }
1736
1737
1738 // INCLUDED BELOW #include "qelg.c"
1739 struct extrapolation_table
1740  {
1741  size_t n;
1742  double rlist2[52];
1743  size_t nres;
1744  double res3la[3];
1745  };
1746
1747 static void
1748  initialise_table (struct extrapolation_table *table);
1749
1750 static void
1751  append_table (struct extrapolation_table *table, double y);
1752
1753 static void
1754 initialise_table (struct extrapolation_table *table)
1755 {
1756  table->n = 0;
1757  table->nres = 0;
1758 }
1759 #ifdef JUNK
1760 static void
1761 initialise_table (struct extrapolation_table *table, double y)
1762 {
1763  table->n = 0;
1764  table->rlist2[0] = y;
1765  table->nres = 0;
1766 }
1767 #endif
1768 static void
1769 append_table (struct extrapolation_table *table, double y)
1770 {
1771  size_t n;
1772  n = table->n;
1773  table->rlist2[n] = y;
1774  table->n++;
1775 }
1776
1777 /* static inline void
1778  qelg (size_t * n, double epstab[],
1779  double * result, double * abserr,
1780  double res3la[], size_t * nres); */
1781
1782 static inline void
1783  qelg (struct extrapolation_table *table, double *result, double *abserr);
1784
1785 static inline void
1786 qelg (struct extrapolation_table *table, double *result, double *abserr)
1787 {
1788  double *epstab = table->rlist2;
1789  double *res3la = table->res3la;
1790  const size_t n = table->n - 1;
1791
1792  const double current = epstab[n];
1793
1794  double absolute = GSL_DBL_MAX;
1795  double relative = 5 * GSL_DBL_EPSILON * fabs (current);
1796
1797  const size_t newelm = n / 2;
1798  const size_t n_orig = n;
1799  size_t n_final = n;
1800  size_t i;
1801
1802  const size_t nres_orig = table->nres;
1803
1804  *result = current;
1805  *abserr = GSL_DBL_MAX;
1806
1807  if (n < 2)
1808  {
1809  *result = current;
1810  *abserr = GSL_MAX_DBL (absolute, relative);
1811  return;
1812  }
1813
1814  epstab[n + 2] = epstab[n];
1815  epstab[n] = GSL_DBL_MAX;
1816
1817  for (i = 0; i < newelm; i++)
1818  {
1819  double res = epstab[n - 2 * i + 2];
1820  double e0 = epstab[n - 2 * i - 2];
1821  double e1 = epstab[n - 2 * i - 1];
1822  double e2 = res;
1823
1824  double e1abs = fabs (e1);
1825  double delta2 = e2 - e1;
1826  double err2 = fabs (delta2);
1827  double tol2 = GSL_MAX_DBL (fabs (e2), e1abs) * GSL_DBL_EPSILON;
1828  double delta3 = e1 - e0;
1829  double err3 = fabs (delta3);
1830  double tol3 = GSL_MAX_DBL (e1abs, fabs (e0)) * GSL_DBL_EPSILON;
1831
1832  double e3, delta1, err1, tol1, ss;
1833
1834  if (err2 <= tol2 && err3 <= tol3)
1835  {
1836  /* If e0, e1 and e2 are equal to within machine accuracy,
1837  convergence is assumed. */
1838
1839  *result = res;
1840  absolute = err2 + err3;
1841  relative = 5 * GSL_DBL_EPSILON * fabs (res);
1842  *abserr = GSL_MAX_DBL (absolute, relative);
1843  return;
1844  }
1845
1846  e3 = epstab[n - 2 * i];
1847  epstab[n - 2 * i] = e1;
1848  delta1 = e1 - e3;
1849  err1 = fabs (delta1);
1850  tol1 = GSL_MAX_DBL (e1abs, fabs (e3)) * GSL_DBL_EPSILON;
1851
1852  /* If two elements are very close to each other, omit a part of
1853  the table by adjusting the value of n */
1854
1855  if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3)
1856  {
1857  n_final = 2 * i;
1858  break;
1859  }
1860
1861  ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
1862
1863  /* Test to detect irregular behaviour in the table, and
1864  eventually omit a part of the table by adjusting the value of
1865  n. */
1866
1867  if (fabs (ss * e1) <= 0.0001)
1868  {
1869  n_final = 2 * i;
1870  break;
1871  }
1872
1873  /* Compute a new element and eventually adjust the value of
1874  result. */
1875
1876  res = e1 + 1 / ss;
1877  epstab[n - 2 * i] = res;
1878
1879  {
1880  const double error = err2 + fabs (res - e2) + err3;
1881
1882  if (error <= *abserr)
1883  {
1884  *abserr = error;
1885  *result = res;
1886  }
1887  }
1888  }
1889
1890  /* Shift the table */
1891
1892  {
1893  const size_t limexp = 50 - 1;
1894
1895  if (n_final == limexp)
1896  {
1897  n_final = 2 * (limexp / 2);
1898  }
1899  }
1900
1901  if (n_orig % 2 == 1)
1902  {
1903  for (i = 0; i <= newelm; i++)
1904  {
1905  epstab[1 + i * 2] = epstab[i * 2 + 3];
1906  }
1907  }
1908  else
1909  {
1910  for (i = 0; i <= newelm; i++)
1911  {
1912  epstab[i * 2] = epstab[i * 2 + 2];
1913  }
1914  }
1915
1916  if (n_orig != n_final)
1917  {
1918  for (i = 0; i <= n_final; i++)
1919  {
1920  epstab[i] = epstab[n_orig - n_final + i];
1921  }
1922  }
1923
1924  table->n = n_final + 1;
1925
1926  if (nres_orig < 3)
1927  {
1928  res3la[nres_orig] = *result;
1929  *abserr = GSL_DBL_MAX;
1930  }
1931  else
1932  { /* Compute error estimate */
1933  *abserr = (fabs (*result - res3la[2]) + fabs (*result - res3la[1])
1934  + fabs (*result - res3la[0]));
1935
1936  res3la[0] = res3la[1];
1937  res3la[1] = res3la[2];
1938  res3la[2] = *result;
1939  }
1940
1941  /* In QUADPACK the variable table->nres is incremented at the top of
1942  qelg, so it increases on every call. This leads to the array
1943  res3la being accessed when its elements are still undefined, so I
1944  have moved the update to this point so that its value more
1945  useful. */
1946
1947  table->nres = nres_orig + 1;
1948
1949  *abserr = GSL_MAX_DBL (*abserr, 5 * GSL_DBL_EPSILON * fabs (*result));
1950
1951  return;
1952 }
1953
1954
1955 // INCLUDED BELOW #include "positivity.c"
1956 /* Compare the integral of f(x) with the integral of |f(x)|
1957  to determine if f(x) covers both positive and negative values */
1958
1959 static inline int
1960 test_positivity (double result, double resabs);
1961
1962 static inline int
1963 test_positivity (double result, double resabs)
1964 {
1965  int status = (fabs (result) >= (1 - 50 * GSL_DBL_EPSILON) * resabs);
1966
1967  return status;
1968 }
1969
1970 static int qags (const gsl_function * f, const double a, const double
1971  b, const double epsabs, const double epsrel, const size_t limit,
1972  gsl_integration_workspace * workspace, double *result, double *abserr,
1974
1975 int
1977  double a, double b,
1978  double epsabs, double epsrel, size_t limit,
1979  gsl_integration_workspace * workspace,
1980  double * result, double * abserr)
1981 {
1982  int status = qags (f, a, b, epsabs, epsrel, limit,
1983  workspace,
1984  result, abserr,
1986  return status ;
1987 }
1988
1989 /* QAGI: evaluate an integral over an infinite range using the
1990  transformation
1991
1992  integrate(f(x),-Inf,Inf) = integrate((f((1-t)/t) + f(-(1-t)/t))/t^2,0,1)
1993
1994  */
1995
1996 static double i_transform (double t, void *params);
1997
1998 int
2000  double epsabs, double epsrel, size_t limit,
2001  gsl_integration_workspace * workspace,
2002  double *result, double *abserr)
2003 {
2004  int status;
2005
2006  gsl_function f_transform;
2007
2008  f_transform.function = &i_transform;
2009  f_transform.params = f;
2010
2011  status = qags (&f_transform, 0.0, 1.0,
2012  epsabs, epsrel, limit,
2013  workspace,
2014  result, abserr,
2016
2017  return status;
2018 }
2019
2020 static double
2021 i_transform (double t, void *params)
2022 {
2023  gsl_function *f = (gsl_function *) params;
2024  double x = (1 - t) / t;
2025  double y = GSL_FN_EVAL (f, x) + GSL_FN_EVAL (f, -x);
2026  return (y / t) / t;
2027 }
2028
2029
2030 /* QAGIL: Evaluate an integral over an infinite range using the
2031  transformation,
2032
2033  integrate(f(x),-Inf,b) = integrate(f(b-(1-t)/t)/t^2,0,1)
2034
2035  */
2036
2037 struct il_params { double b ; gsl_function * f ; } ;
2038
2039 static double il_transform (double t, void *params);
2040
2041 int
2043  double b,
2044  double epsabs, double epsrel, size_t limit,
2045  gsl_integration_workspace * workspace,
2046  double *result, double *abserr)
2047 {
2048  int status;
2049
2050  gsl_function f_transform;
2051  struct il_params transform_params ;
2052
2053  transform_params.b = b ;
2054  transform_params.f = f ;
2055
2056  f_transform.function = &il_transform;
2057  f_transform.params = &transform_params;
2058
2059  status = qags (&f_transform, 0.0, 1.0,
2060  epsabs, epsrel, limit,
2061  workspace,
2062  result, abserr,
2064
2065  return status;
2066 }
2067
2068 static double
2069 il_transform (double t, void *params)
2070 {
2071  struct il_params *p = (struct il_params *) params;
2072  double b = p->b;
2073  gsl_function * f = p->f;
2074  double x = b - (1 - t) / t;
2075  double y = GSL_FN_EVAL (f, x);
2076  return (y / t) / t;
2077 }
2078
2079 /* QAGIU: Evaluate an integral over an infinite range using the
2080  transformation
2081
2082  integrate(f(x),a,Inf) = integrate(f(a+(1-t)/t)/t^2,0,1)
2083
2084  */
2085
2086 struct iu_params { double a ; gsl_function * f ; } ;
2087
2088 static double iu_transform (double t, void *params);
2089
2090 int
2092  double a,
2093  double epsabs, double epsrel, size_t limit,
2094  gsl_integration_workspace * workspace,
2095  double *result, double *abserr)
2096 {
2097  int status;
2098
2099  gsl_function f_transform;
2100  struct iu_params transform_params ;
2101
2102  transform_params.a = a ;
2103  transform_params.f = f ;
2104
2105  f_transform.function = &iu_transform;
2106  f_transform.params = &transform_params;
2107
2108  status = qags (&f_transform, 0.0, 1.0,
2109  epsabs, epsrel, limit,
2110  workspace,
2111  result, abserr,
2113
2114  return status;
2115 }
2116
2117 static double
2118 iu_transform (double t, void *params)
2119 {
2120  struct iu_params *p = (struct iu_params *) params;
2121  double a = p->a;
2122  gsl_function * f = p->f;
2123  double x = a + (1 - t) / t;
2124  double y = GSL_FN_EVAL (f, x);
2125  return (y / t) / t;
2126 }
2127
2128 /* Main integration function */
2129
2130 static int
2132  const double a, const double b,
2133  const double epsabs, const double epsrel,
2134  const size_t limit,
2135  gsl_integration_workspace * workspace,
2136  double *result, double *abserr,
2138 {
2139  double area, errsum;
2140  double res_ext, err_ext;
2141  double result0, abserr0, resabs0, resasc0;
2142  double tolerance;
2143
2144  double ertest = 0;
2145  double error_over_large_intervals = 0;
2146  double reseps = 0, abseps = 0, correc = 0;
2147  size_t ktmin = 0;
2148  int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
2149  int error_type = 0, error_type2 = 0;
2150
2151  size_t iteration = 0;
2152
2153  int positive_integrand = 0;
2154  int extrapolate = 0;
2155  int disallow_extrapolation = 0;
2156
2157  struct extrapolation_table table;
2158
2159  /* Initialize results */
2160
2161  initialise (workspace, a, b);
2162
2163  *result = 0;
2164  *abserr = 0;
2165
2166  if (limit > workspace->limit)
2167  {
2168  GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
2169  }
2170
2171  /* Test on accuracy */
2172
2173  if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
2174  {
2175  GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
2176  GSL_EBADTOL);
2177  }
2178
2179  /* Perform the first integration */
2180
2181  q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
2182
2183  set_initial_result (workspace, result0, abserr0);
2184
2185  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
2186
2187  if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance)
2188  {
2189  *result = result0;
2190  *abserr = abserr0;
2191
2192  GSL_ERROR ("cannot reach tolerance because of roundoff error"
2193  "on first attempt", GSL_EROUND);
2194  }
2195  else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
2196  {
2197  *result = result0;
2198  *abserr = abserr0;
2199
2200  return GSL_SUCCESS;
2201  }
2202  else if (limit == 1)
2203  {
2204  *result = result0;
2205  *abserr = abserr0;
2206
2207  GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
2208  }
2209
2210  /* Initialization */
2211
2212  initialise_table (&table);
2213  append_table (&table, result0);
2214
2215  area = result0;
2216  errsum = abserr0;
2217
2218  res_ext = result0;
2219  err_ext = GSL_DBL_MAX;
2220
2221  positive_integrand = test_positivity (result0, resabs0);
2222
2223  iteration = 1;
2224
2225  do
2226  {
2227  size_t current_level;
2228  double a1, b1, a2, b2;
2229  double a_i, b_i, r_i, e_i;
2230  double area1 = 0, area2 = 0, area12 = 0;
2231  double error1 = 0, error2 = 0, error12 = 0;
2232  double resasc1, resasc2;
2233  double resabs1, resabs2;
2234  double last_e_i;
2235
2236  /* Bisect the subinterval with the largest error estimate */
2237
2238  retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
2239
2240  current_level = workspace->level[workspace->i] + 1;
2241
2242  a1 = a_i;
2243  b1 = 0.5 * (a_i + b_i);
2244  a2 = b1;
2245  b2 = b_i;
2246
2247  iteration++;
2248
2249  q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
2250  q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
2251
2252  area12 = area1 + area2;
2253  error12 = error1 + error2;
2254  last_e_i = e_i;
2255
2256  /* Improve previous approximations to the integral and test for
2257  accuracy.
2258
2259  We write these expressions in the same way as the original
2260  QUADPACK code so that the rounding errors are the same, which
2261  makes testing easier. */
2262
2263  errsum = errsum + error12 - e_i;
2264  area = area + area12 - r_i;
2265
2266  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
2267
2268  if (resasc1 != error1 && resasc2 != error2)
2269  {
2270  double delta = r_i - area12;
2271
2272  if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
2273  {
2274  if (!extrapolate)
2275  {
2276  roundoff_type1++;
2277  }
2278  else
2279  {
2280  roundoff_type2++;
2281  }
2282  }
2283  if (iteration > 10 && error12 > e_i)
2284  {
2285  roundoff_type3++;
2286  }
2287  }
2288
2289  /* Test for roundoff and eventually set error flag */
2290
2291  if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
2292  {
2293  error_type = 2; /* round off error */
2294  }
2295
2296  if (roundoff_type2 >= 5)
2297  {
2298  error_type2 = 1;
2299  }
2300
2301  /* set error flag in the case of bad integrand behaviour at
2302  a point of the integration range */
2303
2304  if (subinterval_too_small (a1, a2, b2))
2305  {
2306  error_type = 4;
2307  }
2308
2309  /* append the newly-created intervals to the list */
2310
2311  update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
2312
2313  if (errsum <= tolerance)
2314  {
2315  goto compute_result;
2316  }
2317
2318  if (error_type)
2319  {
2320  break;
2321  }
2322
2323  if (iteration >= limit - 1)
2324  {
2325  error_type = 1;
2326  break;
2327  }
2328
2329  if (iteration == 2) /* set up variables on first iteration */
2330  {
2331  error_over_large_intervals = errsum;
2332  ertest = tolerance;
2333  append_table (&table, area);
2334  continue;
2335  }
2336
2337  if (disallow_extrapolation)
2338  {
2339  continue;
2340  }
2341
2342  error_over_large_intervals += -last_e_i;
2343
2344  if (current_level < workspace->maximum_level)
2345  {
2346  error_over_large_intervals += error12;
2347  }
2348
2349  if (!extrapolate)
2350  {
2351  /* test whether the interval to be bisected next is the
2352  smallest interval. */
2353
2354  if (large_interval (workspace))
2355  continue;
2356
2357  extrapolate = 1;
2358  workspace->nrmax = 1;
2359  }
2360
2361  if (!error_type2 && error_over_large_intervals > ertest)
2362  {
2363  if (increase_nrmax (workspace))
2364  continue;
2365  }
2366
2367  /* Perform extrapolation */
2368
2369  append_table (&table, area);
2370
2371  qelg (&table, &reseps, &abseps);
2372
2373  ktmin++;
2374
2375  if (ktmin > 5 && err_ext < 0.001 * errsum)
2376  {
2377  error_type = 5;
2378  }
2379
2380  if (abseps < err_ext)
2381  {
2382  ktmin = 0;
2383  err_ext = abseps;
2384  res_ext = reseps;
2385  correc = error_over_large_intervals;
2386  ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
2387  if (err_ext <= ertest)
2388  break;
2389  }
2390
2391  /* Prepare bisection of the smallest interval. */
2392
2393  if (table.n == 1)
2394  {
2395  disallow_extrapolation = 1;
2396  }
2397
2398  if (error_type == 5)
2399  {
2400  break;
2401  }
2402
2403  /* work on interval with largest error */
2404
2405  reset_nrmax (workspace);
2406  extrapolate = 0;
2407  error_over_large_intervals = errsum;
2408
2409  }
2410  while (iteration < limit);
2411
2412  *result = res_ext;
2413  *abserr = err_ext;
2414
2415  if (err_ext == GSL_DBL_MAX)
2416  goto compute_result;
2417
2418  if (error_type || error_type2)
2419  {
2420  if (error_type2)
2421  {
2422  err_ext += correc;
2423  }
2424
2425 // if (error_type == 0)
2426 // error_type = 3;
2427
2428  if (res_ext != 0.0 && area != 0.0)
2429  {
2430  if (err_ext / fabs (res_ext) > errsum / fabs (area))
2431  goto compute_result;
2432  }
2433  else if (err_ext > errsum)
2434  {
2435  goto compute_result;
2436  }
2437  else if (area == 0.0)
2438  {
2439  goto return_error;
2440  }
2441  }
2442
2443  /* Test on divergence. */
2444
2445  {
2446  double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
2447
2448  if (!positive_integrand && max_area < 0.01 * resabs0)
2449  goto return_error;
2450  }
2451
2452  {
2453  double ratio = res_ext / area;
2454
2455  if (ratio < 0.01 || ratio > 100.0 || errsum > fabs (area))
2456  error_type = 6;
2457  }
2458
2459  goto return_error;
2460
2461 compute_result:
2462
2463  *result = sum_results (workspace);
2464  *abserr = errsum;
2465
2466 return_error:
2467
2468  if (error_type > 2)
2469  error_type--;
2470
2471
2472
2473  if (error_type == 0)
2474  {
2475  return GSL_SUCCESS;
2476  }
2477  else if (error_type == 1)
2478  {
2479  GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
2480  }
2481  else if (error_type == 2)
2482  {
2483  GSL_ERROR ("cannot reach tolerance because of roundoff error",
2484  GSL_EROUND);
2485  }
2486  else if (error_type == 3)
2487  {
2488  GSL_ERROR ("bad integrand behavior found in the integration interval",
2489  GSL_ESING);
2490  }
2491  else if (error_type == 4)
2492  {
2493  GSL_ERROR ("roundoff error detected in the extrapolation table",
2494  GSL_EROUND);
2495  }
2496  else if (error_type == 5)
2497  {
2498  GSL_ERROR ("integral is divergent, or slowly convergent",
2499  GSL_EDIVERGE);
2500  }
2501
2502  GSL_ERROR ("could not integrate function", GSL_EFAILED);
2503 }
static double sum_results(const gsl_integration_workspace *workspace)
static const double xgkC[16]
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
Int_t getCatIndex(const char *name, Int_t defVal=0, Bool_t verbose=kFALSE) const
Get index value of a RooAbsCategory stored in set with given name.
Definition: RooArgSet.cxx:613
#define coutE(a)
Definition: RooMsgService.h:34
float xmin
Definition: THbookFile.cxx:93
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
static double il_transform(double t, void *params)
double gsl_coerce_double(const double x)
static const double wgF[15]
int gsl_integration_qagiu(gsl_function *f, double a, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
virtual Bool_t setIndex(Int_t index, Bool_t printError=kTRUE)
Set value by specifying the index code of the desired state.
RooNumIntFactory is a factory to instantiate numeric integrators from a given function binding and a ...
static const double wgkB[11]
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
static const double wgkD[21]
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
STL namespace.
void gsl_integration_qk51(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
int gsl_integration_qag(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, int key, gsl_integration_workspace *workspace, double *result, double *abserr)
static void qelg(struct extrapolation_table *table, double *result, double *abserr)
Bool_t isValid() const
#define malloc
Definition: civetweb.c:818
virtual Bool_t checkLimits() const
Check that our integration range is finite and otherwise return kFALSE.
static void retrieve(const gsl_integration_workspace *workspace, double *a, double *b, double *r, double *e)
Float_t delta
static void append_table(struct extrapolation_table *table, double y)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:628
static const double wgkC[16]
static const double wgkA[8]
Double_t x[n]
Definition: legend1.C:17
Double_t integrand(const Double_t x[]) const
static const double wgA[4]
void gsl_integration_rule(const gsl_function *f, double a, double b, double *result, double *abserr, double *defabs, double *resabs)
int gsl_integration_qags(const gsl_function *f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
#define GSL_ERROR(a, b)
#define GSL_FN_EVAL(F, x)
static int subinterval_too_small(double a1, double a2, double b2)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:146
static int increase_nrmax(gsl_integration_workspace *workspace)
#define GSL_ERROR_VAL(reason, gsl_errno, value)
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
#define F(x, y, z)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
TRandom2 r(17)
static void reset_nrmax(gsl_integration_workspace *workspace)
static int large_interval(gsl_integration_workspace *workspace)
void gsl_integration_workspace_free(gsl_integration_workspace *w)
UInt_t getDimension() const
Definition: RooAbsFunc.h:29
void gsl_integration_qcheb(gsl_function *f, double a, double b, double *cheb12, double *cheb24)
static void update(gsl_integration_workspace *workspace, double a1, double b1, double area1, double error1, double a2, double b2, double area2, double error2)
static void qpsrt(gsl_integration_workspace *workspace)
unsigned int UInt_t
Definition: RtypesCore.h:42
const RooAbsFunc * _function
static const double wgkE[26]
void gsl_integration_qk(const int n, const double xgk[], const double wg[], const double wgk[], double fv1[], double fv2[], const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
virtual Double_t getMinLimit(UInt_t dimension) const =0
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
void gsl_integration_qk41(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
void gsl_integration_qk61(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
float xmax
Definition: THbookFile.cxx:93
static void initialise_table(struct extrapolation_table *table)
RooCategory represents a fundamental (non-derived) discrete value object.
Definition: RooCategory.h:24
double GSL_MAX_DBL(double a, double b)
static void initialise(gsl_integration_workspace *workspace, double a, double b)
static void set_initial_result(gsl_integration_workspace *workspace, double result, double error)
const Bool_t kFALSE
Definition: RtypesCore.h:92
const RooAbsFunc * integrand() const
static const double wgB[5]
#define ClassImp(name)
Definition: Rtypes.h:336
double f(double x)
struct gsl_function_struct gsl_function
virtual Double_t getMaxLimit(UInt_t dimension) const =0
double Double_t
Definition: RtypesCore.h:55
static double iu_transform(double t, void *params)
#define free
Definition: civetweb.c:821
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Virtual constructor.
static int test_positivity(double result, double resabs)
Double_t y[n]
Definition: legend1.C:17
static int qag(const gsl_function *f, const double a, const double b, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr, gsl_integration_rule *q)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
virtual Double_t integral(const Double_t *yvec=0)
Calculate and return integral at at given parameter values.
static const double wgE[13]
int gsl_integration_qagil(gsl_function *f, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
static const double wgC[8]
static const double xgkD[21]
static int qags(const gsl_function *f, const double a, const double b, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr, gsl_integration_rule *q)
static const double xgkE[26]
static const double xgkB[11]
void gsl_integration_qk21(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
static void registerIntegrator(RooNumIntFactory &fact)
Register this class with RooNumIntConfig as a possible choice of numeric integrator for one-dimension...
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Definition: RooArgSet.cxx:527
void gsl_integration_qk31(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
friend double RooAdaptiveGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
Glue function interacing to GSL code.
static const double xgkF[31]
int gsl_integration_qagi(gsl_function *f, double epsabs, double epsrel, size_t limit, gsl_integration_workspace *workspace, double *result, double *abserr)
Bool_t defineType(const char *label)
Define a state with given name, the lowest available positive integer is assigned as index...
static const double xgkA[8]
double result[121]
static double i_transform(double t, void *params)
void gsl_integration_qk15(const gsl_function *f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
RooAdaptiveGaussKronrodIntegrator1D implements the Gauss-Kronrod integration algorithm.
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:364
static const double wgD[11]
float * q
Definition: THbookFile.cxx:87
const Bool_t kTRUE
Definition: RtypesCore.h:91
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
const Int_t n
Definition: legend1.C:16
gsl_integration_workspace * gsl_integration_workspace_alloc(const size_t n)
static const double wgkF[31]
static double rescale_error(double err, const double result_abs, const double result_asc)
Bool_t storeProtoIntegrator(RooAbsIntegrator *proto, const RooArgSet &defConfig, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
RooAdaptiveGaussKronrodIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
Bool_t initialize()
Initialize integrator allocate buffers and setup GSL workspace.