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