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