TMB Documentation  v1.9.11
integrate.cpp
1 /*
2  * R : A Computer Language for Statistical Data Analysis
3  * Copyright (C) 2001-2014 The R Core Team
4  *
5  * This program is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation; either version 2 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, a copy is available at
17  * https://www.R-project.org/Licenses/
18  *
19  * Most of this file is C translations of Fortran routines in
20  * QUADPACK: the latter is part of SLATEC 'and therefore in the public
21  * domain' (https://en.wikipedia.org/wiki/QUADPACK).
22  *
23  *
24  */
25 
26 #ifdef HAVE_CONFIG_H
27 #include <config.h>
28 #endif
29 
30 #include <math.h>
31 #include <float.h>
32 //#include <Rmath.h> /* for fmax2, fmin2, imin2 */
33 //#include <R_ext/Applic.h> /* exporting the API , particularly */
34 /*--- typedef void integr_fn(double *x, int n, void *ex) ---
35  * vectorizing function f(x[1:n], ...) -> x[] {overwriting x[]}.
36  * Vectorization can be used to speed up the integrand
37  * instead of calling it n times.
38 */
39 
40 /* f2c-ed translations + modifications of QUADPACK functions */
41 
42 template<class Float, class integr_fn> static void rdqagie(integr_fn f, void *ex,
43  Float *, int *, Float * , Float *, int *,
44  Float *, Float *, int *,
45  int *, Float *, Float *, Float *, Float *,
46  int *, int *);
47 
48 template<class Float, class integr_fn> static void rdqk15i(integr_fn f, void *ex,
49  Float *, int *, Float * , Float *,
50  Float *, Float *, Float *, Float *);
51 
52 template<class Float, class integr_fn> static void rdqagse(integr_fn f, void *ex, Float *, Float *,
53  Float *, Float *, int *, Float *, Float *,
54  int *, int *, Float *, Float *, Float *,
55  Float *, int *, int *);
56 
57 template<class Float, class integr_fn> static void rdqk21(integr_fn f, void *ex,
58  Float *, Float *, Float *, Float *, Float *, Float *);
59 
60 template<class Float> static void rdqpsrt(int *, int *, int *, Float *, Float *, int *, int *);
61 
62 template<class Float> static void rdqelg(int *, Float *, Float *, Float *, Float *, int *);
63 
64 template<class Float, class integr_fn>
65 void Rdqagi(integr_fn f, void *ex, Float *bound, int *inf,
66  Float *epsabs, Float *epsrel,
67  Float *result, Float *abserr, int *neval, int *ier,
68  int *limit, int *lenw, int *last,
69  int *iwork, Float *work)
70 {
71  int l1, l2, l3;
72 
73 /*
74 ***begin prologue dqagi
75 ***date written 800101 (yymmdd)
76 ***revision date 830518 (yymmdd)
77 ***category no. h2a3a1,h2a4a1
78 ***keywords automatic integrator, infinite intervals,
79  general-purpose, transformation, extrapolation,
80  globally adaptive
81 ***author piessens,robert,appl. math. & progr. div. - k.u.leuven
82  de doncker,elise,appl. math. & progr. div. -k.u.leuven
83 ***purpose the routine calculates an approximation result to a given
84  integral i = integral of f over (bound,+infinity)
85  or i = integral of f over (-infinity,bound)
86  or i = integral of f over (-infinity,+infinity)
87  hopefully satisfying following claim for accuracy
88  abs(i-result) <= max(epsabs,epsrel*abs(i)).
89 ***description
90 
91  integration over infinite intervals
92  standard fortran subroutine
93 
94  parameters
95  on entry
96  f - double precision
97  function subprogram defining the integrand
98  function f(x). the actual name for f needs to be
99  declared e x t e r n a l in the driver program.
100 
101  bound - double precision
102  finite bound of integration range
103  (has no meaning if interval is doubly-infinite)
104 
105  inf - int
106  indicating the kind of integration range involved
107  inf = 1 corresponds to (bound,+infinity),
108  inf = -1 to (-infinity,bound),
109  inf = 2 to (-infinity,+infinity).
110 
111  epsabs - double precision
112  absolute accuracy requested
113  epsrel - double precision
114  relative accuracy requested
115  if epsabs <= 0
116  and epsrel < max(50*rel.mach.acc.,0.5d-28),
117  the routine will end with ier = 6.
118 
119 
120  on return
121  result - double precision
122  approximation to the integral
123 
124  abserr - double precision
125  estimate of the modulus of the absolute error,
126  which should equal or exceed abs(i-result)
127 
128  neval - int
129  number of integrand evaluations
130 
131  ier - int
132  ier = 0 normal and reliable termination of the
133  routine. it is assumed that the requested
134  accuracy has been achieved.
135  - ier > 0 abnormal termination of the routine. the
136  estimates for result and error are less
137  reliable. it is assumed that the requested
138  accuracy has not been achieved.
139  error messages
140  ier = 1 maximum number of subdivisions allowed
141  has been achieved. one can allow more
142  subdivisions by increasing the value of
143  limit (and taking the according dimension
144  adjustments into account). however, if
145  this yields no improvement it is advised
146  to analyze the integrand in order to
147  determine the integration difficulties. if
148  the position of a local difficulty can be
149  determined (e.g. singularity,
150  discontinuity within the interval) one
151  will probably gain from splitting up the
152  interval at this point and calling the
153  integrator on the subranges. if possible,
154  an appropriate special-purpose integrator
155  should be used, which is designed for
156  handling the type of difficulty involved.
157  = 2 the occurrence of roundoff error is
158  detected, which prevents the requested
159  tolerance from being achieved.
160  the error may be under-estimated.
161  = 3 extremely bad integrand behaviour occurs
162  at some points of the integration
163  interval.
164  = 4 the algorithm does not converge.
165  roundoff error is detected in the
166  extrapolation table.
167  it is assumed that the requested tolerance
168  cannot be achieved, and that the returned
169  result is the best which can be obtained.
170  = 5 the integral is probably divergent, or
171  slowly convergent. it must be noted that
172  divergence can occur with any other value
173  of ier.
174  = 6 the input is invalid, because
175  (epsabs <= 0 and
176  epsrel < max(50*rel.mach.acc.,0.5d-28))
177  or limit < 1 or leniw < limit*4.
178  result, abserr, neval, last are set to
179  zero. exept when limit or leniw is
180  invalid, iwork(1), work(limit*2+1) and
181  work(limit*3+1) are set to zero, work(1)
182  is set to a and work(limit+1) to b.
183 
184  dimensioning parameters
185  limit - int
186  dimensioning parameter for iwork
187  limit determines the maximum number of subintervals
188  in the partition of the given integration interval
189  (a,b), limit >= 1.
190  if limit < 1, the routine will end with ier = 6.
191 
192  lenw - int
193  dimensioning parameter for work
194  lenw must be at least limit*4.
195  if lenw < limit*4, the routine will end
196  with ier = 6.
197 
198  last - int
199  on return, last equals the number of subintervals
200  produced in the subdivision process, which
201  determines the number of significant elements
202  actually in the work arrays.
203 
204  work arrays
205  iwork - int
206  vector of dimension at least limit, the first
207  k elements of which contain pointers
208  to the error estimates over the subintervals,
209  such that work(limit*3+iwork(1)),... ,
210  work(limit*3+iwork(k)) form a decreasing
211  sequence, with k = last if last <= (limit/2+2), and
212  k = limit+1-last otherwise
213 
214  work - double precision
215  vector of dimension at least lenw
216  on return
217  work(1), ..., work(last) contain the left
218  end points of the subintervals in the
219  partition of (a,b),
220  work(limit+1), ..., work(limit+last) contain
221  the right end points,
222  work(limit*2+1), ...,work(limit*2+last) contain the
223  integral approximations over the subintervals,
224  work(limit*3+1), ..., work(limit*3)
225  contain the error estimates.
226 
227 ***routines called dqagie
228 ***end prologue dqagi */
229 
230  *ier = 6;
231  *neval = 0;
232  *last = 0;
233  *result = 0.;
234  *abserr = 0.;
235  if (*limit < 1 || *lenw < *limit << 2) return;
236 
237  l1 = *limit;
238  l2 = *limit + l1;
239  l3 = *limit + l2;
240 
241  rdqagie(f, ex, bound, inf, epsabs, epsrel, limit, result, abserr, neval, ier,
242  work, &work[l1], &work[l2], &work[l3], iwork, last);
243 
244  return;
245 } /* Rdqagi */
246 
247 template<class Float, class integr_fn> static
248 void rdqagie(integr_fn f, void *ex, Float *bound, int *inf, Float *
249  epsabs, Float *epsrel, int *limit, Float *result,
250  Float *abserr, int *neval, int *ier, Float *alist,
251  Float *blist, Float *rlist, Float *elist, int *
252  iord, int *last)
253 {
254  /* Local variables */
255  Float area, dres;
256  int ksgn;
257  Float boun;
258  int nres;
259  Float area1, area2, area12;
260  int k;
261  Float small = 0.0, erro12;
262  int ierro;
263  Float a1, a2, b1, b2, defab1, defab2, oflow;
264  int ktmin, nrmax;
265  Float uflow;
266  bool noext;
267  int iroff1, iroff2, iroff3;
268  Float res3la[3], error1, error2;
269  int id;
270  Float rlist2[52];
271  int numrl2;
272  Float defabs, epmach, erlarg = 0.0, abseps, correc = 0.0, errbnd, resabs;
273  int jupbnd;
274  Float erlast, errmax;
275  int maxerr;
276  Float reseps;
277  bool extrap;
278  Float ertest = 0.0, errsum;
279 
280 /* begin prologue dqagie
281 ***date written 800101 (yymmdd)
282 ***revision date 830518 (yymmdd)
283 ***category no. h2a3a1,h2a4a1
284 ***keywords automatic integrator, infinite intervals,
285  general-purpose, transformation, extrapolation,
286  globally adaptive
287 ***author piessens,robert,appl. math & progr. div - k.u.leuven
288  de doncker,elise,appl. math & progr. div - k.u.leuven
289 ***purpose the routine calculates an approximation result to a given
290  integral i = integral of f over (bound,+infinity)
291  or i = integral of f over (-infinity,bound)
292  or i = integral of f over (-infinity,+infinity),
293  hopefully satisfying following claim for accuracy
294  abs(i-result) <= max(epsabs,epsrel*abs(i))
295 ***description
296 
297 integration over infinite intervals
298 standard fortran subroutine
299 
300  f - double precision
301  function subprogram defining the integrand
302  function f(x). the actual name for f needs to be
303  declared e x t e r n a l in the driver program.
304 
305  bound - double precision
306  finite bound of integration range
307  (has no meaning if interval is doubly-infinite)
308 
309  inf - double precision
310  indicating the kind of integration range involved
311  inf = 1 corresponds to (bound,+infinity),
312  inf = -1 to (-infinity,bound),
313  inf = 2 to (-infinity,+infinity).
314 
315  epsabs - double precision
316  absolute accuracy requested
317  epsrel - double precision
318  relative accuracy requested
319  if epsabs <= 0
320  and epsrel < max(50*rel.mach.acc.,0.5d-28),
321  the routine will end with ier = 6.
322 
323  limit - int
324  gives an upper bound on the number of subintervals
325  in the partition of (a,b), limit >= 1
326 
327  on return
328  result - double precision
329  approximation to the integral
330 
331  abserr - double precision
332  estimate of the modulus of the absolute error,
333  which should equal or exceed abs(i-result)
334 
335  neval - int
336  number of integrand evaluations
337 
338  ier - int
339  ier = 0 normal and reliable termination of the
340  routine. it is assumed that the requested
341  accuracy has been achieved.
342  - ier > 0 abnormal termination of the routine. the
343  estimates for result and error are less
344  reliable. it is assumed that the requested
345  accuracy has not been achieved.
346  error messages
347  ier = 1 maximum number of subdivisions allowed
348  has been achieved. one can allow more
349  subdivisions by increasing the value of
350  limit (and taking the according dimension
351  adjustments into account). however,if
352  this yields no improvement it is advised
353  to analyze the integrand in order to
354  determine the integration difficulties.
355  if the position of a local difficulty can
356  be determined (e.g. singularity,
357  discontinuity within the interval) one
358  will probably gain from splitting up the
359  interval at this point and calling the
360  integrator on the subranges. if possible,
361  an appropriate special-purpose integrator
362  should be used, which is designed for
363  handling the type of difficulty involved.
364  = 2 the occurrence of roundoff error is
365  detected, which prevents the requested
366  tolerance from being achieved.
367  the error may be under-estimated.
368  = 3 extremely bad integrand behaviour occurs
369  at some points of the integration
370  interval.
371  = 4 the algorithm does not converge.
372  roundoff error is detected in the
373  extrapolation table.
374  it is assumed that the requested tolerance
375  cannot be achieved, and that the returned
376  result is the best which can be obtained.
377  = 5 the integral is probably divergent, or
378  slowly convergent. it must be noted that
379  divergence can occur with any other value
380  of ier.
381  = 6 the input is invalid, because
382  (epsabs <= 0 and
383  epsrel < max(50*rel.mach.acc.,0.5d-28),
384  result, abserr, neval, last, rlist(1),
385  elist(1) and iord(1) are set to zero.
386  alist(1) and blist(1) are set to 0
387  and 1 respectively.
388 
389  alist - double precision
390  vector of dimension at least limit, the first
391  last elements of which are the left
392  end points of the subintervals in the partition
393  of the transformed integration range (0,1).
394 
395  blist - double precision
396  vector of dimension at least limit, the first
397  last elements of which are the right
398  end points of the subintervals in the partition
399  of the transformed integration range (0,1).
400 
401  rlist - double precision
402  vector of dimension at least limit, the first
403  last elements of which are the integral
404  approximations on the subintervals
405 
406  elist - double precision
407  vector of dimension at least limit, the first
408  last elements of which are the moduli of the
409  absolute error estimates on the subintervals
410 
411  iord - int
412  vector of dimension limit, the first k
413  elements of which are pointers to the
414  error estimates over the subintervals,
415  such that elist(iord(1)), ..., elist(iord(k))
416  form a decreasing sequence, with k = last
417  if last <= (limit/2+2), and k = limit+1-last
418  otherwise
419 
420  last - int
421  number of subintervals actually produced
422  in the subdivision process
423 
424 ***routines called dqelg,dqk15i,dqpsrt
425 ***end prologue dqagie
426 
427 
428  the dimension of rlist2 is determined by the value of
429  limexp in subroutine dqelg.
430 
431  list of major variables
432  -----------------------
433 
434  alist - list of left end points of all subintervals
435  considered up to now
436  blist - list of right end points of all subintervals
437  considered up to now
438  rlist(i) - approximation to the integral over
439  (alist(i),blist(i))
440  rlist2 - array of dimension at least (limexp+2),
441  containing the part of the epsilon table
442  wich is still needed for further computations
443  elist(i) - error estimate applying to rlist(i)
444  maxerr - pointer to the interval with largest error
445  estimate
446  errmax - elist(maxerr)
447  erlast - error on the interval currently subdivided
448  (before that subdivision has taken place)
449  area - sum of the integrals over the subintervals
450  errsum - sum of the errors over the subintervals
451  errbnd - requested accuracy max(epsabs,epsrel*
452  abs(result))
453  *****1 - variable for the left subinterval
454  *****2 - variable for the right subinterval
455  last - index for subdivision
456  nres - number of calls to the extrapolation routine
457  numrl2 - number of elements currently in rlist2. if an
458  appropriate approximation to the compounded
459  integral has been obtained, it is put in
460  rlist2(numrl2) after numrl2 has been increased
461  by one.
462  small - length of the smallest interval considered up
463  to now, multiplied by 1.5
464  erlarg - sum of the errors over the intervals larger
465  than the smallest interval considered up to now
466  extrap - logical variable denoting that the routine
467  is attempting to perform extrapolation. i.e.
468  before subdividing the smallest interval we
469  try to decrease the value of erlarg.
470  noext - logical variable denoting that extrapolation
471  is no longer allowed (true-value)
472 
473  machine dependent constants
474  ---------------------------
475 
476  epmach is the largest relative spacing.
477  uflow is the smallest positive magnitude.
478  oflow is the largest positive magnitude. */
479 
480 /* ***first executable statement dqagie */
481  /* Parameter adjustments */
482  --iord;
483  --elist;
484  --rlist;
485  --blist;
486  --alist;
487 
488  /* Function Body */
489  epmach = DBL_EPSILON;
490 
491 /* test on validity of parameters */
492 /* ----------------------------- */
493 
494  *ier = 0;
495  *neval = 0;
496  *last = 0;
497  *result = 0.;
498  *abserr = 0.;
499  alist[1] = 0.;
500  blist[1] = 1.;
501  rlist[1] = 0.;
502  elist[1] = 0.;
503  iord[1] = 0;
504  if (*epsabs <= 0. && (*epsrel < fmax2(epmach * 50., 5e-29))) *ier = 6;
505  if (*ier == 6) return;
506 
507 /* first approximation to the integral */
508 /* ----------------------------------- */
509 
510 /* determine the interval to be mapped onto (0,1).
511  if inf = 2 the integral is computed as i = i1+i2, where
512  i1 = integral of f over (-infinity,0),
513  i2 = integral of f over (0,+infinity). */
514 
515  boun = *bound;
516  if (*inf == 2) {
517  boun = 0.;
518  }
519 
520  /* Table of constant values */
521  static Float c_b6 = 0.;
522  static Float c_b7 = 1.;
523 
524  rdqk15i(f, ex, &boun, inf, &c_b6, &c_b7, result, abserr, &defabs, &resabs);
525 
526 /* test on accuracy */
527 
528  *last = 1;
529  rlist[1] = *result;
530  elist[1] = *abserr;
531  iord[1] = 1;
532  dres = fabs(*result);
533  errbnd = fmax2(*epsabs, *epsrel * dres);
534  if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) *ier = 2;
535  if (*limit == 1) *ier = 1;
536  if (*ier != 0 || (*abserr <= errbnd && *abserr != resabs)
537  || *abserr == 0.) goto L130;
538 
539 /* initialization */
540 /* -------------- */
541 
542  uflow = DBL_MIN;
543  oflow = DBL_MAX;
544  rlist2[0] = *result;
545  errmax = *abserr;
546  maxerr = 1;
547  area = *result;
548  errsum = *abserr;
549  *abserr = oflow;
550  nrmax = 1;
551  nres = 0;
552  ktmin = 0;
553  numrl2 = 2;
554  extrap = FALSE;
555  noext = FALSE;
556  ierro = 0;
557  iroff1 = 0;
558  iroff2 = 0;
559  iroff3 = 0;
560  ksgn = -1;
561  if (dres >= (1. - epmach * 50.) * defabs) {
562  ksgn = 1;
563  }
564 
565 /* main do-loop */
566 /* ------------ */
567 
568  for (*last = 2; *last <= *limit; ++(*last)) {
569 
570 /* bisect the subinterval with nrmax-th largest error estimate. */
571 
572  a1 = alist[maxerr];
573  b1 = (alist[maxerr] + blist[maxerr]) * .5;
574  a2 = b1;
575  b2 = blist[maxerr];
576  erlast = errmax;
577  rdqk15i(f, ex, &boun, inf, &a1, &b1, &area1, &error1, &resabs, &defab1);
578  rdqk15i(f, ex, &boun, inf, &a2, &b2, &area2, &error2, &resabs, &defab2);
579 
580 /* improve previous approximations to integral
581  and error and test for accuracy. */
582 
583  area12 = area1 + area2;
584  erro12 = error1 + error2;
585  errsum = errsum + erro12 - errmax;
586  area = area + area12 - rlist[maxerr];
587  if (!(defab1 == error1 || defab2 == error2)) {
588  if (fabs(rlist[maxerr] - area12) <= fabs(area12) * 1e-5 &&
589  erro12 >= errmax * .99) {
590  if (extrap)
591  ++iroff2;
592  else /* if (! extrap) */
593  ++iroff1;
594  }
595  if (*last > 10 && erro12 > errmax)
596  ++iroff3;
597  }
598 
599  rlist[maxerr] = area1;
600  rlist[*last] = area2;
601  errbnd = fmax2(*epsabs, *epsrel * fabs(area));
602 
603 /* test for roundoff error and eventually set error flag. */
604 
605  if (iroff1 + iroff2 >= 10 || iroff3 >= 20)
606  *ier = 2;
607  if (iroff2 >= 5)
608  ierro = 3;
609 
610 /* set error flag in the case that the number of
611  subintervals equals limit. */
612 
613  if (*last == *limit)
614  *ier = 1;
615 
616 /* set error flag in the case of bad integrand behaviour
617  at some points of the integration range. */
618 
619  if (fmax2(fabs(a1), fabs(b2)) <=
620  (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
621  {
622  *ier = 4;
623  }
624 
625 /* append the newly-created intervals to the list. */
626 
627  if (error2 <= error1) {
628  alist[*last] = a2;
629  blist[maxerr] = b1;
630  blist[*last] = b2;
631  elist[maxerr] = error1;
632  elist[*last] = error2;
633  }
634  else {
635  alist[maxerr] = a2;
636  alist[*last] = a1;
637  blist[*last] = b1;
638  rlist[maxerr] = area2;
639  rlist[*last] = area1;
640  elist[maxerr] = error2;
641  elist[*last] = error1;
642  }
643 
644 /* call subroutine dqpsrt to maintain the descending ordering
645  in the list of error estimates and select the subinterval
646  with nrmax-th largest error estimate (to be bisected next). */
647 
648  rdqpsrt(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
649  if (errsum <= errbnd) {
650  goto L115;
651  }
652  if (*ier != 0) break;
653  if (*last == 2) { /* L80: */
654  small = .375;
655  erlarg = errsum;
656  ertest = errbnd;
657  rlist2[1] = area; continue;
658  }
659  if (noext) continue;
660 
661  erlarg -= erlast;
662  if (fabs(b1 - a1) > small) {
663  erlarg += erro12;
664  }
665  if (!extrap) {
666 
667 /* test whether the interval to be bisected next is the
668  smallest interval. */
669 
670  if (fabs(blist[maxerr] - alist[maxerr]) > small) {
671  continue;
672  }
673  extrap = TRUE;
674  nrmax = 2;
675  }
676 
677  if (ierro != 3 && erlarg > ertest) {
678 
679 /* the smallest interval has the largest error.
680  before bisecting decrease the sum of the errors over the
681  larger intervals (erlarg) and perform extrapolation. */
682 
683  id = nrmax;
684  jupbnd = *last;
685  if (*last > *limit / 2 + 2) {
686  jupbnd = *limit + 3 - *last;
687  }
688  for (k = id; k <= jupbnd; ++k) {
689  maxerr = iord[nrmax];
690  errmax = elist[maxerr];
691  if (fabs(blist[maxerr] - alist[maxerr]) > small) {
692  goto L90;
693  }
694  ++nrmax;
695  /* L50: */
696  }
697  }
698 /* perform extrapolation. L60: */
699  ++numrl2;
700  rlist2[numrl2 - 1] = area;
701  rdqelg(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
702  ++ktmin;
703  if (ktmin > 5 && *abserr < errsum * .001) {
704  *ier = 5;
705  }
706  if (abseps >= *abserr) {
707  goto L70;
708  }
709  ktmin = 0;
710  *abserr = abseps;
711  *result = reseps;
712  correc = erlarg;
713  ertest = fmax2(*epsabs, *epsrel * fabs(reseps));
714  if (*abserr <= ertest) {
715  break;
716  }
717 
718 /* prepare bisection of the smallest interval. */
719 
720 L70:
721  if (numrl2 == 1) {
722  noext = TRUE;
723  }
724  if (*ier == 5) {
725  break;
726  }
727  maxerr = iord[1];
728  errmax = elist[maxerr];
729  nrmax = 1;
730  extrap = FALSE;
731  small *= .5;
732  erlarg = errsum;
733  L90:
734  ;
735  }
736 
737 /* L100: set final result and error estimate. */
738 /* ------------------------------------ */
739 
740  if (*abserr == oflow) {
741  goto L115;
742  }
743  if (*ier + ierro == 0) {
744  goto L110;
745  }
746  if (ierro == 3) {
747  *abserr += correc;
748  }
749  if (*ier == 0) {
750  *ier = 3;
751  }
752  if (*result == 0. || area == 0.) {
753  if (*abserr > errsum)
754  goto L115;
755 
756  if (area == 0.)
757  goto L130;
758  }
759  else { /* L105: */
760  if (*abserr / fabs(*result) > errsum / fabs(area)) {
761  goto L115;
762  }
763  }
764 
765 /* test on divergence */
766 L110:
767  if (ksgn == -1 && fmax2(fabs(*result), fabs(area)) <= defabs * .01) {
768  goto L130;
769  }
770  if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
771  *ier = 6;
772  }
773  goto L130;
774 
775 /* compute global integral sum. */
776 
777 L115:
778  *result = 0.;
779  for (k = 1; k <= *last; ++k)
780  *result += rlist[k];
781 
782  *abserr = errsum;
783 L130:
784  *neval = *last * 30 - 15;
785  if (*inf == 2) {
786  *neval <<= 1;
787  }
788  if (*ier > 2) {
789  --(*ier);
790  }
791  return;
792 } /* rdqagie_ */
793 
794 template<class Float, class integr_fn>
795 void Rdqags(integr_fn f, void *ex, Float *a, Float *b,
796  Float *epsabs, Float *epsrel,
797  Float *result, Float *abserr, int *neval, int *ier,
798  int *limit, int *lenw, int *last, int *iwork, Float *work)
799 {
800  int l1, l2, l3;
801 
802 /*
803 ***begin prologue dqags
804 ***date written 800101 (yymmdd)
805 ***revision date 830518 (yymmdd)
806 ***category no. h2a1a1
807 ***keywords automatic integrator, general-purpose,
808  (end-point) singularities, extrapolation,
809  globally adaptive
810 ***author piessens,robert,appl. math. & progr. div. - k.u.leuven
811  de doncker,elise,appl. math. & prog. div. - k.u.leuven
812 ***purpose the routine calculates an approximation result to a given
813  definite integral i = integral of f over (a,b),
814  hopefully satisfying following claim for accuracy
815  abs(i-result) <= max(epsabs,epsrel*abs(i)).
816 ***description
817 
818  computation of a definite integral
819  standard fortran subroutine
820  double precision version
821 
822 
823  parameters
824  on entry
825  f - double precision
826  function subprogram defining the integrand
827  function f(x). the actual name for f needs to be
828  declared e x t e r n a l in the driver program.
829 
830  a - double precision
831  lower limit of integration
832 
833  b - double precision
834  upper limit of integration
835 
836  epsabs - double precision
837  absolute accuracy requested
838  epsrel - double precision
839  relative accuracy requested
840  if epsabs <= 0
841  and epsrel < max(50*rel.mach.acc.,0.5d-28),
842  the routine will end with ier = 6.
843 
844  on return
845  result - double precision
846  approximation to the integral
847 
848  abserr - double precision
849  estimate of the modulus of the absolute error,
850  which should equal or exceed abs(i-result)
851 
852  neval - int
853  number of integrand evaluations
854 
855  ier - int
856  ier = 0 normal and reliable termination of the
857  routine. it is assumed that the requested
858  accuracy has been achieved.
859  ier > 0 abnormal termination of the routine
860  the estimates for integral and error are
861  less reliable. it is assumed that the
862  requested accuracy has not been achieved.
863  error messages
864  ier = 1 maximum number of subdivisions allowed
865  has been achieved. one can allow more sub-
866  divisions by increasing the value of limit
867  (and taking the according dimension
868  adjustments into account. however, if
869  this yields no improvement it is advised
870  to analyze the integrand in order to
871  determine the integration difficulties. if
872  the position of a local difficulty can be
873  determined (e.g. singularity,
874  discontinuity within the interval) one
875  will probably gain from splitting up the
876  interval at this point and calling the
877  integrator on the subranges. if possible,
878  an appropriate special-purpose integrator
879  should be used, which is designed for
880  handling the type of difficulty involved.
881  = 2 the occurrence of roundoff error is detec-
882  ted, which prevents the requested
883  tolerance from being achieved.
884  the error may be under-estimated.
885  = 3 extremely bad integrand behaviour
886  occurs at some points of the integration
887  interval.
888  = 4 the algorithm does not converge.
889  roundoff error is detected in the
890  extrapolation table. it is presumed that
891  the requested tolerance cannot be
892  achieved, and that the returned result is
893  the best which can be obtained.
894  = 5 the integral is probably divergent, or
895  slowly convergent. it must be noted that
896  divergence can occur with any other value
897  of ier.
898  = 6 the input is invalid, because
899  (epsabs <= 0 and
900  epsrel < max(50*rel.mach.acc.,0.5d-28)
901  or limit < 1 or lenw < limit*4.
902  result, abserr, neval, last are set to
903  zero.except when limit or lenw is invalid,
904  iwork(1), work(limit*2+1) and
905  work(limit*3+1) are set to zero, work(1)
906  is set to a and work(limit+1) to b.
907 
908  dimensioning parameters
909  limit - int
910  dimensioning parameter for iwork
911  limit determines the maximum number of subintervals
912  in the partition of the given integration interval
913  (a,b), limit >= 1.
914  if limit < 1, the routine will end with ier = 6.
915 
916  lenw - int
917  dimensioning parameter for work
918  lenw must be at least limit*4.
919  if lenw < limit*4, the routine will end
920  with ier = 6.
921 
922  last - int
923  on return, last equals the number of subintervals
924  produced in the subdivision process, detemines the
925  number of significant elements actually in the work
926  arrays.
927 
928  work arrays
929  iwork - int
930  vector of dimension at least limit, the first k
931  elements of which contain pointers
932  to the error estimates over the subintervals
933  such that work(limit*3+iwork(1)),... ,
934  work(limit*3+iwork(k)) form a decreasing
935  sequence, with k = last if last <= (limit/2+2),
936  and k = limit+1-last otherwise
937 
938  work - double precision
939  vector of dimension at least lenw
940  on return
941  work(1), ..., work(last) contain the left
942  end-points of the subintervals in the
943  partition of (a,b),
944  work(limit+1), ..., work(limit+last) contain
945  the right end-points,
946  work(limit*2+1), ..., work(limit*2+last) contain
947  the integral approximations over the subintervals,
948  work(limit*3+1), ..., work(limit*3+last)
949  contain the error estimates.
950 
951 ***routines called dqagse
952 ***end prologue dqags */
953 
954 /* check validity of limit and lenw. */
955 
956  *ier = 6;
957  *neval = 0;
958  *last = 0;
959  *result = 0.;
960  *abserr = 0.;
961  if (*limit < 1 || *lenw < *limit *4) return;
962 
963 /* prepare call for dqagse. */
964 
965  l1 = *limit;
966  l2 = *limit + l1;
967  l3 = *limit + l2;
968 
969  rdqagse(f, ex, a, b, epsabs, epsrel, limit, result, abserr, neval, ier,
970  work, &work[l1], &work[l2], &work[l3], iwork, last);
971 
972  return;
973 } /* rdqags_ */
974 
975 template<class Float, class integr_fn> static
976 void rdqagse(integr_fn f, void *ex, Float *a, Float *b, Float *
977  epsabs, Float *epsrel, int *limit, Float *result,
978  Float *abserr, int *neval, int *ier, Float *alist,
979  Float *blist, Float *rlist, Float *elist, int *
980  iord, int *last)
981 {
982  /* Local variables */
983  bool noext, extrap;
984  int k,ksgn, nres;
985  int ierro;
986  int ktmin, nrmax;
987  int iroff1, iroff2, iroff3;
988  int id;
989  int numrl2;
990  int jupbnd;
991  int maxerr;
992  Float res3la[3];
993  Float rlist2[52];
994  Float abseps, area, area1, area2, area12, dres, epmach;
995  Float a1, a2, b1, b2, defabs, defab1, defab2, oflow, uflow, resabs, reseps;
996  Float error1, error2, erro12, errbnd, erlast, errmax, errsum;
997 
998  Float correc = 0.0, erlarg = 0.0, ertest = 0.0, small = 0.0;
999 /*
1000 ***begin prologue dqagse
1001 ***date written 800101 (yymmdd)
1002 ***revision date 830518 (yymmdd)
1003 ***category no. h2a1a1
1004 ***keywords automatic integrator, general-purpose,
1005  (end point) singularities, extrapolation,
1006  globally adaptive
1007 ***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1008  de doncker,elise,appl. math. & progr. div. - k.u.leuven
1009 ***purpose the routine calculates an approximation result to a given
1010  definite integral i = integral of f over (a,b),
1011  hopefully satisfying following claim for accuracy
1012  abs(i-result) <= max(epsabs,epsrel*abs(i)).
1013 ***description
1014 
1015  computation of a definite integral
1016  standard fortran subroutine
1017  double precision version
1018 
1019  parameters
1020  on entry
1021  f - double precision
1022  function subprogram defining the integrand
1023  function f(x). the actual name for f needs to be
1024  declared e x t e r n a l in the driver program.
1025 
1026  a - double precision
1027  lower limit of integration
1028 
1029  b - double precision
1030  upper limit of integration
1031 
1032  epsabs - double precision
1033  absolute accuracy requested
1034  epsrel - double precision
1035  relative accuracy requested
1036  if epsabs <= 0
1037  and epsrel < max(50*rel.mach.acc.,0.5d-28),
1038  the routine will end with ier = 6.
1039 
1040  limit - int
1041  gives an upperbound on the number of subintervals
1042  in the partition of (a,b)
1043 
1044  on return
1045  result - double precision
1046  approximation to the integral
1047 
1048  abserr - double precision
1049  estimate of the modulus of the absolute error,
1050  which should equal or exceed abs(i-result)
1051 
1052  neval - int
1053  number of integrand evaluations
1054 
1055  ier - int
1056  ier = 0 normal and reliable termination of the
1057  routine. it is assumed that the requested
1058  accuracy has been achieved.
1059  ier > 0 abnormal termination of the routine
1060  the estimates for integral and error are
1061  less reliable. it is assumed that the
1062  requested accuracy has not been achieved.
1063  error messages
1064  = 1 maximum number of subdivisions allowed
1065  has been achieved. one can allow more sub-
1066  divisions by increasing the value of limit
1067  (and taking the according dimension
1068  adjustments into account). however, if
1069  this yields no improvement it is advised
1070  to analyze the integrand in order to
1071  determine the integration difficulties. if
1072  the position of a local difficulty can be
1073  determined (e.g. singularity,
1074  discontinuity within the interval) one
1075  will probably gain from splitting up the
1076  interval at this point and calling the
1077  integrator on the subranges. if possible,
1078  an appropriate special-purpose integrator
1079  should be used, which is designed for
1080  handling the type of difficulty involved.
1081  = 2 the occurrence of roundoff error is detec-
1082  ted, which prevents the requested
1083  tolerance from being achieved.
1084  the error may be under-estimated.
1085  = 3 extremely bad integrand behaviour
1086  occurs at some points of the integration
1087  interval.
1088  = 4 the algorithm does not converge.
1089  roundoff error is detected in the
1090  extrapolation table.
1091  it is presumed that the requested
1092  tolerance cannot be achieved, and that the
1093  returned result is the best which can be
1094  obtained.
1095  = 5 the integral is probably divergent, or
1096  slowly convergent. it must be noted that
1097  divergence can occur with any other value
1098  of ier.
1099  = 6 the input is invalid, because
1100  epsabs <= 0 and
1101  epsrel < max(50*rel.mach.acc.,0.5d-28).
1102  result, abserr, neval, last, rlist(1),
1103  iord(1) and elist(1) are set to zero.
1104  alist(1) and blist(1) are set to a and b
1105  respectively.
1106 
1107  alist - double precision
1108  vector of dimension at least limit, the first
1109  last elements of which are the left end points
1110  of the subintervals in the partition of the
1111  given integration range (a,b)
1112 
1113  blist - double precision
1114  vector of dimension at least limit, the first
1115  last elements of which are the right end points
1116  of the subintervals in the partition of the given
1117  integration range (a,b)
1118 
1119  rlist - double precision
1120  vector of dimension at least limit, the first
1121  last elements of which are the integral
1122  approximations on the subintervals
1123 
1124  elist - double precision
1125  vector of dimension at least limit, the first
1126  last elements of which are the moduli of the
1127  absolute error estimates on the subintervals
1128 
1129  iord - int
1130  vector of dimension at least limit, the first k
1131  elements of which are pointers to the
1132  error estimates over the subintervals,
1133  such that elist(iord(1)), ..., elist(iord(k))
1134  form a decreasing sequence, with k = last
1135  if last <= (limit/2+2), and k = limit+1-last
1136  otherwise
1137 
1138  last - int
1139  number of subintervals actually produced in the
1140  subdivision process
1141 
1142 ***references (none)
1143 ***routines called dqelg,dqk21,dqpsrt
1144 ***end prologue dqagse
1145 
1146 
1147 
1148  the dimension of rlist2 is determined by the value of
1149  limexp in subroutine dqelg (rlist2 should be of dimension
1150  (limexp+2) at least).
1151 
1152  list of major variables
1153  -----------------------
1154 
1155  alist - list of left end points of all subintervals
1156  considered up to now
1157  blist - list of right end points of all subintervals
1158  considered up to now
1159  rlist(i) - approximation to the integral over
1160  (alist(i),blist(i))
1161  rlist2 - array of dimension at least limexp+2 containing
1162  the part of the epsilon table which is still
1163  needed for further computations
1164  elist(i) - error estimate applying to rlist(i)
1165  maxerr - pointer to the interval with largest error
1166  estimate
1167  errmax - elist(maxerr)
1168  erlast - error on the interval currently subdivided
1169  (before that subdivision has taken place)
1170  area - sum of the integrals over the subintervals
1171  errsum - sum of the errors over the subintervals
1172  errbnd - requested accuracy max(epsabs,epsrel*
1173  abs(result))
1174  *****1 - variable for the left interval
1175  *****2 - variable for the right interval
1176  last - index for subdivision
1177  nres - number of calls to the extrapolation routine
1178  numrl2 - number of elements currently in rlist2. if an
1179  appropriate approximation to the compounded
1180  integral has been obtained it is put in
1181  rlist2(numrl2) after numrl2 has been increased
1182  by one.
1183  small - length of the smallest interval considered up
1184  to now, multiplied by 1.5
1185  erlarg - sum of the errors over the intervals larger
1186  than the smallest interval considered up to now
1187  extrap - logical variable denoting that the routine is
1188  attempting to perform extrapolation i.e. before
1189  subdividing the smallest interval we try to
1190  decrease the value of erlarg.
1191  noext - logical variable denoting that extrapolation
1192  is no longer allowed (true value)
1193 
1194  machine dependent constants
1195  ---------------------------
1196 
1197  epmach is the largest relative spacing.
1198  uflow is the smallest positive magnitude.
1199  oflow is the largest positive magnitude. */
1200 
1201 /* ***first executable statement dqagse */
1202  /* Parameter adjustments */
1203  --iord;
1204  --elist;
1205  --rlist;
1206  --blist;
1207  --alist;
1208 
1209  /* Function Body */
1210  epmach = DBL_EPSILON;
1211 
1212 /* test on validity of parameters */
1213 /* ------------------------------ */
1214  *ier = 0;
1215  *neval = 0;
1216  *last = 0;
1217  *result = 0.;
1218  *abserr = 0.;
1219  alist[1] = *a;
1220  blist[1] = *b;
1221  rlist[1] = 0.;
1222  elist[1] = 0.;
1223  if (*epsabs <= 0. && *epsrel < fmax2(epmach * 50., 5e-29)) {
1224  *ier = 6;
1225  return;
1226  }
1227 
1228 /* first approximation to the integral */
1229 /* ----------------------------------- */
1230 
1231  uflow = DBL_MIN;
1232  oflow = DBL_MAX;
1233  ierro = 0;
1234  rdqk21(f, ex, a, b, result, abserr, &defabs, &resabs);
1235 
1236 /* test on accuracy. */
1237 
1238  dres = fabs(*result);
1239  errbnd = fmax2(*epsabs, *epsrel * dres);
1240  *last = 1;
1241  rlist[1] = *result;
1242  elist[1] = *abserr;
1243  iord[1] = 1;
1244  if (*abserr <= epmach * 100. * defabs && *abserr > errbnd)
1245  *ier = 2;
1246  if (*limit == 1)
1247  *ier = 1;
1248  if (*ier != 0 || (*abserr <= errbnd && *abserr != resabs)
1249  || *abserr == 0.) goto L140;
1250 
1251 /* initialization */
1252 /* -------------- */
1253 
1254  rlist2[0] = *result;
1255  errmax = *abserr;
1256  maxerr = 1;
1257  area = *result;
1258  errsum = *abserr;
1259  *abserr = oflow;
1260  nrmax = 1;
1261  nres = 0;
1262  numrl2 = 2;
1263  ktmin = 0;
1264  extrap = FALSE;
1265  noext = FALSE;
1266  iroff1 = 0;
1267  iroff2 = 0;
1268  iroff3 = 0;
1269  ksgn = -1;
1270  if (dres >= (1. - epmach * 50.) * defabs) {
1271  ksgn = 1;
1272  }
1273 
1274 /* main do-loop */
1275 /* ------------ */
1276 
1277  for (*last = 2; *last <= *limit; ++(*last)) {
1278 
1279 /* bisect the subinterval with the nrmax-th largest error estimate. */
1280 
1281  a1 = alist[maxerr];
1282  b1 = (alist[maxerr] + blist[maxerr]) * .5;
1283  a2 = b1;
1284  b2 = blist[maxerr];
1285  erlast = errmax;
1286  rdqk21(f, ex, &a1, &b1, &area1, &error1, &resabs, &defab1);
1287  rdqk21(f, ex, &a2, &b2, &area2, &error2, &resabs, &defab2);
1288 
1289 /* improve previous approximations to integral
1290  and error and test for accuracy. */
1291 
1292  area12 = area1 + area2;
1293  erro12 = error1 + error2;
1294  errsum = errsum + erro12 - errmax;
1295  area = area + area12 - rlist[maxerr];
1296  if (!(defab1 == error1 || defab2 == error2)) {
1297 
1298  if (fabs(rlist[maxerr] - area12) <= fabs(area12) * 1e-5 &&
1299  erro12 >= errmax * .99) {
1300  if (extrap)
1301  ++iroff2;
1302  else /* if(! extrap) */
1303  ++iroff1;
1304  }
1305  if (*last > 10 && erro12 > errmax)
1306  ++iroff3;
1307  }
1308  rlist[maxerr] = area1;
1309  rlist[*last] = area2;
1310  errbnd = fmax2(*epsabs, *epsrel * fabs(area));
1311 
1312 /* test for roundoff error and eventually set error flag. */
1313 
1314  if (iroff1 + iroff2 >= 10 || iroff3 >= 20)
1315  *ier = 2;
1316  if (iroff2 >= 5)
1317  ierro = 3;
1318 
1319 /* set error flag in the case that the number of subintervals equals limit. */
1320  if (*last == *limit)
1321  *ier = 1;
1322 
1323 /* set error flag in the case of bad integrand behaviour
1324  at a point of the integration range. */
1325 
1326  if (fmax2(fabs(a1), fabs(b2)) <=
1327  (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3)) {
1328  *ier = 4;
1329  }
1330 
1331 /* append the newly-created intervals to the list. */
1332 
1333  if (error2 > error1) {
1334  alist[maxerr] = a2;
1335  alist[*last] = a1;
1336  blist[*last] = b1;
1337  rlist[maxerr] = area2;
1338  rlist[*last] = area1;
1339  elist[maxerr] = error2;
1340  elist[*last] = error1;
1341  } else {
1342  alist[*last] = a2;
1343  blist[maxerr] = b1;
1344  blist[*last] = b2;
1345  elist[maxerr] = error1;
1346  elist[*last] = error2;
1347  }
1348 
1349 /* call subroutine dqpsrt to maintain the descending ordering
1350  in the list of error estimates and select the subinterval
1351  with nrmax-th largest error estimate (to be bisected next). */
1352 
1353 /*L30:*/
1354  rdqpsrt(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
1355 
1356  if (errsum <= errbnd) goto L115;/* ***jump out of do-loop */
1357  if (*ier != 0) break;
1358  if (*last == 2) { /* L80: */
1359  small = fabs(*b - *a) * .375;
1360  erlarg = errsum;
1361  ertest = errbnd;
1362  rlist2[1] = area; continue;
1363  }
1364  if (noext) continue;
1365 
1366  erlarg -= erlast;
1367  if (fabs(b1 - a1) > small) {
1368  erlarg += erro12;
1369  }
1370  if (!extrap) {
1371 
1372 /* test whether the interval to be bisected next is the
1373  smallest interval. */
1374 
1375  if (fabs(blist[maxerr] - alist[maxerr]) > small) {
1376  continue;
1377  }
1378  extrap = TRUE;
1379  nrmax = 2;
1380  }
1381 
1382  if (ierro != 3 && erlarg > ertest) {
1383 
1384 /* the smallest interval has the largest error.
1385  before bisecting decrease the sum of the errors over the
1386  larger intervals (erlarg) and perform extrapolation. */
1387 
1388  id = nrmax;
1389  jupbnd = *last;
1390  if (*last > *limit / 2 + 2) {
1391  jupbnd = *limit + 3 - *last;
1392  }
1393  for (k = id; k <= jupbnd; ++k) {
1394  maxerr = iord[nrmax];
1395  errmax = elist[maxerr];
1396  if (fabs(blist[maxerr] - alist[maxerr]) > small) {
1397  goto L90;
1398  }
1399  ++nrmax;
1400  /* L50: */
1401  }
1402  }
1403 /* perform extrapolation. L60: */
1404 
1405  ++numrl2;
1406  rlist2[numrl2 - 1] = area;
1407  rdqelg(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
1408  ++ktmin;
1409  if (ktmin > 5 && *abserr < errsum * .001) {
1410  *ier = 5;
1411  }
1412  if (abseps < *abserr) {
1413  ktmin = 0;
1414  *abserr = abseps;
1415  *result = reseps;
1416  correc = erlarg;
1417  ertest = fmax2(*epsabs, *epsrel * fabs(reseps));
1418  if (*abserr <= ertest) {
1419  break;
1420  }
1421  }
1422 
1423 /* prepare bisection of the smallest interval. L70: */
1424 
1425  if (numrl2 == 1) {
1426  noext = TRUE;
1427  }
1428  if (*ier == 5) {
1429  break;
1430  }
1431  maxerr = iord[1];
1432  errmax = elist[maxerr];
1433  nrmax = 1;
1434  extrap = FALSE;
1435  small *= .5;
1436  erlarg = errsum;
1437  L90:
1438  ;
1439  }
1440 
1441 
1442 /* L100: set final result and error estimate. */
1443 /* ------------------------------------ */
1444 
1445  if (*abserr == oflow) goto L115;
1446  if (*ier + ierro == 0) goto L110;
1447  if (ierro == 3)
1448  *abserr += correc;
1449  if (*ier == 0)
1450  *ier = 3;
1451  if (*result == 0. || area == 0.) {
1452  if (*abserr > errsum) goto L115;
1453  if (area == 0.) goto L130;
1454  }
1455  else { /* L105:*/
1456  if (*abserr / fabs(*result) > errsum / fabs(area))
1457  goto L115;
1458  }
1459 
1460 L110:/* test on divergence. */
1461  if (ksgn == -1 && fmax2(fabs(*result), fabs(area)) <= defabs * .01) {
1462  goto L130;
1463  }
1464  if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
1465  *ier = 5;
1466  }
1467  goto L130;
1468 
1469 L115:/* compute global integral sum. */
1470  *result = 0.;
1471  for (k = 1; k <= *last; ++k)
1472  *result += rlist[k];
1473  *abserr = errsum;
1474 L130:
1475  if (*ier > 2)
1476 L140:
1477  *neval = *last * 42 - 21;
1478  return;
1479 } /* rdqagse_ */
1480 
1481 
1482 template<class Float, class integr_fn> static void rdqk15i(integr_fn f, void *ex,
1483  Float *boun, int *inf, Float *a, Float *b,
1484  Float *result,
1485  Float *abserr, Float *resabs, Float *resasc)
1486 {
1487  /* Initialized data */
1488 
1489  static double wg[8] = {
1490  0., .129484966168869693270611432679082,
1491  0., .27970539148927666790146777142378,
1492  0., .381830050505118944950369775488975,
1493  0., .417959183673469387755102040816327 };
1494  static double xgk[8] = {
1495  .991455371120812639206854697526329,
1496  .949107912342758524526189684047851,
1497  .864864423359769072789712788640926,
1498  .741531185599394439863864773280788,
1499  .58608723546769113029414483825873,
1500  .405845151377397166906606412076961,
1501  .207784955007898467600689403773245, 0. };
1502  static double wgk[8] = {
1503  .02293532201052922496373200805897,
1504  .063092092629978553290700663189204,
1505  .104790010322250183839876322541518,
1506  .140653259715525918745189590510238,
1507  .16900472663926790282658342659855,
1508  .190350578064785409913256402421014,
1509  .204432940075298892414161999234649,
1510  .209482141084727828012999174891714 };
1511 
1512  /* Local variables */
1513  Float absc, dinf, resg, resk, fsum, absc1, absc2, fval1, fval2;
1514  int j;
1515  Float hlgth, centr, reskh, uflow;
1516  Float tabsc1, tabsc2, fc, epmach;
1517  Float fv1[7], fv2[7], vec[15], vec2[15];
1518 
1519 /*
1520 ***begin prologue dqk15i
1521 ***date written 800101 (yymmdd)
1522 ***revision date 830518 (yymmdd)
1523 ***category no. h2a3a2,h2a4a2
1524 ***keywords 15-point transformed gauss-kronrod rules
1525 ***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1526  de doncker,elise,appl. math. & progr. div. - k.u.leuven
1527 ***purpose the original (infinite integration range is mapped
1528  onto the interval (0,1) and (a,b) is a part of (0,1).
1529  it is the purpose to compute
1530  i = integral of transformed integrand over (a,b),
1531  j = integral of abs(transformed integrand) over (a,b).
1532 ***description
1533 
1534  integration rule
1535  standard fortran subroutine
1536  double precision version
1537 
1538  parameters
1539  on entry
1540  f - double precision
1541  fuction subprogram defining the integrand
1542  function f(x). the actual name for f needs to be
1543  declared e x t e r n a l in the calling program.
1544 
1545  boun - double precision
1546  finite bound of original integration
1547  range (set to zero if inf = +2)
1548 
1549  inf - int
1550  if inf = -1, the original interval is
1551  (-infinity,bound),
1552  if inf = +1, the original interval is
1553  (bound,+infinity),
1554  if inf = +2, the original interval is
1555  (-infinity,+infinity) and
1556  the integral is computed as the sum of two
1557  integrals, one over (-infinity,0) and one over
1558  (0,+infinity).
1559 
1560  a - double precision
1561  lower limit for integration over subrange
1562  of (0,1)
1563 
1564  b - double precision
1565  upper limit for integration over subrange
1566  of (0,1)
1567 
1568  on return
1569  result - double precision
1570  approximation to the integral i
1571  result is computed by applying the 15-point
1572  kronrod rule(resk) obtained by optimal addition
1573  of abscissae to the 7-point gauss rule(resg).
1574 
1575  abserr - double precision
1576  estimate of the modulus of the absolute error,
1577  which should equal or exceed abs(i-result)
1578 
1579  resabs - double precision
1580  approximation to the integral j
1581 
1582  resasc - double precision
1583  approximation to the integral of
1584  abs((transformed integrand)-i/(b-a)) over (a,b)
1585 
1586 ***references (none)
1587 ***end prologue dqk15i
1588 
1589 
1590  the abscissae and weights are supplied for the interval
1591  (-1,1). because of symmetry only the positive abscissae and
1592  their corresponding weights are given.
1593 
1594  xgk - abscissae of the 15-point kronrod rule
1595  xgk(2), xgk(4), ... abscissae of the 7-point
1596  gauss rule
1597  xgk(1), xgk(3), ... abscissae which are optimally
1598  added to the 7-point gauss rule
1599 
1600  wgk - weights of the 15-point kronrod rule
1601 
1602  wg - weights of the 7-point gauss rule, corresponding
1603  to the abscissae xgk(2), xgk(4), ...
1604  wg(1), wg(3), ... are set to zero.
1605 
1606 
1607 
1608 
1609 
1610  list of major variables
1611  -----------------------
1612 
1613  centr - mid point of the interval
1614  hlgth - half-length of the interval
1615  absc* - abscissa
1616  tabsc* - transformed abscissa
1617  fval* - function value
1618  resg - result of the 7-point gauss formula
1619  resk - result of the 15-point kronrod formula
1620  reskh - approximation to the mean value of the transformed
1621  integrand over (a,b), i.e. to i/(b-a)
1622 
1623  machine dependent constants
1624  ---------------------------
1625 
1626  epmach is the largest relative spacing.
1627  uflow is the smallest positive magnitude.
1628 */
1629 
1630 /* ***first executable statement dqk15i */
1631  epmach = DBL_EPSILON;
1632  uflow = DBL_MIN;
1633  dinf = (double) imin2(1, *inf);
1634 
1635  centr = (*a + *b) * .5;
1636  hlgth = (*b - *a) * .5;
1637  tabsc1 = *boun + dinf * (1. - centr) / centr;
1638  vec[0] = tabsc1;
1639  if (*inf == 2) {
1640  vec2[0] = -tabsc1;
1641  }
1642  for (j = 1; j <= 7; ++j) {
1643  absc = hlgth * xgk[j - 1];
1644  absc1 = centr - absc;
1645  absc2 = centr + absc;
1646  tabsc1 = *boun + dinf * (1. - absc1) / absc1;
1647  tabsc2 = *boun + dinf * (1. - absc2) / absc2;
1648  vec[(j << 1) - 1] = tabsc1;
1649  vec[j * 2] = tabsc2;
1650  if (*inf == 2) {
1651  vec2[(j << 1) - 1] = -tabsc1;
1652  vec2[j * 2] = -tabsc2;
1653  }
1654 /* L5: */
1655  }
1656  f(vec, 15, ex); /* -> new vec[] overwriting old vec[] */
1657  if (*inf == 2) f(vec2, 15, ex);
1658  fval1 = vec[0];
1659  if (*inf == 2) fval1 += vec2[0];
1660  fc = fval1 / centr / centr;
1661 
1662 /* compute the 15-point kronrod approximation to
1663  the integral, and estimate the error. */
1664 
1665  resg = wg[7] * fc;
1666  resk = wgk[7] * fc;
1667  *resabs = fabs(resk);
1668  for (j = 1; j <= 7; ++j) {
1669  absc = hlgth * xgk[j - 1];
1670  absc1 = centr - absc;
1671  absc2 = centr + absc;
1672  tabsc1 = *boun + dinf * (1. - absc1) / absc1;
1673  tabsc2 = *boun + dinf * (1. - absc2) / absc2;
1674  fval1 = vec[(j << 1) - 1];
1675  fval2 = vec[j * 2];
1676  if (*inf == 2) {
1677  fval1 += vec2[(j << 1) - 1];
1678  }
1679  if (*inf == 2) {
1680  fval2 += vec2[j * 2];
1681  }
1682  fval1 = fval1 / absc1 / absc1;
1683  fval2 = fval2 / absc2 / absc2;
1684  fv1[j - 1] = fval1;
1685  fv2[j - 1] = fval2;
1686  fsum = fval1 + fval2;
1687  resg += wg[j - 1] * fsum;
1688  resk += wgk[j - 1] * fsum;
1689  *resabs += wgk[j - 1] * (fabs(fval1) + fabs(fval2));
1690 /* L10: */
1691  }
1692  reskh = resk * .5;
1693  *resasc = wgk[7] * fabs(fc - reskh);
1694  for (j = 1; j <= 7; ++j) {
1695  *resasc += wgk[j - 1] * (fabs(fv1[j - 1] - reskh) +
1696  fabs(fv2[j - 1] - reskh));
1697 /* L20: */
1698  }
1699  *result = resk * hlgth;
1700  *resasc *= hlgth;
1701  *resabs *= hlgth;
1702  *abserr = fabs((resk - resg) * hlgth);
1703  if (*resasc != 0. && *abserr != 0.) {
1704  *abserr = *resasc * fmin2(1., pow(*abserr * 200. / *resasc, 1.5));
1705  }
1706  if (*resabs > uflow / (epmach * 50.)) {
1707  *abserr = fmax2(epmach * 50. * *resabs, *abserr);
1708  }
1709  return;
1710 } /* rdqk15i_ */
1711 
1712 template<class Float> static void rdqelg(int *n, Float *epstab, Float *
1713  result, Float *abserr, Float *res3la, int *nres)
1714 {
1715  /* Local variables */
1716  int i__, indx, ib, ib2, ie, k1, k2, k3, num, newelm, limexp;
1717  Float delta1, delta2, delta3, e0, e1, e1abs, e2, e3, epmach, epsinf;
1718  Float oflow, ss, res;
1719  Float errA, err1, err2, err3, tol1, tol2, tol3;
1720 
1721 /* ***begin prologue dqelg
1722 ***refer to dqagie,dqagoe,dqagpe,dqagse
1723 ***revision date 830518 (yymmdd)
1724 ***keywords epsilon algorithm, convergence acceleration,
1725  extrapolation
1726 ***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1727  de doncker,elise,appl. math & progr. div. - k.u.leuven
1728 ***purpose the routine determines the limit of a given sequence of
1729  approximations, by means of the epsilon algorithm of
1730  p.wynn. an estimate of the absolute error is also given.
1731  the condensed epsilon table is computed. only those
1732  elements needed for the computation of the next diagonal
1733  are preserved.
1734 ***description
1735 
1736  epsilon algorithm
1737  standard fortran subroutine
1738  double precision version
1739 
1740  parameters
1741  n - int
1742  epstab(n) contains the new element in the
1743  first column of the epsilon table.
1744 
1745  epstab - double precision
1746  vector of dimension 52 containing the elements
1747  of the two lower diagonals of the triangular
1748  epsilon table. the elements are numbered
1749  starting at the right-hand corner of the
1750  triangle.
1751 
1752  result - double precision
1753  resulting approximation to the integral
1754 
1755  abserr - double precision
1756  estimate of the absolute error computed from
1757  result and the 3 previous results
1758 
1759  res3la - double precision
1760  vector of dimension 3 containing the last 3
1761  results
1762 
1763  nres - int
1764  number of calls to the routine
1765  (should be zero at first call)
1766 
1767 ***end prologue dqelg
1768 
1769 
1770  list of major variables
1771  -----------------------
1772 
1773  e0 - the 4 elements on which the computation of a new
1774  e1 element in the epsilon table is based
1775  e2
1776  e3 e0
1777  e3 e1 new
1778  e2
1779 
1780  newelm - number of elements to be computed in the new diagonal
1781  errA - errA = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
1782  result - the element in the new diagonal with least value of errA
1783 
1784  machine dependent constants
1785  ---------------------------
1786 
1787  epmach is the largest relative spacing.
1788  oflow is the largest positive magnitude.
1789  limexp is the maximum number of elements the epsilon
1790  table can contain. if this number is reached, the upper
1791  diagonal of the epsilon table is deleted. */
1792 
1793 /* ***first executable statement dqelg */
1794  /* Parameter adjustments */
1795  --res3la;
1796  --epstab;
1797 
1798  /* Function Body */
1799  epmach = DBL_EPSILON;
1800  oflow = DBL_MAX;
1801  ++(*nres);
1802  *abserr = oflow;
1803  *result = epstab[*n];
1804  if (*n < 3) {
1805  goto L100;
1806  }
1807  limexp = 50;
1808  epstab[*n + 2] = epstab[*n];
1809  newelm = (*n - 1) / 2;
1810  epstab[*n] = oflow;
1811  num = *n;
1812  k1 = *n;
1813  for (i__ = 1; i__ <= newelm; ++i__) {
1814  k2 = k1 - 1;
1815  k3 = k1 - 2;
1816  res = epstab[k1 + 2];
1817  e0 = epstab[k3];
1818  e1 = epstab[k2];
1819  e2 = res;
1820  e1abs = fabs(e1);
1821  delta2 = e2 - e1;
1822  err2 = fabs(delta2);
1823  tol2 = fmax2(fabs(e2), e1abs) * epmach;
1824  delta3 = e1 - e0;
1825  err3 = fabs(delta3);
1826  tol3 = fmax2(e1abs, fabs(e0)) * epmach;
1827  if (err2 <= tol2 && err3 <= tol3) {
1828  /* if e0, e1 and e2 are equal to within machine
1829  accuracy, convergence is assumed. */
1830  *result = res;/* result = e2 */
1831  *abserr = err2 + err3;/* abserr = fabs(e1-e0)+fabs(e2-e1) */
1832 
1833  goto L100; /* ***jump out of do-loop */
1834  }
1835 
1836  e3 = epstab[k1];
1837  epstab[k1] = e1;
1838  delta1 = e1 - e3;
1839  err1 = fabs(delta1);
1840  tol1 = fmax2(e1abs, fabs(e3)) * epmach;
1841 
1842 /* if two elements are very close to each other, omit
1843  a part of the table by adjusting the value of n */
1844 
1845  if (err1 > tol1 && err2 > tol2 && err3 > tol3) {
1846  ss = 1. / delta1 + 1. / delta2 - 1. / delta3;
1847  epsinf = fabs(ss * e1);
1848 
1849 /* test to detect irregular behaviour in the table, and
1850  eventually omit a part of the table adjusting the value of n. */
1851 
1852  if (epsinf > 1e-4) {
1853  goto L30;
1854  }
1855  }
1856 
1857  *n = i__ + i__ - 1;
1858  goto L50;/* ***jump out of do-loop */
1859 
1860 
1861 L30:/* compute a new element and eventually adjust the value of result. */
1862 
1863  res = e1 + 1. / ss;
1864  epstab[k1] = res;
1865  k1 += -2;
1866  errA = err2 + fabs(res - e2) + err3;
1867  if (errA <= *abserr) {
1868  *abserr = errA;
1869  *result = res;
1870  }
1871  }
1872 
1873 /* shift the table. */
1874 
1875 L50:
1876  if (*n == limexp) {
1877  *n = (limexp / 2 << 1) - 1;
1878  }
1879 
1880  if (num / 2 << 1 == num) ib = 2; else ib = 1;
1881  ie = newelm + 1;
1882  for (i__ = 1; i__ <= ie; ++i__) {
1883  ib2 = ib + 2;
1884  epstab[ib] = epstab[ib2];
1885  ib = ib2;
1886  }
1887  if (num != *n) {
1888  indx = num - *n + 1;
1889  for (i__ = 1; i__ <= *n; ++i__) {
1890  epstab[i__] = epstab[indx];
1891  ++indx;
1892  }
1893  }
1894  /*L80:*/
1895  if (*nres >= 4) {
1896  /* L90: */
1897  *abserr = fabs(*result - res3la[3]) +
1898  fabs(*result - res3la[2]) +
1899  fabs(*result - res3la[1]);
1900  res3la[1] = res3la[2];
1901  res3la[2] = res3la[3];
1902  res3la[3] = *result;
1903  } else {
1904  res3la[*nres] = *result;
1905  *abserr = oflow;
1906  }
1907 
1908 L100:/* compute error estimate */
1909  *abserr = fmax2(*abserr, epmach * 5. * fabs(*result));
1910  return;
1911 } /* rdqelg_ */
1912 
1913 template<class Float, class integr_fn> static void rdqk21(integr_fn f, void *ex, Float *a, Float *b, Float *result,
1914  Float *abserr, Float *resabs, Float *resasc)
1915 {
1916  /* Initialized data */
1917 
1918  static double wg[5] = { .066671344308688137593568809893332,
1919  .149451349150580593145776339657697,
1920  .219086362515982043995534934228163,
1921  .269266719309996355091226921569469,
1922  .295524224714752870173892994651338 };
1923  static double xgk[11] = { .995657163025808080735527280689003,
1924  .973906528517171720077964012084452,
1925  .930157491355708226001207180059508,
1926  .865063366688984510732096688423493,
1927  .780817726586416897063717578345042,
1928  .679409568299024406234327365114874,
1929  .562757134668604683339000099272694,
1930  .433395394129247190799265943165784,
1931  .294392862701460198131126603103866,
1932  .14887433898163121088482600112972,0. };
1933  static double wgk[11] = { .011694638867371874278064396062192,
1934  .03255816230796472747881897245939,
1935  .05475589657435199603138130024458,
1936  .07503967481091995276704314091619,
1937  .093125454583697605535065465083366,
1938  .109387158802297641899210590325805,
1939  .123491976262065851077958109831074,
1940  .134709217311473325928054001771707,
1941  .142775938577060080797094273138717,
1942  .147739104901338491374841515972068,
1943  .149445554002916905664936468389821 };
1944 
1945 
1946  /* Local variables */
1947  Float fv1[10], fv2[10], vec[21];
1948  Float absc, resg, resk, fsum, fval1, fval2;
1949  Float hlgth, centr, reskh, uflow;
1950  Float fc, epmach, dhlgth;
1951  int j, jtw, jtwm1;
1952 
1953 /* ***begin prologue dqk21
1954 ***date written 800101 (yymmdd)
1955 ***revision date 830518 (yymmdd)
1956 ***category no. h2a1a2
1957 ***keywords 21-point gauss-kronrod rules
1958 ***author piessens,robert,appl. math. & progr. div. - k.u.leuven
1959  de doncker,elise,appl. math. & progr. div. - k.u.leuven
1960 ***purpose to compute i = integral of f over (a,b), with error
1961  estimate
1962  j = integral of abs(f) over (a,b)
1963 ***description
1964 
1965  integration rules
1966  standard fortran subroutine
1967  double precision version
1968 
1969  parameters
1970  on entry
1971  f - double precision
1972  function subprogram defining the integrand
1973  function f(x). the actual name for f needs to be
1974  declared e x t e r n a l in the driver program.
1975 
1976  a - double precision
1977  lower limit of integration
1978 
1979  b - double precision
1980  upper limit of integration
1981 
1982  on return
1983  result - double precision
1984  approximation to the integral i
1985  result is computed by applying the 21-point
1986  kronrod rule (resk) obtained by optimal addition
1987  of abscissae to the 10-point gauss rule (resg).
1988 
1989  abserr - double precision
1990  estimate of the modulus of the absolute error,
1991  which should not exceed abs(i-result)
1992 
1993  resabs - double precision
1994  approximation to the integral j
1995 
1996  resasc - double precision
1997  approximation to the integral of abs(f-i/(b-a))
1998  over (a,b)
1999 
2000 ***references (none)
2001 ***end prologue dqk21
2002 
2003 
2004 
2005  the abscissae and weights are given for the interval (-1,1).
2006  because of symmetry only the positive abscissae and their
2007  corresponding weights are given.
2008 
2009  xgk - abscissae of the 21-point kronrod rule
2010  xgk(2), xgk(4), ... abscissae of the 10-point
2011  gauss rule
2012  xgk(1), xgk(3), ... abscissae which are optimally
2013  added to the 10-point gauss rule
2014 
2015  wgk - weights of the 21-point kronrod rule
2016 
2017  wg - weights of the 10-point gauss rule
2018 
2019 
2020 gauss quadrature weights and kronron quadrature abscissae and weights
2021 as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
2022 bell labs, nov. 1981.
2023 
2024 
2025 
2026 
2027 
2028  list of major variables
2029  -----------------------
2030 
2031  centr - mid point of the interval
2032  hlgth - half-length of the interval
2033  absc - abscissa
2034  fval* - function value
2035  resg - result of the 10-point gauss formula
2036  resk - result of the 21-point kronrod formula
2037  reskh - approximation to the mean value of f over (a,b),
2038  i.e. to i/(b-a)
2039 
2040 
2041  machine dependent constants
2042  ---------------------------
2043 
2044  epmach is the largest relative spacing.
2045  uflow is the smallest positive magnitude. */
2046 
2047 /* ***first executable statement dqk21 */
2048  epmach = DBL_EPSILON;
2049  uflow = DBL_MIN;
2050 
2051  centr = (*a + *b) * .5;
2052  hlgth = (*b - *a) * .5;
2053  dhlgth = fabs(hlgth);
2054 
2055 /* compute the 21-point kronrod approximation to
2056  the integral, and estimate the absolute error. */
2057 
2058  resg = 0.;
2059  vec[0] = centr;
2060  for (j = 1; j <= 5; ++j) {
2061  jtw = j << 1;
2062  absc = hlgth * xgk[jtw - 1];
2063  vec[(j << 1) - 1] = centr - absc;
2064 /* L5: */
2065  vec[j * 2] = centr + absc;
2066  }
2067  for (j = 1; j <= 5; ++j) {
2068  jtwm1 = (j << 1) - 1;
2069  absc = hlgth * xgk[jtwm1 - 1];
2070  vec[(j << 1) + 9] = centr - absc;
2071  vec[(j << 1) + 10] = centr + absc;
2072  }
2073  f(vec, 21, ex);
2074  fc = vec[0];
2075  resk = wgk[10] * fc;
2076  *resabs = fabs(resk);
2077  for (j = 1; j <= 5; ++j) {
2078  jtw = j << 1;
2079  absc = hlgth * xgk[jtw - 1];
2080  fval1 = vec[(j << 1) - 1];
2081  fval2 = vec[j * 2];
2082  fv1[jtw - 1] = fval1;
2083  fv2[jtw - 1] = fval2;
2084  fsum = fval1 + fval2;
2085  resg += wg[j - 1] * fsum;
2086  resk += wgk[jtw - 1] * fsum;
2087  *resabs += wgk[jtw - 1] * (fabs(fval1) + fabs(fval2));
2088 /* L10: */
2089  }
2090  for (j = 1; j <= 5; ++j) {
2091  jtwm1 = (j << 1) - 1;
2092  absc = hlgth * xgk[jtwm1 - 1];
2093  fval1 = vec[(j << 1) + 9];
2094  fval2 = vec[(j << 1) + 10];
2095  fv1[jtwm1 - 1] = fval1;
2096  fv2[jtwm1 - 1] = fval2;
2097  fsum = fval1 + fval2;
2098  resk += wgk[jtwm1 - 1] * fsum;
2099  *resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
2100 /* L15: */
2101  }
2102  reskh = resk * .5;
2103  *resasc = wgk[10] * fabs(fc - reskh);
2104  for (j = 1; j <= 10; ++j) {
2105  *resasc += wgk[j - 1] * (fabs(fv1[j - 1] - reskh) +
2106  fabs(fv2[j - 1] - reskh));
2107 /* L20: */
2108  }
2109  *result = resk * hlgth;
2110  *resabs *= dhlgth;
2111  *resasc *= dhlgth;
2112  *abserr = fabs((resk - resg) * hlgth);
2113  if (*resasc != 0. && *abserr != 0.) {
2114  *abserr = *resasc * fmin2(1., pow(*abserr * 200. / *resasc, 1.5));
2115  }
2116  if (*resabs > uflow / (epmach * 50.)) {
2117  *abserr = fmax2(epmach * 50. * *resabs, *abserr);
2118  }
2119  return;
2120 } /* rdqk21_ */
2121 
2122 template<class Float> static void rdqpsrt(int *limit, int *last, int *maxerr,
2123  Float *ermax, Float *elist, int *iord, int *nrmax)
2124 {
2125  /* Local variables */
2126  int i, j, k, ido, jbnd, isucc, jupbn;
2127  Float errmin, errmax;
2128 
2129 /* ***begin prologue dqpsrt
2130  ***refer to dqage,dqagie,dqagpe,dqawse
2131  ***routines called (none)
2132  ***revision date 810101 (yymmdd)
2133  ***keywords sequential sorting
2134  ***author piessens,robert,appl. math. & progr. div. - k.u.leuven
2135  de doncker,elise,appl. math. & progr. div. - k.u.leuven
2136  ***purpose this routine maintains the descending ordering in the
2137  list of the local error estimated resulting from the
2138  interval subdivision process. at each call two error
2139  estimates are inserted using the sequential search
2140  method, top-down for the largest error estimate and
2141  bottom-up for the smallest error estimate.
2142  ***description
2143 
2144  ordering routine
2145  standard fortran subroutine
2146  double precision version
2147 
2148  parameters (meaning at output)
2149  limit - int
2150  maximum number of error estimates the list
2151  can contain
2152 
2153  last - int
2154  number of error estimates currently in the list
2155 
2156  maxerr - int
2157  maxerr points to the nrmax-th largest error
2158  estimate currently in the list
2159 
2160  ermax - double precision
2161  nrmax-th largest error estimate
2162  ermax = elist(maxerr)
2163 
2164  elist - double precision
2165  vector of dimension last containing
2166  the error estimates
2167 
2168  iord - int
2169  vector of dimension last, the first k elements
2170  of which contain pointers to the error
2171  estimates, such that
2172  elist(iord(1)),..., elist(iord(k))
2173  form a decreasing sequence, with
2174  k = last if last <= (limit/2+2), and
2175  k = limit+1-last otherwise
2176 
2177  nrmax - int
2178  maxerr = iord(nrmax)
2179 
2180 ***end prologue dqpsrt
2181 */
2182 
2183 
2184  /* Parameter adjustments */
2185  --iord;
2186  --elist;
2187 
2188  /* Function Body */
2189 
2190 /* check whether the list contains more than
2191  two error estimates. */
2192  if (*last <= 2) {
2193  iord[1] = 1;
2194  iord[2] = 2;
2195  goto Last;
2196  }
2197 /* this part of the routine is only executed if, due to a
2198  difficult integrand, subdivision increased the error
2199  estimate. in the normal case the insert procedure should
2200  start after the nrmax-th largest error estimate. */
2201 
2202  errmax = elist[*maxerr];
2203  if (*nrmax > 1) {
2204  ido = *nrmax - 1;
2205  for (i = 1; i <= ido; ++i) {
2206  isucc = iord[*nrmax - 1];
2207  if (errmax <= elist[isucc])
2208  break; /* out of for-loop */
2209  iord[*nrmax] = isucc;
2210  --(*nrmax);
2211  /* L20: */
2212  }
2213  }
2214 
2215 /*L30: compute the number of elements in the list to be maintained
2216  in descending order. this number depends on the number of
2217  subdivisions still allowed. */
2218  if (*last > *limit / 2 + 2)
2219  jupbn = *limit + 3 - *last;
2220  else
2221  jupbn = *last;
2222 
2223  errmin = elist[*last];
2224 
2225 /* insert errmax by traversing the list top-down,
2226  starting comparison from the element elist(iord(nrmax+1)). */
2227 
2228  jbnd = jupbn - 1;
2229  for (i = *nrmax + 1; i <= jbnd; ++i) {
2230  isucc = iord[i];
2231  if (errmax >= elist[isucc]) {/* ***jump out of do-loop */
2232  /* L60: insert errmin by traversing the list bottom-up. */
2233  iord[i - 1] = *maxerr;
2234  for (j = i, k = jbnd; j <= jbnd; j++, k--) {
2235  isucc = iord[k];
2236  if (errmin < elist[isucc]) {
2237  /* goto L80; ***jump out of do-loop */
2238  iord[k + 1] = *last;
2239  goto Last;
2240  }
2241  iord[k + 1] = isucc;
2242  }
2243  iord[i] = *last;
2244  goto Last;
2245  }
2246  iord[i - 1] = isucc;
2247  }
2248 
2249  iord[jbnd] = *maxerr;
2250  iord[jupbn] = *last;
2251 
2252 Last:/* set maxerr and ermax. */
2253 
2254  *maxerr = iord[*nrmax];
2255  *ermax = elist[*maxerr];
2256  return;
2257 } /* rdqpsrt_ */
License: GPL v2