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