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