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