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