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