cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
thirdparty_quadpack.cpp
Go to the documentation of this file.
1 #include "cddefines.h"
2 #include "thirdparty_quadpack.h"
3 
4 /* translated by f2c (version 20100827). */
5 
6 STATIC long xerror_(const char *mesg, const long *nmesg,
7  const long *nerr, const long *level, long)
8 {
9  DEBUG_ENTRY("xerror()");
10  if (*level > 0)
11  {
12  fprintf(ioQQQ,"Error [%ld]: %*s\n",*nerr,int(*nmesg),mesg);
14  }
15  else
16  {
17  fprintf(ioQQQ,"Warning [%ld]: %*s\n",*nerr,int(*nmesg),mesg);
18  }
19  return 0;
20 }
21 static sys_float r1mach(long i)
22 {
23  DEBUG_ENTRY("r1mach()");
24  // DOUBLE-PRECISION MACHINE CONSTANTS
25  // D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
26  // D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
27  // D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
28  // D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
29  // D1MACH( 5) = LOG10(B)
30  switch (i)
31  {
32  case 1:
33  return FLT_MIN;
34  case 2:
35  return FLT_MAX;
36  case 3:
37  return FLT_EPSILON/FLT_RADIX;
38  case 4:
39  return FLT_EPSILON;
40  case 5:
41  return log10(double(FLT_RADIX));
42  default:
43  fprintf(stderr,"Error in input to r1mach");
45  }
46 }
47 static double d1mach(long i)
48 {
49  DEBUG_ENTRY("d1mach()");
50  // DOUBLE-PRECISION MACHINE CONSTANTS
51  // D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
52  // D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
53  // D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
54  // D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
55  // D1MACH( 5) = LOG10(B)
56  switch (i)
57  {
58  case 1:
59  return DBL_MIN;
60  case 2:
61  return DBL_MAX;
62  case 3:
63  return DBL_EPSILON/FLT_RADIX;
64  case 4:
65  return DBL_EPSILON;
66  case 5:
67  return log10(double(FLT_RADIX));
68  default:
70  }
71 }
72 
73 /* Table of constant values */
74 
75 const static long c__4 = 4;
76 const static long c__1 = 1;
77 const static long c__26 = 26;
78 const static double c_b20 = 0.;
79 const static double c_b21 = 1.;
80 const static long c__2 = 2;
81 const static long c__0 = 0;
82 const static double c_b270 = 1.5;
83 const static double c_b320 = 2.;
84 const static sys_float c_b390 = 0.f;
85 const static sys_float c_b391 = 1.f;
86 
87 static void sgtsl_(const long *, sys_float*, sys_float*,
88  sys_float*, sys_float *, long *);
89 static void dgtsl_(const long *, double*, double*, double*,
90  double *, long *);
91 
92 void dqage_(const D_fp& f, const double *a, const double *b, const double *epsabs,
93  const double *epsrel, const long *key, const long *limit,
94  double *result, double *abserr, long *neval, long *ier, double *
95  alist__, double *blist, double *rlist, double *elist,
96  long *iord, long *last)
97 {
98  /* System generated locals */
99  long i__1;
100  double d__1, d__2;
101 
102  /* Local variables */
103  long k;
104  double a1, a2, b1, b2, area;
105  long keyf;
106  double area1, area2, area12, erro12, defab1, defab2;
107  long nrmax;
108  double uflow;
109  long iroff1, iroff2;
110  double error1, error2, defabs, epmach, errbnd, resabs, errmax;
111  long maxerr;
112  double errsum;
113 
114  /* ***begin prologue dqage */
115  /* ***date written 800101 (yymmdd) */
116  /* ***revision date 830518 (yymmdd) */
117  /* ***category no. h2a1a1 */
118  /* ***keywords automatic integrator, general-purpose, */
119  /* integrand examinator, globally adaptive, */
120  /* gauss-kronrod */
121  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
122  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
123  /* ***purpose the routine calculates an approximation result to a given */
124  /* definite integral i = integral of f over (a,b), */
125  /* hopefully satisfying following claim for accuracy */
126  /* abs(i-reslt).le.max(epsabs,epsrel*abs(i)). */
127  /* ***description */
128 
129  /* computation of a definite integral */
130  /* standard fortran subroutine */
131  /* double precision version */
132 
133  /* parameters */
134  /* on entry */
135  /* f - double precision */
136  /* function subprogram defining the integrand */
137  /* function f(x). the actual name for f needs to be */
138  /* declared e x t e r n a l in the driver program. */
139 
140  /* a - double precision */
141  /* lower limit of integration */
142 
143  /* b - double precision */
144  /* upper limit of integration */
145 
146  /* epsabs - double precision */
147  /* absolute accuracy requested */
148  /* epsrel - double precision */
149  /* relative accuracy requested */
150  /* if epsabs.le.0 */
151  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
152  /* the routine will end with ier = 6. */
153 
154  /* key - long */
155  /* key for choice of local integration rule */
156  /* a gauss-kronrod pair is used with */
157  /* 7 - 15 points if key.lt.2, */
158  /* 10 - 21 points if key = 2, */
159  /* 15 - 31 points if key = 3, */
160  /* 20 - 41 points if key = 4, */
161  /* 25 - 51 points if key = 5, */
162  /* 30 - 61 points if key.gt.5. */
163 
164  /* limit - long */
165  /* gives an upperbound on the number of subintervals */
166  /* in the partition of (a,b), limit.ge.1. */
167 
168  /* on return */
169  /* result - double precision */
170  /* approximation to the integral */
171 
172  /* abserr - double precision */
173  /* estimate of the modulus of the absolute error, */
174  /* which should equal or exceed abs(i-result) */
175 
176  /* neval - long */
177  /* number of integrand evaluations */
178 
179  /* ier - long */
180  /* ier = 0 normal and reliable termination of the */
181  /* routine. it is assumed that the requested */
182  /* accuracy has been achieved. */
183  /* ier.gt.0 abnormal termination of the routine */
184  /* the estimates for result and error are */
185  /* less reliable. it is assumed that the */
186  /* requested accuracy has not been achieved. */
187  /* error messages */
188  /* ier = 1 maximum number of subdivisions allowed */
189  /* has been achieved. one can allow more */
190  /* subdivisions by increasing the value */
191  /* of limit. */
192  /* however, if this yields no improvement it */
193  /* is rather advised to analyze the integrand */
194  /* in order to determine the integration */
195  /* difficulties. if the position of a local */
196  /* difficulty can be determined(e.g. */
197  /* singularity, discontinuity within the */
198  /* interval) one will probably gain from */
199  /* splitting up the interval at this point */
200  /* and calling the integrator on the */
201  /* subranges. if possible, an appropriate */
202  /* special-purpose integrator should be used */
203  /* which is designed for handling the type of */
204  /* difficulty involved. */
205  /* = 2 the occurrence of roundoff error is */
206  /* detected, which prevents the requested */
207  /* tolerance from being achieved. */
208  /* = 3 extremely bad integrand behaviour occurs */
209  /* at some points of the integration */
210  /* interval. */
211  /* = 6 the input is invalid, because */
212  /* (epsabs.le.0 and */
213  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
214  /* result, abserr, neval, last, rlist(1) , */
215  /* elist(1) and iord(1) are set to zero. */
216  /* alist(1) and blist(1) are set to a and b */
217  /* respectively. */
218 
219  /* alist - double precision */
220  /* vector of dimension at least limit, the first */
221  /* last elements of which are the left */
222  /* end points of the subintervals in the partition */
223  /* of the given integration range (a,b) */
224 
225  /* blist - double precision */
226  /* vector of dimension at least limit, the first */
227  /* last elements of which are the right */
228  /* end points of the subintervals in the partition */
229  /* of the given integration range (a,b) */
230 
231  /* rlist - double precision */
232  /* vector of dimension at least limit, the first */
233  /* last elements of which are the */
234  /* integral approximations on the subintervals */
235 
236  /* elist - double precision */
237  /* vector of dimension at least limit, the first */
238  /* last elements of which are the moduli of the */
239  /* absolute error estimates on the subintervals */
240 
241  /* iord - int */
242  /* vector of dimension at least limit, the first k */
243  /* elements of which are pointers to the */
244  /* error estimates over the subintervals, */
245  /* such that elist(iord(1)), ..., */
246  /* elist(iord(k)) form a decreasing sequence, */
247  /* with k = last if last.le.(limit/2+2), and */
248  /* k = limit+1-last otherwise */
249 
250  /* last - int */
251  /* number of subintervals actually produced in the */
252  /* subdivision process */
253 
254  /* ***references (none) */
255  /* ***routines called d1mach,dqk15,dqk21,dqk31, */
256  /* dqk41,dqk51,dqk61,dqpsrt */
257  /* ***end prologue dqage */
258 
259 
260 
261 
262  /* list of major variables */
263  /* ----------------------- */
264 
265  /* alist - list of left end points of all subintervals */
266  /* considered up to now */
267  /* blist - list of right end points of all subintervals */
268  /* considered up to now */
269  /* rlist(i) - approximation to the integral over */
270  /* (alist(i),blist(i)) */
271  /* elist(i) - error estimate applying to rlist(i) */
272  /* maxerr - pointer to the interval with largest */
273  /* error estimate */
274  /* errmax - elist(maxerr) */
275  /* area - sum of the integrals over the subintervals */
276  /* errsum - sum of the errors over the subintervals */
277  /* errbnd - requested accuracy max(epsabs,epsrel* */
278  /* abs(result)) */
279  /* *****1 - variable for the left subinterval */
280  /* *****2 - variable for the right subinterval */
281  /* last - index for subdivision */
282 
283 
284  /* machine dependent constants */
285  /* --------------------------- */
286 
287  /* epmach is the largest relative spacing. */
288  /* uflow is the smallest positive magnitude. */
289 
290  /* ***first executable statement dqage */
291  /* Parameter adjustments */
292  --iord;
293  --elist;
294  --rlist;
295  --blist;
296  --alist__;
297 
298  /* Function Body */
299  epmach = d1mach(c__4);
300  uflow = d1mach(c__1);
301 
302  /* test on validity of parameters */
303  /* ------------------------------ */
304 
305  *ier = 0;
306  *neval = 0;
307  *last = 0;
308  *result = 0.;
309  *abserr = 0.;
310  alist__[1] = *a;
311  blist[1] = *b;
312  rlist[1] = 0.;
313  elist[1] = 0.;
314  iord[1] = 0;
315  /* Computing MAX */
316  d__1 = epmach * 50.;
317  if (*epsabs <= 0. && *epsrel < max(d__1,5e-29)) {
318  *ier = 6;
319  }
320  if (*ier == 6) {
321  goto L999;
322  }
323 
324  /* first approximation to the integral */
325  /* ----------------------------------- */
326 
327  keyf = *key;
328  if (*key <= 0) {
329  keyf = 1;
330  }
331  if (*key >= 7) {
332  keyf = 6;
333  }
334  *neval = 0;
335  if (keyf == 1) {
336  dqk15_(f, a, b, result, abserr, &defabs, &resabs);
337  }
338  if (keyf == 2) {
339  dqk21_(f, a, b, result, abserr, &defabs, &resabs);
340  }
341  if (keyf == 3) {
342  dqk31_(f, a, b, result, abserr, &defabs, &resabs);
343  }
344  if (keyf == 4) {
345  dqk41_(f, a, b, result, abserr, &defabs, &resabs);
346  }
347  if (keyf == 5) {
348  dqk51_(f, a, b, result, abserr, &defabs, &resabs);
349  }
350  if (keyf == 6) {
351  dqk61_(f, a, b, result, abserr, &defabs, &resabs);
352  }
353  *last = 1;
354  rlist[1] = *result;
355  elist[1] = *abserr;
356  iord[1] = 1;
357 
358  /* test on accuracy. */
359 
360  /* Computing MAX */
361  d__1 = *epsabs, d__2 = *epsrel * fabs(*result);
362  errbnd = max(d__1,d__2);
363  if (*abserr <= epmach * 50. * defabs && *abserr > errbnd) {
364  *ier = 2;
365  }
366  if (*limit == 1) {
367  *ier = 1;
368  }
369  if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) || *abserr == 0.)
370  {
371  goto L60;
372  }
373 
374  /* initialization */
375  /* -------------- */
376 
377 
378  errmax = *abserr;
379  maxerr = 1;
380  area = *result;
381  errsum = *abserr;
382  nrmax = 1;
383  iroff1 = 0;
384  iroff2 = 0;
385 
386  /* main do-loop */
387  /* ------------ */
388 
389  i__1 = *limit;
390  for (*last = 2; *last <= i__1; ++(*last)) {
391 
392  /* bisect the subinterval with the largest error estimate. */
393 
394  a1 = alist__[maxerr];
395  b1 = (alist__[maxerr] + blist[maxerr]) * .5;
396  a2 = b1;
397  b2 = blist[maxerr];
398  if (keyf == 1) {
399  dqk15_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
400  }
401  if (keyf == 2) {
402  dqk21_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
403  }
404  if (keyf == 3) {
405  dqk31_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
406  }
407  if (keyf == 4) {
408  dqk41_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
409  }
410  if (keyf == 5) {
411  dqk51_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
412  }
413  if (keyf == 6) {
414  dqk61_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
415  }
416  if (keyf == 1) {
417  dqk15_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
418  }
419  if (keyf == 2) {
420  dqk21_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
421  }
422  if (keyf == 3) {
423  dqk31_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
424  }
425  if (keyf == 4) {
426  dqk41_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
427  }
428  if (keyf == 5) {
429  dqk51_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
430  }
431  if (keyf == 6) {
432  dqk61_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
433  }
434 
435  /* improve previous approximations to integral */
436  /* and error and test for accuracy. */
437 
438  ++(*neval);
439  area12 = area1 + area2;
440  erro12 = error1 + error2;
441  errsum = errsum + erro12 - errmax;
442  area = area + area12 - rlist[maxerr];
443  if (defab1 == error1 || defab2 == error2) {
444  goto L5;
445  }
446  if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) <= fabs(area12) * 1e-5
447  && erro12 >= errmax * .99) {
448  ++iroff1;
449  }
450  if (*last > 10 && erro12 > errmax) {
451  ++iroff2;
452  }
453  L5:
454  rlist[maxerr] = area1;
455  rlist[*last] = area2;
456  /* Computing MAX */
457  d__1 = *epsabs, d__2 = *epsrel * fabs(area);
458  errbnd = max(d__1,d__2);
459  if (errsum <= errbnd) {
460  goto L8;
461  }
462 
463  /* test for roundoff error and eventually set error flag. */
464 
465  if (iroff1 >= 6 || iroff2 >= 20) {
466  *ier = 2;
467  }
468 
469  /* set error flag in the case that the number of subintervals */
470  /* equals limit. */
471 
472  if (*last == *limit) {
473  *ier = 1;
474  }
475 
476  /* set error flag in the case of bad integrand behaviour */
477  /* at a point of the integration range. */
478 
479  /* Computing MAX */
480  d__1 = fabs(a1), d__2 = fabs(b2);
481  if (max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
482  {
483  *ier = 3;
484  }
485 
486  /* append the newly-created intervals to the list. */
487 
488  L8:
489  if (error2 > error1) {
490  goto L10;
491  }
492  alist__[*last] = a2;
493  blist[maxerr] = b1;
494  blist[*last] = b2;
495  elist[maxerr] = error1;
496  elist[*last] = error2;
497  goto L20;
498  L10:
499  alist__[maxerr] = a2;
500  alist__[*last] = a1;
501  blist[*last] = b1;
502  rlist[maxerr] = area2;
503  rlist[*last] = area1;
504  elist[maxerr] = error2;
505  elist[*last] = error1;
506 
507  /* call subroutine dqpsrt to maintain the descending ordering */
508  /* in the list of error estimates and select the subinterval */
509  /* with the largest error estimate (to be bisected next). */
510 
511  L20:
512  dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
513  /* ***jump out of do-loop */
514  if (*ier != 0 || errsum <= errbnd) {
515  goto L40;
516  }
517  /* L30: */
518  }
519 
520  /* compute final result. */
521  /* --------------------- */
522 
523 L40:
524  *result = 0.;
525  i__1 = *last;
526  for (k = 1; k <= i__1; ++k) {
527  *result += rlist[k];
528  /* L50: */
529  }
530  *abserr = errsum;
531 L60:
532  if (keyf != 1) {
533  *neval = (keyf * 10 + 1) * ((*neval << 1) + 1);
534  }
535  if (keyf == 1) {
536  *neval = *neval * 30 + 15;
537  }
538 L999:
539  return;
540 } /* dqage_ */
541 
542 void dqag_(const D_fp& f, const double *a, const double *b, const double *epsabs,
543  const double *epsrel, const long *key, double *result,
544  double *abserr, long *neval, long *ier, long *limit,
545  const long *lenw, long *last, long *iwork, double *work)
546 {
547  long l1, l2, l3, lvl;
548 
549  /* ***begin prologue dqag */
550  /* ***date written 800101 (yymmdd) */
551  /* ***revision date 830518 (yymmdd) */
552  /* ***category no. h2a1a1 */
553  /* ***keywords automatic integrator, general-purpose, */
554  /* integrand examinator, globally adaptive, */
555  /* gauss-kronrod */
556  /* ***author piessens,robert,appl. math. & progr. div - k.u.leuven */
557  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
558  /* ***purpose the routine calculates an approximation result to a given */
559  /* definite integral i = integral of f over (a,b), */
560  /* hopefully satisfying following claim for accuracy */
561  /* abs(i-result)le.max(epsabs,epsrel*abs(i)). */
562  /* ***description */
563 
564  /* computation of a definite integral */
565  /* standard fortran subroutine */
566  /* double precision version */
567 
568  /* f - double precision */
569  /* function subprogam defining the integrand */
570  /* function f(x). the actual name for f needs to be */
571  /* declared e x t e r n a l in the driver program. */
572 
573  /* a - double precision */
574  /* lower limit of integration */
575 
576  /* b - double precision */
577  /* upper limit of integration */
578 
579  /* epsabs - double precision */
580  /* absolute accoracy requested */
581  /* epsrel - double precision */
582  /* relative accuracy requested */
583  /* if epsabs.le.0 */
584  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
585  /* the routine will end with ier = 6. */
586 
587  /* key - int */
588  /* key for choice of local integration rule */
589  /* a gauss-kronrod pair is used with */
590  /* 7 - 15 points if key.lt.2, */
591  /* 10 - 21 points if key = 2, */
592  /* 15 - 31 points if key = 3, */
593  /* 20 - 41 points if key = 4, */
594  /* 25 - 51 points if key = 5, */
595  /* 30 - 61 points if key.gt.5. */
596 
597  /* on return */
598  /* result - double precision */
599  /* approximation to the integral */
600 
601  /* abserr - double precision */
602  /* estimate of the modulus of the absolute error, */
603  /* which should equal or exceed abs(i-result) */
604 
605  /* neval - int */
606  /* number of integrand evaluations */
607 
608  /* ier - long */
609  /* ier = 0 normal and reliable termination of the */
610  /* routine. it is assumed that the requested */
611  /* accuracy has been achieved. */
612  /* ier.gt.0 abnormal termination of the routine */
613  /* the estimates for result and error are */
614  /* less reliable. it is assumed that the */
615  /* requested accuracy has not been achieved. */
616  /* error messages */
617  /* ier = 1 maximum number of subdivisions allowed */
618  /* has been achieved. one can allow more */
619  /* subdivisions by increasing the value of */
620  /* limit (and taking the according dimension */
621  /* adjustments into account). however, if */
622  /* this yield no improvement it is advised */
623  /* to analyze the integrand in order to */
624  /* determine the integration difficulaties. */
625  /* if the position of a local difficulty can */
626  /* be determined (i.e.singularity, */
627  /* discontinuity within the interval) one */
628  /* will probably gain from splitting up the */
629  /* interval at this point and calling the */
630  /* integrator on the subranges. if possible, */
631  /* an appropriate special-purpose integrator */
632  /* should be used which is designed for */
633  /* handling the type of difficulty involved. */
634  /* = 2 the occurrence of roundoff error is */
635  /* detected, which prevents the requested */
636  /* tolerance from being achieved. */
637  /* = 3 extremely bad integrand behaviour occurs */
638  /* at some points of the integration */
639  /* interval. */
640  /* = 6 the input is invalid, because */
641  /* (epsabs.le.0 and */
642  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
643  /* or limit.lt.1 or lenw.lt.limit*4. */
644  /* result, abserr, neval, last are set */
645  /* to zero. */
646  /* except when lenw is invalid, iwork(1), */
647  /* work(limit*2+1) and work(limit*3+1) are */
648  /* set to zero, work(1) is set to a and */
649  /* work(limit+1) to b. */
650 
651  /* dimensioning parameters */
652  /* limit - int */
653  /* dimensioning parameter for iwork */
654  /* limit determines the maximum number of subintervals */
655  /* in the partition of the given integration interval */
656  /* (a,b), limit.ge.1. */
657  /* if limit.lt.1, the routine will end with ier = 6. */
658 
659  /* lenw - long */
660  /* dimensioning parameter for work */
661  /* lenw must be at least limit*4. */
662  /* if lenw.lt.limit*4, the routine will end with */
663  /* ier = 6. */
664 
665  /* last - int */
666  /* on return, last equals the number of subintervals */
667  /* produced in the subdiviosion process, which */
668  /* determines the number of significant elements */
669  /* actually in the work arrays. */
670 
671  /* work arrays */
672  /* iwork - int */
673  /* vector of dimension at least limit, the first k */
674  /* elements of which contain pointers to the error */
675  /* estimates over the subintervals, such that */
676  /* work(limit*3+iwork(1)),... , work(limit*3+iwork(k)) */
677  /* form a decreasing sequence with k = last if */
678  /* last.le.(limit/2+2), and k = limit+1-last otherwise */
679 
680  /* work - double precision */
681  /* vector of dimension at least lenw */
682  /* on return */
683  /* work(1), ..., work(last) contain the left end */
684  /* points of the subintervals in the partition of */
685  /* (a,b), */
686  /* work(limit+1), ..., work(limit+last) contain the */
687  /* right end points, */
688  /* work(limit*2+1), ..., work(limit*2+last) contain */
689  /* the integral approximations over the subintervals, */
690  /* work(limit*3+1), ..., work(limit*3+last) contain */
691  /* the error estimates. */
692 
693  /* ***references (none) */
694  /* ***routines called dqage,xerror */
695  /* ***end prologue dqag */
696 
697 
698 
699  /* check validity of lenw. */
700 
701  /* ***first executable statement dqag */
702  /* Parameter adjustments */
703  --iwork;
704  --work;
705 
706  /* Function Body */
707  *ier = 6;
708  *neval = 0;
709  *last = 0;
710  *result = 0.;
711  *abserr = 0.;
712  if (*limit < 1 || *lenw < *limit << 2) {
713  goto L10;
714  }
715 
716  /* prepare call for dqage. */
717 
718  l1 = *limit + 1;
719  l2 = *limit + l1;
720  l3 = *limit + l2;
721 
722  dqage_(f, a, b, epsabs, epsrel, key, limit, result, abserr, neval,
723  ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
724 
725  /* call error handler if necessary. */
726 
727  lvl = 0;
728 L10:
729  if (*ier == 6) {
730  lvl = 1;
731  }
732  if (*ier != 0) {
733  xerror_("abnormal return from dqag ", &c__26, ier, &lvl, 26);
734  }
735  return;
736 } /* dqag_ */
737 
738 void dqagie_(const D_fp& f, const double *bound, const long *inf,
739  const double *epsabs, const double *epsrel, const long *limit,
740  double *result, double *abserr, long *neval, long *ier,
741  double *alist__, double *blist, double *rlist, double *elist,
742  long *iord, long *last)
743 {
744  /* System generated locals */
745  long i__1, i__2;
746  double d__1, d__2;
747 
748  /* Local variables */
749  long k;
750  double a1, a2, b1, b2;
751  long id;
752  double area, dres;
753  long ksgn;
754  double boun;
755  long nres;
756  double area1, area2, area12;
757  double small=0., erro12;
758  long ierro;
759  double defab1, defab2;
760  long ktmin, nrmax;
761  double oflow, uflow;
762  bool noext;
763  long iroff1, iroff2, iroff3;
764  double res3la[3], error1, error2, rlist2[52];
765  long numrl2;
766  double defabs, epmach, erlarg=0., abseps, correc=0., errbnd, resabs;
767  long jupbnd;
768  double erlast, errmax;
769  long maxerr;
770  double reseps;
771  bool extrap;
772  double ertest=0., errsum;
773 
774  /* ***begin prologue dqagie */
775  /* ***date written 800101 (yymmdd) */
776  /* ***revision date 830518 (yymmdd) */
777  /* ***category no. h2a3a1,h2a4a1 */
778  /* ***keywords automatic integrator, infinite intervals, */
779  /* general-purpose, transformation, extrapolation, */
780  /* globally adaptive */
781  /* ***author piessens,robert,appl. math & progr. div - k.u.leuven */
782  /* de doncker,elise,appl. math & progr. div - k.u.leuven */
783  /* ***purpose the routine calculates an approximation result to a given */
784  /* integral i = integral of f over (bound,+infinity) */
785  /* or i = integral of f over (-infinity,bound) */
786  /* or i = integral of f over (-infinity,+infinity), */
787  /* hopefully satisfying following claim for accuracy */
788  /* abs(i-result).le.max(epsabs,epsrel*abs(i)) */
789  /* ***description */
790 
791  /* integration over infinite intervals */
792  /* standard fortran subroutine */
793 
794  /* f - double precision */
795  /* function subprogram defining the integrand */
796  /* function f(x). the actual name for f needs to be */
797  /* declared e x t e r n a l in the driver program. */
798 
799  /* bound - double precision */
800  /* finite bound of integration range */
801  /* (has no meaning if interval is doubly-infinite) */
802 
803  /* inf - double precision */
804  /* indicating the kind of integration range involved */
805  /* inf = 1 corresponds to (bound,+infinity), */
806  /* inf = -1 to (-infinity,bound), */
807  /* inf = 2 to (-infinity,+infinity). */
808 
809  /* epsabs - double precision */
810  /* absolute accuracy requested */
811  /* epsrel - double precision */
812  /* relative accuracy requested */
813  /* if epsabs.le.0 */
814  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
815  /* the routine will end with ier = 6. */
816 
817  /* limit - long */
818  /* gives an upper bound on the number of subintervals */
819  /* in the partition of (a,b), limit.ge.1 */
820 
821  /* on return */
822  /* result - double precision */
823  /* approximation to the integral */
824 
825  /* abserr - double precision */
826  /* estimate of the modulus of the absolute error, */
827  /* which should equal or exceed abs(i-result) */
828 
829  /* neval - long */
830  /* number of integrand evaluations */
831 
832  /* ier - long */
833  /* ier = 0 normal and reliable termination of the */
834  /* routine. it is assumed that the requested */
835  /* accuracy has been achieved. */
836  /* - ier.gt.0 abnormal termination of the routine. the */
837  /* estimates for result and error are less */
838  /* reliable. it is assumed that the requested */
839  /* accuracy has not been achieved. */
840  /* error messages */
841  /* ier = 1 maximum number of subdivisions allowed */
842  /* has been achieved. one can allow more */
843  /* subdivisions by increasing the value of */
844  /* limit (and taking the according dimension */
845  /* adjustments into account). however,if */
846  /* this yields no improvement it is advised */
847  /* to analyze the integrand in order to */
848  /* determine the integration difficulties. */
849  /* if the position of a local difficulty can */
850  /* be determined (e.g. singularity, */
851  /* discontinuity within the interval) one */
852  /* will probably gain from splitting up the */
853  /* interval at this point and calling the */
854  /* integrator on the subranges. if possible, */
855  /* an appropriate special-purpose integrator */
856  /* should be used, which is designed for */
857  /* handling the type of difficulty involved. */
858  /* = 2 the occurrence of roundoff error is */
859  /* detected, which prevents the requested */
860  /* tolerance from being achieved. */
861  /* the error may be under-estimated. */
862  /* = 3 extremely bad integrand behaviour occurs */
863  /* at some points of the integration */
864  /* interval. */
865  /* = 4 the algorithm does not converge. */
866  /* roundoff error is detected in the */
867  /* extrapolation table. */
868  /* it is assumed that the requested tolerance */
869  /* cannot be achieved, and that the returned */
870  /* result is the best which can be obtained. */
871  /* = 5 the integral is probably divergent, or */
872  /* slowly convergent. it must be noted that */
873  /* divergence can occur with any other value */
874  /* of ier. */
875  /* = 6 the input is invalid, because */
876  /* (epsabs.le.0 and */
877  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
878  /* result, abserr, neval, last, rlist(1), */
879  /* elist(1) and iord(1) are set to zero. */
880  /* alist(1) and blist(1) are set to 0 */
881  /* and 1 respectively. */
882 
883  /* alist - double precision */
884  /* vector of dimension at least limit, the first */
885  /* last elements of which are the left */
886  /* end points of the subintervals in the partition */
887  /* of the transformed integration range (0,1). */
888 
889  /* blist - double precision */
890  /* vector of dimension at least limit, the first */
891  /* last elements of which are the right */
892  /* end points of the subintervals in the partition */
893  /* of the transformed integration range (0,1). */
894 
895  /* rlist - double precision */
896  /* vector of dimension at least limit, the first */
897  /* last elements of which are the integral */
898  /* approximations on the subintervals */
899 
900  /* elist - double precision */
901  /* vector of dimension at least limit, the first */
902  /* last elements of which are the moduli of the */
903  /* absolute error estimates on the subintervals */
904 
905  /* iord - long */
906  /* vector of dimension limit, the first k */
907  /* elements of which are pointers to the */
908  /* error estimates over the subintervals, */
909  /* such that elist(iord(1)), ..., elist(iord(k)) */
910  /* form a decreasing sequence, with k = last */
911  /* if last.le.(limit/2+2), and k = limit+1-last */
912  /* otherwise */
913 
914  /* last - long */
915  /* number of subintervals actually produced */
916  /* in the subdivision process */
917 
918  /* ***references (none) */
919  /* ***routines called d1mach,dqelg,dqk15i,dqpsrt */
920  /* ***end prologue dqagie */
921 
922 
923 
924  /* the dimension of rlist2 is determined by the value of */
925  /* limexp in subroutine dqelg. */
926 
927 
928  /* list of major variables */
929  /* ----------------------- */
930 
931  /* alist - list of left end points of all subintervals */
932  /* considered up to now */
933  /* blist - list of right end points of all subintervals */
934  /* considered up to now */
935  /* rlist(i) - approximation to the integral over */
936  /* (alist(i),blist(i)) */
937  /* rlist2 - array of dimension at least (limexp+2), */
938  /* containing the part of the epsilon table */
939  /* wich is still needed for further computations */
940  /* elist(i) - error estimate applying to rlist(i) */
941  /* maxerr - pointer to the interval with largest error */
942  /* estimate */
943  /* errmax - elist(maxerr) */
944  /* erlast - error on the interval currently subdivided */
945  /* (before that subdivision has taken place) */
946  /* area - sum of the integrals over the subintervals */
947  /* errsum - sum of the errors over the subintervals */
948  /* errbnd - requested accuracy max(epsabs,epsrel* */
949  /* abs(result)) */
950  /* *****1 - variable for the left subinterval */
951  /* *****2 - variable for the right subinterval */
952  /* last - index for subdivision */
953  /* nres - number of calls to the extrapolation routine */
954  /* numrl2 - number of elements currently in rlist2. if an */
955  /* appropriate approximation to the compounded */
956  /* integral has been obtained, it is put in */
957  /* rlist2(numrl2) after numrl2 has been increased */
958  /* by one. */
959  /* small - length of the smallest interval considered up */
960  /* to now, multiplied by 1.5 */
961  /* erlarg - sum of the errors over the intervals larger */
962  /* than the smallest interval considered up to now */
963  /* extrap - bool variable denoting that the routine */
964  /* is attempting to perform extrapolation. i.e. */
965  /* before subdividing the smallest interval we */
966  /* try to decrease the value of erlarg. */
967  /* noext - bool variable denoting that extrapolation */
968  /* is no longer allowed (true-value) */
969 
970  /* machine dependent constants */
971  /* --------------------------- */
972 
973  /* epmach is the largest relative spacing. */
974  /* uflow is the smallest positive magnitude. */
975  /* oflow is the largest positive magnitude. */
976 
977  /* ***first executable statement dqagie */
978  /* Parameter adjustments */
979  --iord;
980  --elist;
981  --rlist;
982  --blist;
983  --alist__;
984 
985  /* Function Body */
986  epmach = d1mach(c__4);
987 
988  /* test on validity of parameters */
989  /* ----------------------------- */
990 
991  *ier = 0;
992  *neval = 0;
993  *last = 0;
994  *result = 0.;
995  *abserr = 0.;
996  alist__[1] = 0.;
997  blist[1] = 1.;
998  rlist[1] = 0.;
999  elist[1] = 0.;
1000  iord[1] = 0;
1001  /* Computing MAX */
1002  d__1 = epmach * 50.;
1003  if (*epsabs <= 0. && *epsrel < max(d__1,5e-29)) {
1004  *ier = 6;
1005  }
1006  if (*ier == 6) {
1007  goto L999;
1008  }
1009 
1010 
1011  /* first approximation to the integral */
1012  /* ----------------------------------- */
1013 
1014  /* determine the interval to be mapped onto (0,1). */
1015  /* if inf = 2 the integral is computed as i = i1+i2, where */
1016  /* i1 = integral of f over (-infinity,0), */
1017  /* i2 = integral of f over (0,+infinity). */
1018 
1019  boun = *bound;
1020  if (*inf == 2) {
1021  boun = 0.;
1022  }
1023  dqk15i_(f, &boun, inf, &c_b20, &c_b21, result, abserr, &defabs, &
1024  resabs);
1025 
1026  /* test on accuracy */
1027 
1028  *last = 1;
1029  rlist[1] = *result;
1030  elist[1] = *abserr;
1031  iord[1] = 1;
1032  dres = fabs(*result);
1033  /* Computing MAX */
1034  d__1 = *epsabs, d__2 = *epsrel * dres;
1035  errbnd = max(d__1,d__2);
1036  if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) {
1037  *ier = 2;
1038  }
1039  if (*limit == 1) {
1040  *ier = 1;
1041  }
1042  if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) || *abserr == 0.)
1043  {
1044  goto L130;
1045  }
1046 
1047  /* initialization */
1048  /* -------------- */
1049 
1050  uflow = d1mach(c__1);
1051  oflow = d1mach(c__2);
1052  rlist2[0] = *result;
1053  errmax = *abserr;
1054  maxerr = 1;
1055  area = *result;
1056  errsum = *abserr;
1057  *abserr = oflow;
1058  nrmax = 1;
1059  nres = 0;
1060  ktmin = 0;
1061  numrl2 = 2;
1062  extrap = false;
1063  noext = false;
1064  ierro = 0;
1065  iroff1 = 0;
1066  iroff2 = 0;
1067  iroff3 = 0;
1068  ksgn = -1;
1069  if (dres >= (1. - epmach * 50.) * defabs) {
1070  ksgn = 1;
1071  }
1072 
1073  /* main do-loop */
1074  /* ------------ */
1075 
1076  i__1 = *limit;
1077  for (*last = 2; *last <= i__1; ++(*last)) {
1078 
1079  /* bisect the subinterval with nrmax-th largest error estimate. */
1080 
1081  a1 = alist__[maxerr];
1082  b1 = (alist__[maxerr] + blist[maxerr]) * .5;
1083  a2 = b1;
1084  b2 = blist[maxerr];
1085  erlast = errmax;
1086  dqk15i_(f, &boun, inf, &a1, &b1, &area1, &error1, &resabs, &
1087  defab1);
1088  dqk15i_(f, &boun, inf, &a2, &b2, &area2, &error2, &resabs, &
1089  defab2);
1090 
1091  /* improve previous approximations to integral */
1092  /* and error and test for accuracy. */
1093 
1094  area12 = area1 + area2;
1095  erro12 = error1 + error2;
1096  errsum = errsum + erro12 - errmax;
1097  area = area + area12 - rlist[maxerr];
1098  if (defab1 == error1 || defab2 == error2) {
1099  goto L15;
1100  }
1101  if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
1102  erro12 < errmax * .99) {
1103  goto L10;
1104  }
1105  if (extrap) {
1106  ++iroff2;
1107  }
1108  if (! extrap) {
1109  ++iroff1;
1110  }
1111  L10:
1112  if (*last > 10 && erro12 > errmax) {
1113  ++iroff3;
1114  }
1115  L15:
1116  rlist[maxerr] = area1;
1117  rlist[*last] = area2;
1118  /* Computing MAX */
1119  d__1 = *epsabs, d__2 = *epsrel * fabs(area);
1120  errbnd = max(d__1,d__2);
1121 
1122  /* test for roundoff error and eventually set error flag. */
1123 
1124  if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
1125  *ier = 2;
1126  }
1127  if (iroff2 >= 5) {
1128  ierro = 3;
1129  }
1130 
1131  /* set error flag in the case that the number of */
1132  /* subintervals equals limit. */
1133 
1134  if (*last == *limit) {
1135  *ier = 1;
1136  }
1137 
1138  /* set error flag in the case of bad integrand behaviour */
1139  /* at some points of the integration range. */
1140 
1141  /* Computing MAX */
1142  d__1 = fabs(a1), d__2 = fabs(b2);
1143  if (max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
1144  {
1145  *ier = 4;
1146  }
1147 
1148  /* append the newly-created intervals to the list. */
1149 
1150  if (error2 > error1) {
1151  goto L20;
1152  }
1153  alist__[*last] = a2;
1154  blist[maxerr] = b1;
1155  blist[*last] = b2;
1156  elist[maxerr] = error1;
1157  elist[*last] = error2;
1158  goto L30;
1159  L20:
1160  alist__[maxerr] = a2;
1161  alist__[*last] = a1;
1162  blist[*last] = b1;
1163  rlist[maxerr] = area2;
1164  rlist[*last] = area1;
1165  elist[maxerr] = error2;
1166  elist[*last] = error1;
1167 
1168  /* call subroutine dqpsrt to maintain the descending ordering */
1169  /* in the list of error estimates and select the subinterval */
1170  /* with nrmax-th largest error estimate (to be bisected next). */
1171 
1172  L30:
1173  dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
1174  if (errsum <= errbnd) {
1175  goto L115;
1176  }
1177  if (*ier != 0) {
1178  goto L100;
1179  }
1180  if (*last == 2) {
1181  goto L80;
1182  }
1183  if (noext) {
1184  goto L90;
1185  }
1186  erlarg -= erlast;
1187  if ((d__1 = b1 - a1, fabs(d__1)) > small) {
1188  erlarg += erro12;
1189  }
1190  if (extrap) {
1191  goto L40;
1192  }
1193 
1194  /* test whether the interval to be bisected next is the */
1195  /* smallest interval. */
1196 
1197  if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
1198  goto L90;
1199  }
1200  extrap = true;
1201  nrmax = 2;
1202  L40:
1203  if (ierro == 3 || erlarg <= ertest) {
1204  goto L60;
1205  }
1206 
1207  /* the smallest interval has the largest error. */
1208  /* before bisecting decrease the sum of the errors over the */
1209  /* larger intervals (erlarg) and perform extrapolation. */
1210 
1211  id = nrmax;
1212  jupbnd = *last;
1213  if (*last > *limit / 2 + 2) {
1214  jupbnd = *limit + 3 - *last;
1215  }
1216  i__2 = jupbnd;
1217  for (k = id; k <= i__2; ++k) {
1218  maxerr = iord[nrmax];
1219  errmax = elist[maxerr];
1220  if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
1221  goto L90;
1222  }
1223  ++nrmax;
1224  /* L50: */
1225  }
1226 
1227  /* perform extrapolation. */
1228 
1229  L60:
1230  ++numrl2;
1231  rlist2[numrl2 - 1] = area;
1232  dqelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
1233  ++ktmin;
1234  if (ktmin > 5 && *abserr < errsum * .001) {
1235  *ier = 5;
1236  }
1237  if (abseps >= *abserr) {
1238  goto L70;
1239  }
1240  ktmin = 0;
1241  *abserr = abseps;
1242  *result = reseps;
1243  correc = erlarg;
1244  /* Computing MAX */
1245  d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
1246  ertest = max(d__1,d__2);
1247  if (*abserr <= ertest) {
1248  goto L100;
1249  }
1250 
1251  /* prepare bisection of the smallest interval. */
1252 
1253  L70:
1254  if (numrl2 == 1) {
1255  noext = true;
1256  }
1257  if (*ier == 5) {
1258  goto L100;
1259  }
1260  maxerr = iord[1];
1261  errmax = elist[maxerr];
1262  nrmax = 1;
1263  extrap = false;
1264  small *= .5;
1265  erlarg = errsum;
1266  goto L90;
1267  L80:
1268  small = .375;
1269  erlarg = errsum;
1270  ertest = errbnd;
1271  rlist2[1] = area;
1272  L90:
1273  ;
1274  }
1275 
1276  /* set final result and error estimate. */
1277  /* ------------------------------------ */
1278 
1279 L100:
1280  if (*abserr == oflow) {
1281  goto L115;
1282  }
1283  if (*ier + ierro == 0) {
1284  goto L110;
1285  }
1286  if (ierro == 3) {
1287  *abserr += correc;
1288  }
1289  if (*ier == 0) {
1290  *ier = 3;
1291  }
1292  if (*result != 0. && area != 0.) {
1293  goto L105;
1294  }
1295  if (*abserr > errsum) {
1296  goto L115;
1297  }
1298  if (area == 0.) {
1299  goto L130;
1300  }
1301  goto L110;
1302 L105:
1303  if (*abserr / fabs(*result) > errsum / fabs(area)) {
1304  goto L115;
1305  }
1306 
1307  /* test on divergence */
1308 
1309 L110:
1310  /* Computing MAX */
1311  d__1 = fabs(*result), d__2 = fabs(area);
1312  if (ksgn == -1 && max(d__1,d__2) <= defabs * .01) {
1313  goto L130;
1314  }
1315  if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
1316  *ier = 6;
1317  }
1318  goto L130;
1319 
1320  /* compute global integral sum. */
1321 
1322 L115:
1323  *result = 0.;
1324  i__1 = *last;
1325  for (k = 1; k <= i__1; ++k) {
1326  *result += rlist[k];
1327  /* L120: */
1328  }
1329  *abserr = errsum;
1330 L130:
1331  *neval = *last * 30 - 15;
1332  if (*inf == 2) {
1333  *neval <<= 1;
1334  }
1335  if (*ier > 2) {
1336  --(*ier);
1337  }
1338 L999:
1339  return;
1340 } /* dqagie_ */
1341 
1342 void dqagi_(const D_fp& f, const double *bound, const long *inf,
1343  const double *epsabs, const double *epsrel, double *result,
1344  double *abserr, long *neval, long *ier, long *limit,
1345  const long *lenw, long *last, long *iwork, double *work)
1346 {
1347  long l1, l2, l3, lvl;
1348 
1349  /* ***begin prologue dqagi */
1350  /* ***date written 800101 (yymmdd) */
1351  /* ***revision date 830518 (yymmdd) */
1352  /* ***category no. h2a3a1,h2a4a1 */
1353  /* ***keywords automatic integrator, infinite intervals, */
1354  /* general-purpose, transformation, extrapolation, */
1355  /* globally adaptive */
1356  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
1357  /* de doncker,elise,appl. math. & progr. div. -k.u.leuven */
1358  /* ***purpose the routine calculates an approximation result to a given */
1359  /* integral i = integral of f over (bound,+infinity) */
1360  /* or i = integral of f over (-infinity,bound) */
1361  /* or i = integral of f over (-infinity,+infinity) */
1362  /* hopefully satisfying following claim for accuracy */
1363  /* fabs(i-result).le.max(epsabs,epsrel*fabs(i)). */
1364  /* ***description */
1365 
1366  /* integration over infinite intervals */
1367  /* standard fortran subroutine */
1368 
1369  /* parameters */
1370  /* on entry */
1371  /* f - double precision */
1372  /* function subprogram defining the integrand */
1373  /* function f(x). the actual name for f needs to be */
1374  /* declared e x t e r n a l in the driver program. */
1375 
1376  /* bound - double precision */
1377  /* finite bound of integration range */
1378  /* (has no meaning if interval is doubly-infinite) */
1379 
1380  /* inf - long */
1381  /* indicating the kind of integration range involved */
1382  /* inf = 1 corresponds to (bound,+infinity), */
1383  /* inf = -1 to (-infinity,bound), */
1384  /* inf = 2 to (-infinity,+infinity). */
1385 
1386  /* epsabs - double precision */
1387  /* absolute accuracy requested */
1388  /* epsrel - double precision */
1389  /* relative accuracy requested */
1390  /* if epsabs.le.0 */
1391  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
1392  /* the routine will end with ier = 6. */
1393 
1394 
1395  /* on return */
1396  /* result - double precision */
1397  /* approximation to the integral */
1398 
1399  /* abserr - double precision */
1400  /* estimate of the modulus of the absolute error, */
1401  /* which should equal or exceed fabs(i-result) */
1402 
1403  /* neval - long */
1404  /* number of integrand evaluations */
1405 
1406  /* ier - long */
1407  /* ier = 0 normal and reliable termination of the */
1408  /* routine. it is assumed that the requested */
1409  /* accuracy has been achieved. */
1410  /* - ier.gt.0 abnormal termination of the routine. the */
1411  /* estimates for result and error are less */
1412  /* reliable. it is assumed that the requested */
1413  /* accuracy has not been achieved. */
1414  /* error messages */
1415  /* ier = 1 maximum number of subdivisions allowed */
1416  /* has been achieved. one can allow more */
1417  /* subdivisions by increasing the value of */
1418  /* limit (and taking the according dimension */
1419  /* adjustments into account). however, if */
1420  /* this yields no improvement it is advised */
1421  /* to analyze the integrand in order to */
1422  /* determine the integration difficulties. if */
1423  /* the position of a local difficulty can be */
1424  /* determined (e.g. singularity, */
1425  /* discontinuity within the interval) one */
1426  /* will probably gain from splitting up the */
1427  /* interval at this point and calling the */
1428  /* integrator on the subranges. if possible, */
1429  /* an appropriate special-purpose integrator */
1430  /* should be used, which is designed for */
1431  /* handling the type of difficulty involved. */
1432  /* = 2 the occurrence of roundoff error is */
1433  /* detected, which prevents the requested */
1434  /* tolerance from being achieved. */
1435  /* the error may be under-estimated. */
1436  /* = 3 extremely bad integrand behaviour occurs */
1437  /* at some points of the integration */
1438  /* interval. */
1439  /* = 4 the algorithm does not converge. */
1440  /* roundoff error is detected in the */
1441  /* extrapolation table. */
1442  /* it is assumed that the requested tolerance */
1443  /* cannot be achieved, and that the returned */
1444  /* result is the best which can be obtained. */
1445  /* = 5 the integral is probably divergent, or */
1446  /* slowly convergent. it must be noted that */
1447  /* divergence can occur with any other value */
1448  /* of ier. */
1449  /* = 6 the input is invalid, because */
1450  /* (epsabs.le.0 and */
1451  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
1452  /* or limit.lt.1 or leniw.lt.limit*4. */
1453  /* result, abserr, neval, last are set to */
1454  /* zero. exept when limit or leniw is */
1455  /* invalid, iwork(1), work(limit*2+1) and */
1456  /* work(limit*3+1) are set to zero, work(1) */
1457  /* is set to a and work(limit+1) to b. */
1458 
1459  /* dimensioning parameters */
1460  /* limit - long */
1461  /* dimensioning parameter for iwork */
1462  /* limit determines the maximum number of subintervals */
1463  /* in the partition of the given integration interval */
1464  /* (a,b), limit.ge.1. */
1465  /* if limit.lt.1, the routine will end with ier = 6. */
1466 
1467  /* lenw - long */
1468  /* dimensioning parameter for work */
1469  /* lenw must be at least limit*4. */
1470  /* if lenw.lt.limit*4, the routine will end */
1471  /* with ier = 6. */
1472 
1473  /* last - long */
1474  /* on return, last equals the number of subintervals */
1475  /* produced in the subdivision process, which */
1476  /* determines the number of significant elements */
1477  /* actually in the work arrays. */
1478 
1479  /* work arrays */
1480  /* iwork - long */
1481  /* vector of dimension at least limit, the first */
1482  /* k elements of which contain pointers */
1483  /* to the error estimates over the subintervals, */
1484  /* such that work(limit*3+iwork(1)),... , */
1485  /* work(limit*3+iwork(k)) form a decreasing */
1486  /* sequence, with k = last if last.le.(limit/2+2), and */
1487  /* k = limit+1-last otherwise */
1488 
1489  /* work - double precision */
1490  /* vector of dimension at least lenw */
1491  /* on return */
1492  /* work(1), ..., work(last) contain the left */
1493  /* end points of the subintervals in the */
1494  /* partition of (a,b), */
1495  /* work(limit+1), ..., work(limit+last) contain */
1496  /* the right end points, */
1497  /* work(limit*2+1), ...,work(limit*2+last) contain the */
1498  /* integral approximations over the subintervals, */
1499  /* work(limit*3+1), ..., work(limit*3) */
1500  /* contain the error estimates. */
1501  /* ***references (none) */
1502  /* ***routines called dqagie,xerror */
1503  /* ***end prologue dqagi */
1504 
1505 
1506 
1507 
1508  /* check validity of limit and lenw. */
1509 
1510  /* ***first executable statement dqagi */
1511  /* Parameter adjustments */
1512  --iwork;
1513  --work;
1514 
1515  /* Function Body */
1516  *ier = 6;
1517  *neval = 0;
1518  *last = 0;
1519  *result = 0.;
1520  *abserr = 0.;
1521  if (*limit < 1 || *lenw < *limit << 2) {
1522  goto L10;
1523  }
1524 
1525  /* prepare call for dqagie. */
1526 
1527  l1 = *limit + 1;
1528  l2 = *limit + l1;
1529  l3 = *limit + l2;
1530 
1531  dqagie_(f, bound, inf, epsabs, epsrel, limit, result, abserr, neval,
1532  ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
1533 
1534  /* call error handler if necessary. */
1535 
1536  lvl = 0;
1537 L10:
1538  if (*ier == 6) {
1539  lvl = 1;
1540  }
1541  if (*ier != 0) {
1542  xerror_("abnormal return from dqagi", &c__26, ier, &lvl, 26);
1543  }
1544  return;
1545 } /* dqagi_ */
1546 
1547 void dqagpe_(const D_fp& f, const double *a, const double *b, const long *npts2,
1548  const double *points, const double *epsabs, const double *epsrel,
1549  const long *limit, double *result, double *abserr, long *
1550  neval, long *ier, double *alist__, double *blist,
1551  double *rlist, double *elist, double *pts, long *iord,
1552  long *level, long *ndin, long *last)
1553 {
1554  /* System generated locals */
1555  int i__1, i__2;
1556  double d__1, d__2;
1557 
1558  /* Local variables */
1559  long i__, j, k=0;
1560  double a1, a2, b1, b2;
1561  long id, ip1, ind1, ind2;
1562  double area;
1563  double resa, dres, sign;
1564  long ksgn;
1565  double temp;
1566  long nres, nint, jlow, npts;
1567  double area1, area2, area12;
1568  double erro12;
1569  long ierro;
1570  double defab1, defab2;
1571  long ktmin, nrmax;
1572  double oflow, uflow;
1573  bool noext;
1574  long iroff1, iroff2, iroff3;
1575  double res3la[3];
1576  long nintp1;
1577  double error1, error2, rlist2[52];
1578  long numrl2;
1579  double defabs, epmach, erlarg, abseps, correc=0., errbnd, resabs;
1580  long jupbnd;
1581  double erlast;
1582  long levmax;
1583  double errmax;
1584  long maxerr, levcur;
1585  double reseps;
1586  bool extrap;
1587  double ertest, errsum;
1588 
1589  /* ***begin prologue dqagpe */
1590  /* ***date written 800101 (yymmdd) */
1591  /* ***revision date 830518 (yymmdd) */
1592  /* ***category no. h2a2a1 */
1593  /* ***keywords automatic integrator, general-purpose, */
1594  /* singularities at user specified points, */
1595  /* extrapolation, globally adaptive. */
1596  /* ***author piessens,robert ,appl. math. & progr. div. - k.u.leuven */
1597  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
1598  /* ***purpose the routine calculates an approximation result to a given */
1599  /* definite integral i = integral of f over (a,b), hopefully */
1600  /* satisfying following claim for accuracy fabs(i-result).le. */
1601  /* max(epsabs,epsrel*fabs(i)). break points of the integration */
1602  /* interval, where local difficulties of the integrand may */
1603  /* occur(e.g. singularities,discontinuities),provided by user. */
1604  /* ***description */
1605 
1606  /* computation of a definite integral */
1607  /* standard fortran subroutine */
1608  /* double precision version */
1609 
1610  /* parameters */
1611  /* on entry */
1612  /* f - double precision */
1613  /* function subprogram defining the integrand */
1614  /* function f(x). the actual name for f needs to be */
1615  /* declared e x t e r n a l in the driver program. */
1616 
1617  /* a - double precision */
1618  /* lower limit of integration */
1619 
1620  /* b - double precision */
1621  /* upper limit of integration */
1622 
1623  /* npts2 - long */
1624  /* number equal to two more than the number of */
1625  /* user-supplied break points within the integration */
1626  /* range, npts2.ge.2. */
1627  /* if npts2.lt.2, the routine will end with ier = 6. */
1628 
1629  /* points - double precision */
1630  /* vector of dimension npts2, the first (npts2-2) */
1631  /* elements of which are the user provided break */
1632  /* points. if these points do not constitute an */
1633  /* ascending sequence there will be an automatic */
1634  /* sorting. */
1635 
1636  /* epsabs - double precision */
1637  /* absolute accuracy requested */
1638  /* epsrel - double precision */
1639  /* relative accuracy requested */
1640  /* if epsabs.le.0 */
1641  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
1642  /* the routine will end with ier = 6. */
1643 
1644  /* limit - long */
1645  /* gives an upper bound on the number of subintervals */
1646  /* in the partition of (a,b), limit.ge.npts2 */
1647  /* if limit.lt.npts2, the routine will end with */
1648  /* ier = 6. */
1649 
1650  /* on return */
1651  /* result - double precision */
1652  /* approximation to the integral */
1653 
1654  /* abserr - double precision */
1655  /* estimate of the modulus of the absolute error, */
1656  /* which should equal or exceed fabs(i-result) */
1657 
1658  /* neval - long */
1659  /* number of integrand evaluations */
1660 
1661  /* ier - long */
1662  /* ier = 0 normal and reliable termination of the */
1663  /* routine. it is assumed that the requested */
1664  /* accuracy has been achieved. */
1665  /* ier.gt.0 abnormal termination of the routine. */
1666  /* the estimates for integral and error are */
1667  /* less reliable. it is assumed that the */
1668  /* requested accuracy has not been achieved. */
1669  /* error messages */
1670  /* ier = 1 maximum number of subdivisions allowed */
1671  /* has been achieved. one can allow more */
1672  /* subdivisions by increasing the value of */
1673  /* limit (and taking the according dimension */
1674  /* adjustments into account). however, if */
1675  /* this yields no improvement it is advised */
1676  /* to analyze the integrand in order to */
1677  /* determine the integration difficulties. if */
1678  /* the position of a local difficulty can be */
1679  /* determined (i.e. singularity, */
1680  /* discontinuity within the interval), it */
1681  /* should be supplied to the routine as an */
1682  /* element of the vector points. if necessary */
1683  /* an appropriate special-purpose integrator */
1684  /* must be used, which is designed for */
1685  /* handling the type of difficulty involved. */
1686  /* = 2 the occurrence of roundoff error is */
1687  /* detected, which prevents the requested */
1688  /* tolerance from being achieved. */
1689  /* the error may be under-estimated. */
1690  /* = 3 extremely bad integrand behaviour occurs */
1691  /* at some points of the integration */
1692  /* interval. */
1693  /* = 4 the algorithm does not converge. */
1694  /* roundoff error is detected in the */
1695  /* extrapolation table. it is presumed that */
1696  /* the requested tolerance cannot be */
1697  /* achieved, and that the returned result is */
1698  /* the best which can be obtained. */
1699  /* = 5 the integral is probably divergent, or */
1700  /* slowly convergent. it must be noted that */
1701  /* divergence can occur with any other value */
1702  /* of ier.gt.0. */
1703  /* = 6 the input is invalid because */
1704  /* npts2.lt.2 or */
1705  /* break points are specified outside */
1706  /* the integration range or */
1707  /* (epsabs.le.0 and */
1708  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
1709  /* or limit.lt.npts2. */
1710  /* result, abserr, neval, last, rlist(1), */
1711  /* and elist(1) are set to zero. alist(1) and */
1712  /* blist(1) are set to a and b respectively. */
1713 
1714  /* alist - double precision */
1715  /* vector of dimension at least limit, the first */
1716  /* last elements of which are the left end points */
1717  /* of the subintervals in the partition of the given */
1718  /* integration range (a,b) */
1719 
1720  /* blist - double precision */
1721  /* vector of dimension at least limit, the first */
1722  /* last elements of which are the right end points */
1723  /* of the subintervals in the partition of the given */
1724  /* integration range (a,b) */
1725 
1726  /* rlist - double precision */
1727  /* vector of dimension at least limit, the first */
1728  /* last elements of which are the integral */
1729  /* approximations on the subintervals */
1730 
1731  /* elist - double precision */
1732  /* vector of dimension at least limit, the first */
1733  /* last elements of which are the moduli of the */
1734  /* absolute error estimates on the subintervals */
1735 
1736  /* pts - double precision */
1737  /* vector of dimension at least npts2, containing the */
1738  /* integration limits and the break points of the */
1739  /* interval in ascending sequence. */
1740 
1741  /* level - long */
1742  /* vector of dimension at least limit, containing the */
1743  /* subdivision levels of the subinterval, i.e. if */
1744  /* (aa,bb) is a subinterval of (p1,p2) where p1 as */
1745  /* well as p2 is a user-provided break point or */
1746  /* integration limit, then (aa,bb) has level l if */
1747  /* fabs(bb-aa) = fabs(p2-p1)*2**(-l). */
1748 
1749  /* ndin - long */
1750  /* vector of dimension at least npts2, after first */
1751  /* integration over the intervals (pts(i)),pts(i+1), */
1752  /* i = 0,1, ..., npts2-2, the error estimates over */
1753  /* some of the intervals may have been increased */
1754  /* artificially, in order to put their subdivision */
1755  /* forward. if this happens for the subinterval */
1756  /* numbered k, ndin(k) is put to 1, otherwise */
1757  /* ndin(k) = 0. */
1758 
1759  /* iord - long */
1760  /* vector of dimension at least limit, the first k */
1761  /* elements of which are pointers to the */
1762  /* error estimates over the subintervals, */
1763  /* such that elist(iord(1)), ..., elist(iord(k)) */
1764  /* form a decreasing sequence, with k = last */
1765  /* if last.le.(limit/2+2), and k = limit+1-last */
1766  /* otherwise */
1767 
1768  /* last - long */
1769  /* number of subintervals actually produced in the */
1770  /* subdivisions process */
1771 
1772  /* ***references (none) */
1773  /* ***routines called d1mach,dqelg,dqk21,dqpsrt */
1774  /* ***end prologue dqagpe */
1775 
1776 
1777 
1778 
1779  /* the dimension of rlist2 is determined by the value of */
1780  /* limexp in subroutine epsalg (rlist2 should be of dimension */
1781  /* (limexp+2) at least). */
1782 
1783 
1784  /* list of major variables */
1785  /* ----------------------- */
1786 
1787  /* alist - list of left end points of all subintervals */
1788  /* considered up to now */
1789  /* blist - list of right end points of all subintervals */
1790  /* considered up to now */
1791  /* rlist(i) - approximation to the integral over */
1792  /* (alist(i),blist(i)) */
1793  /* rlist2 - array of dimension at least limexp+2 */
1794  /* containing the part of the epsilon table which */
1795  /* is still needed for further computations */
1796  /* elist(i) - error estimate applying to rlist(i) */
1797  /* maxerr - pointer to the interval with largest error */
1798  /* estimate */
1799  /* errmax - elist(maxerr) */
1800  /* erlast - error on the interval currently subdivided */
1801  /* (before that subdivision has taken place) */
1802  /* area - sum of the integrals over the subintervals */
1803  /* errsum - sum of the errors over the subintervals */
1804  /* errbnd - requested accuracy max(epsabs,epsrel* */
1805  /* fabs(result)) */
1806  /* *****1 - variable for the left subinterval */
1807  /* *****2 - variable for the right subinterval */
1808  /* last - index for subdivision */
1809  /* nres - number of calls to the extrapolation routine */
1810  /* numrl2 - number of elements in rlist2. if an appropriate */
1811  /* approximation to the compounded integral has */
1812  /* been obtained, it is put in rlist2(numrl2) after */
1813  /* numrl2 has been increased by one. */
1814  /* erlarg - sum of the errors over the intervals larger */
1815  /* than the smallest interval considered up to now */
1816  /* extrap - bool variable denoting that the routine */
1817  /* is attempting to perform extrapolation. i.e. */
1818  /* before subdividing the smallest interval we */
1819  /* try to decrease the value of erlarg. */
1820  /* noext - bool variable denoting that extrapolation is */
1821  /* no longer allowed (true-value) */
1822 
1823  /* machine dependent constants */
1824  /* --------------------------- */
1825 
1826  /* epmach is the largest relative spacing. */
1827  /* uflow is the smallest positive magnitude. */
1828  /* oflow is the largest positive magnitude. */
1829 
1830  /* ***first executable statement dqagpe */
1831  /* Parameter adjustments */
1832  --ndin;
1833  --pts;
1834  --points;
1835  --level;
1836  --iord;
1837  --elist;
1838  --rlist;
1839  --blist;
1840  --alist__;
1841 
1842  /* Function Body */
1843  epmach = d1mach(c__4);
1844 
1845  /* test on validity of parameters */
1846  /* ----------------------------- */
1847 
1848  *ier = 0;
1849  *neval = 0;
1850  *last = 0;
1851  *result = 0.;
1852  *abserr = 0.;
1853  alist__[1] = *a;
1854  blist[1] = *b;
1855  rlist[1] = 0.;
1856  elist[1] = 0.;
1857  iord[1] = 0;
1858  level[1] = 0;
1859  npts = *npts2 - 2;
1860  /* Computing MAX */
1861  d__1 = epmach * 50.;
1862  if (*npts2 < 2 || *limit <= npts ||
1863  ( *epsabs <= 0. && *epsrel < max(d__1,5e-29 ))) {
1864  *ier = 6;
1865  }
1866  if (*ier == 6) {
1867  goto L999;
1868  }
1869 
1870  /* if any break points are provided, sort them into an */
1871  /* ascending sequence. */
1872 
1873  sign = 1.;
1874  if (*a > *b) {
1875  sign = -1.;
1876  }
1877  pts[1] = min(*a,*b);
1878  if (npts == 0) {
1879  goto L15;
1880  }
1881  i__1 = npts;
1882  for (i__ = 1; i__ <= i__1; ++i__) {
1883  pts[i__ + 1] = points[i__];
1884  /* L10: */
1885  }
1886 L15:
1887  pts[npts + 2] = max(*a,*b);
1888  nint = npts + 1;
1889  a1 = pts[1];
1890  if (npts == 0) {
1891  goto L40;
1892  }
1893  nintp1 = nint + 1;
1894  i__1 = nint;
1895  for (i__ = 1; i__ <= i__1; ++i__) {
1896  ip1 = i__ + 1;
1897  i__2 = nintp1;
1898  for (j = ip1; j <= i__2; ++j) {
1899  if (pts[i__] <= pts[j]) {
1900  goto L20;
1901  }
1902  temp = pts[i__];
1903  pts[i__] = pts[j];
1904  pts[j] = temp;
1905  L20:
1906  ;
1907  }
1908  }
1909  if (pts[1] != min(*a,*b) || pts[nintp1] != max(*a,*b)) {
1910  *ier = 6;
1911  }
1912  if (*ier == 6) {
1913  goto L999;
1914  }
1915 
1916  /* compute first integral and error approximations. */
1917  /* ------------------------------------------------ */
1918 
1919 L40:
1920  resabs = 0.;
1921  i__2 = nint;
1922  for (i__ = 1; i__ <= i__2; ++i__) {
1923  b1 = pts[i__ + 1];
1924  dqk21_(f, &a1, &b1, &area1, &error1, &defabs, &resa);
1925  *abserr += error1;
1926  *result += area1;
1927  ndin[i__] = 0;
1928  if (error1 == resa && error1 != 0.) {
1929  ndin[i__] = 1;
1930  }
1931  resabs += defabs;
1932  level[i__] = 0;
1933  elist[i__] = error1;
1934  alist__[i__] = a1;
1935  blist[i__] = b1;
1936  rlist[i__] = area1;
1937  iord[i__] = i__;
1938  a1 = b1;
1939  /* L50: */
1940  }
1941  errsum = 0.;
1942  i__2 = nint;
1943  for (i__ = 1; i__ <= i__2; ++i__) {
1944  if (ndin[i__] == 1) {
1945  elist[i__] = *abserr;
1946  }
1947  errsum += elist[i__];
1948  /* L55: */
1949  }
1950 
1951  /* test on accuracy. */
1952 
1953  *last = nint;
1954  *neval = nint * 21;
1955  dres = fabs(*result);
1956  /* Computing MAX */
1957  d__1 = *epsabs, d__2 = *epsrel * dres;
1958  errbnd = max(d__1,d__2);
1959  if (*abserr <= epmach * 100. * resabs && *abserr > errbnd) {
1960  *ier = 2;
1961  }
1962  if (nint == 1) {
1963  goto L80;
1964  }
1965  i__2 = npts;
1966  for (i__ = 1; i__ <= i__2; ++i__) {
1967  jlow = i__ + 1;
1968  ind1 = iord[i__];
1969  i__1 = nint;
1970  for (j = jlow; j <= i__1; ++j) {
1971  ind2 = iord[j];
1972  if (elist[ind1] > elist[ind2]) {
1973  goto L60;
1974  }
1975  ind1 = ind2;
1976  k = j;
1977  L60:
1978  ;
1979  }
1980  if (ind1 == iord[i__]) {
1981  goto L70;
1982  }
1983  iord[k] = iord[i__];
1984  iord[i__] = ind1;
1985  L70:
1986  ;
1987  }
1988  if (*limit < *npts2) {
1989  *ier = 1;
1990  }
1991 L80:
1992  if (*ier != 0 || *abserr <= errbnd) {
1993  goto L210;
1994  }
1995 
1996  /* initialization */
1997  /* -------------- */
1998 
1999  rlist2[0] = *result;
2000  maxerr = iord[1];
2001  errmax = elist[maxerr];
2002  area = *result;
2003  nrmax = 1;
2004  nres = 0;
2005  numrl2 = 1;
2006  ktmin = 0;
2007  extrap = false;
2008  noext = false;
2009  erlarg = errsum;
2010  ertest = errbnd;
2011  levmax = 1;
2012  iroff1 = 0;
2013  iroff2 = 0;
2014  iroff3 = 0;
2015  ierro = 0;
2016  uflow = d1mach(c__1);
2017  oflow = d1mach(c__2);
2018  *abserr = oflow;
2019  ksgn = -1;
2020  if (dres >= (1. - epmach * 50.) * resabs) {
2021  ksgn = 1;
2022  }
2023 
2024  /* main do-loop */
2025  /* ------------ */
2026 
2027  i__2 = *limit;
2028  for (*last = *npts2; *last <= i__2; ++(*last)) {
2029 
2030  /* bisect the subinterval with the nrmax-th largest error */
2031  /* estimate. */
2032 
2033  levcur = level[maxerr] + 1;
2034  a1 = alist__[maxerr];
2035  b1 = (alist__[maxerr] + blist[maxerr]) * .5;
2036  a2 = b1;
2037  b2 = blist[maxerr];
2038  erlast = errmax;
2039  dqk21_(f, &a1, &b1, &area1, &error1, &resa, &defab1);
2040  dqk21_(f, &a2, &b2, &area2, &error2, &resa, &defab2);
2041 
2042  /* improve previous approximations to integral */
2043  /* and error and test for accuracy. */
2044 
2045  *neval += 42;
2046  area12 = area1 + area2;
2047  erro12 = error1 + error2;
2048  errsum = errsum + erro12 - errmax;
2049  area = area + area12 - rlist[maxerr];
2050  if (defab1 == error1 || defab2 == error2) {
2051  goto L95;
2052  }
2053  if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
2054  erro12 < errmax * .99) {
2055  goto L90;
2056  }
2057  if (extrap) {
2058  ++iroff2;
2059  }
2060  if (! extrap) {
2061  ++iroff1;
2062  }
2063  L90:
2064  if (*last > 10 && erro12 > errmax) {
2065  ++iroff3;
2066  }
2067  L95:
2068  level[maxerr] = levcur;
2069  level[*last] = levcur;
2070  rlist[maxerr] = area1;
2071  rlist[*last] = area2;
2072  /* Computing MAX */
2073  d__1 = *epsabs, d__2 = *epsrel * fabs(area);
2074  errbnd = max(d__1,d__2);
2075 
2076  /* test for roundoff error and eventually set error flag. */
2077 
2078  if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
2079  *ier = 2;
2080  }
2081  if (iroff2 >= 5) {
2082  ierro = 3;
2083  }
2084 
2085  /* set error flag in the case that the number of */
2086  /* subintervals equals limit. */
2087 
2088  if (*last == *limit) {
2089  *ier = 1;
2090  }
2091 
2092  /* set error flag in the case of bad integrand behaviour */
2093  /* at a point of the integration range */
2094 
2095  /* Computing MAX */
2096  d__1 = fabs(a1), d__2 = fabs(b2);
2097  if (max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
2098  {
2099  *ier = 4;
2100  }
2101 
2102  /* append the newly-created intervals to the list. */
2103 
2104  if (error2 > error1) {
2105  goto L100;
2106  }
2107  alist__[*last] = a2;
2108  blist[maxerr] = b1;
2109  blist[*last] = b2;
2110  elist[maxerr] = error1;
2111  elist[*last] = error2;
2112  goto L110;
2113  L100:
2114  alist__[maxerr] = a2;
2115  alist__[*last] = a1;
2116  blist[*last] = b1;
2117  rlist[maxerr] = area2;
2118  rlist[*last] = area1;
2119  elist[maxerr] = error2;
2120  elist[*last] = error1;
2121 
2122  /* call subroutine dqpsrt to maintain the descending ordering */
2123  /* in the list of error estimates and select the subinterval */
2124  /* with nrmax-th largest error estimate (to be bisected next). */
2125 
2126  L110:
2127  dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
2128  /* ***jump out of do-loop */
2129  if (errsum <= errbnd) {
2130  goto L190;
2131  }
2132  /* ***jump out of do-loop */
2133  if (*ier != 0) {
2134  goto L170;
2135  }
2136  if (noext) {
2137  goto L160;
2138  }
2139  erlarg -= erlast;
2140  if (levcur + 1 <= levmax) {
2141  erlarg += erro12;
2142  }
2143  if (extrap) {
2144  goto L120;
2145  }
2146 
2147  /* test whether the interval to be bisected next is the */
2148  /* smallest interval. */
2149 
2150  if (level[maxerr] + 1 <= levmax) {
2151  goto L160;
2152  }
2153  extrap = true;
2154  nrmax = 2;
2155  L120:
2156  if (ierro == 3 || erlarg <= ertest) {
2157  goto L140;
2158  }
2159 
2160  /* the smallest interval has the largest error. */
2161  /* before bisecting decrease the sum of the errors over */
2162  /* the larger intervals (erlarg) and perform extrapolation. */
2163 
2164  id = nrmax;
2165  jupbnd = *last;
2166  if (*last > *limit / 2 + 2) {
2167  jupbnd = *limit + 3 - *last;
2168  }
2169  i__1 = jupbnd;
2170  for (k = id; k <= i__1; ++k) {
2171  maxerr = iord[nrmax];
2172  errmax = elist[maxerr];
2173  /* ***jump out of do-loop */
2174  if (level[maxerr] + 1 <= levmax) {
2175  goto L160;
2176  }
2177  ++nrmax;
2178  /* L130: */
2179  }
2180 
2181  /* perform extrapolation. */
2182 
2183  L140:
2184  ++numrl2;
2185  rlist2[numrl2 - 1] = area;
2186  if (numrl2 <= 2) {
2187  goto L155;
2188  }
2189  dqelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
2190  ++ktmin;
2191  if (ktmin > 5 && *abserr < errsum * .001) {
2192  *ier = 5;
2193  }
2194  if (abseps >= *abserr) {
2195  goto L150;
2196  }
2197  ktmin = 0;
2198  *abserr = abseps;
2199  *result = reseps;
2200  correc = erlarg;
2201  /* Computing MAX */
2202  d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
2203  ertest = max(d__1,d__2);
2204  /* ***jump out of do-loop */
2205  if (*abserr < ertest) {
2206  goto L170;
2207  }
2208 
2209  /* prepare bisection of the smallest interval. */
2210 
2211  L150:
2212  if (numrl2 == 1) {
2213  noext = true;
2214  }
2215  if (*ier >= 5) {
2216  goto L170;
2217  }
2218  L155:
2219  maxerr = iord[1];
2220  errmax = elist[maxerr];
2221  nrmax = 1;
2222  extrap = false;
2223  ++levmax;
2224  erlarg = errsum;
2225  L160:
2226  ;
2227  }
2228 
2229  /* set the final result. */
2230  /* --------------------- */
2231 
2232 
2233 L170:
2234  if (*abserr == oflow) {
2235  goto L190;
2236  }
2237  if (*ier + ierro == 0) {
2238  goto L180;
2239  }
2240  if (ierro == 3) {
2241  *abserr += correc;
2242  }
2243  if (*ier == 0) {
2244  *ier = 3;
2245  }
2246  if (*result != 0. && area != 0.) {
2247  goto L175;
2248  }
2249  if (*abserr > errsum) {
2250  goto L190;
2251  }
2252  if (area == 0.) {
2253  goto L210;
2254  }
2255  goto L180;
2256 L175:
2257  if (*abserr / fabs(*result) > errsum / fabs(area)) {
2258  goto L190;
2259  }
2260 
2261  /* test on divergence. */
2262 
2263 L180:
2264  /* Computing MAX */
2265  d__1 = fabs(*result), d__2 = fabs(area);
2266  if (ksgn == -1 && max(d__1,d__2) <= resabs * .01) {
2267  goto L210;
2268  }
2269  if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
2270  *ier = 6;
2271  }
2272  goto L210;
2273 
2274  /* compute global integral sum. */
2275 
2276 L190:
2277  *result = 0.;
2278  i__2 = *last;
2279  for (k = 1; k <= i__2; ++k) {
2280  *result += rlist[k];
2281  /* L200: */
2282  }
2283  *abserr = errsum;
2284 L210:
2285  if (*ier > 2) {
2286  --(*ier);
2287  }
2288  *result *= sign;
2289 L999:
2290  return;
2291 } /* dqagpe_ */
2292 
2293 void dqagp_(const D_fp& f, const double *a, const double *b, const long *npts2,
2294  const double *points, const double *epsabs, const double *epsrel,
2295  double *result, double *abserr, long *neval, long *ier,
2296  const long *leniw, const long *lenw, long *last, long *iwork,
2297  double *work)
2298 {
2299  long l1, l2, l3, l4, lvl, limit;
2300 
2301  /* ***begin prologue dqagp */
2302  /* ***date written 800101 (yymmdd) */
2303  /* ***revision date 830518 (yymmdd) */
2304  /* ***category no. h2a2a1 */
2305  /* ***keywords automatic integrator, general-purpose, */
2306  /* singularities at user specified points, */
2307  /* extrapolation, globally adaptive */
2308  /* ***author piessens,robert,appl. math. & progr. div - k.u.leuven */
2309  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
2310  /* ***purpose the routine calculates an approximation result to a given */
2311  /* definite integral i = integral of f over (a,b), */
2312  /* hopefully satisfying following claim for accuracy */
2313  /* break points of the integration interval, where local */
2314  /* difficulties of the integrand may occur (e.g. */
2315  /* singularities, discontinuities), are provided by the user. */
2316  /* ***description */
2317 
2318  /* computation of a definite integral */
2319  /* standard fortran subroutine */
2320  /* double precision version */
2321 
2322  /* parameters */
2323  /* on entry */
2324  /* f - double precision */
2325  /* function subprogram defining the integrand */
2326  /* function f(x). the actual name for f needs to be */
2327  /* declared e x t e r n a l in the driver program. */
2328 
2329  /* a - double precision */
2330  /* lower limit of integration */
2331 
2332  /* b - double precision */
2333  /* upper limit of integration */
2334 
2335  /* npts2 - long */
2336  /* number equal to two more than the number of */
2337  /* user-supplied break points within the integration */
2338  /* range, npts.ge.2. */
2339  /* if npts2.lt.2, the routine will end with ier = 6. */
2340 
2341  /* points - double precision */
2342  /* vector of dimension npts2, the first (npts2-2) */
2343  /* elements of which are the user provided break */
2344  /* points. if these points do not constitute an */
2345  /* ascending sequence there will be an automatic */
2346  /* sorting. */
2347 
2348  /* epsabs - double precision */
2349  /* absolute accuracy requested */
2350  /* epsrel - double precision */
2351  /* relative accuracy requested */
2352  /* if epsabs.le.0 */
2353  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
2354  /* the routine will end with ier = 6. */
2355 
2356  /* on return */
2357  /* result - double precision */
2358  /* approximation to the integral */
2359 
2360  /* abserr - double precision */
2361  /* estimate of the modulus of the absolute error, */
2362  /* which should equal or exceed fabs(i-result) */
2363 
2364  /* neval - long */
2365  /* number of integrand evaluations */
2366 
2367  /* ier - long */
2368  /* ier = 0 normal and reliable termination of the */
2369  /* routine. it is assumed that the requested */
2370  /* accuracy has been achieved. */
2371  /* ier.gt.0 abnormal termination of the routine. */
2372  /* the estimates for integral and error are */
2373  /* less reliable. it is assumed that the */
2374  /* requested accuracy has not been achieved. */
2375  /* error messages */
2376  /* ier = 1 maximum number of subdivisions allowed */
2377  /* has been achieved. one can allow more */
2378  /* subdivisions by increasing the value of */
2379  /* limit (and taking the according dimension */
2380  /* adjustments into account). however, if */
2381  /* this yields no improvement it is advised */
2382  /* to analyze the integrand in order to */
2383  /* determine the integration difficulties. if */
2384  /* the position of a local difficulty can be */
2385  /* determined (i.e. singularity, */
2386  /* discontinuity within the interval), it */
2387  /* should be supplied to the routine as an */
2388  /* element of the vector points. if necessary */
2389  /* an appropriate special-purpose integrator */
2390  /* must be used, which is designed for */
2391  /* handling the type of difficulty involved. */
2392  /* = 2 the occurrence of roundoff error is */
2393  /* detected, which prevents the requested */
2394  /* tolerance from being achieved. */
2395  /* the error may be under-estimated. */
2396  /* = 3 extremely bad integrand behaviour occurs */
2397  /* at some points of the integration */
2398  /* interval. */
2399  /* = 4 the algorithm does not converge. */
2400  /* roundoff error is detected in the */
2401  /* extrapolation table. */
2402  /* it is presumed that the requested */
2403  /* tolerance cannot be achieved, and that */
2404  /* the returned result is the best which */
2405  /* can be obtained. */
2406  /* = 5 the integral is probably divergent, or */
2407  /* slowly convergent. it must be noted that */
2408  /* divergence can occur with any other value */
2409  /* of ier.gt.0. */
2410  /* = 6 the input is invalid because */
2411  /* npts2.lt.2 or */
2412  /* break points are specified outside */
2413  /* the integration range or */
2414  /* (epsabs.le.0 and */
2415  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
2416  /* result, abserr, neval, last are set to */
2417  /* zero. exept when leniw or lenw or npts2 is */
2418  /* invalid, iwork(1), iwork(limit+1), */
2419  /* work(limit*2+1) and work(limit*3+1) */
2420  /* are set to zero. */
2421  /* work(1) is set to a and work(limit+1) */
2422  /* to b (where limit = (leniw-npts2)/2). */
2423 
2424  /* dimensioning parameters */
2425  /* leniw - long */
2426  /* dimensioning parameter for iwork */
2427  /* leniw determines limit = (leniw-npts2)/2, */
2428  /* which is the maximum number of subintervals in the */
2429  /* partition of the given integration interval (a,b), */
2430  /* leniw.ge.(3*npts2-2). */
2431  /* if leniw.lt.(3*npts2-2), the routine will end with */
2432  /* ier = 6. */
2433 
2434  /* lenw - long */
2435  /* dimensioning parameter for work */
2436  /* lenw must be at least leniw*2-npts2. */
2437  /* if lenw.lt.leniw*2-npts2, the routine will end */
2438  /* with ier = 6. */
2439 
2440  /* last - long */
2441  /* on return, last equals the number of subintervals */
2442  /* produced in the subdivision process, which */
2443  /* determines the number of significant elements */
2444  /* actually in the work arrays. */
2445 
2446  /* work arrays */
2447  /* iwork - long */
2448  /* vector of dimension at least leniw. on return, */
2449  /* the first k elements of which contain */
2450  /* pointers to the error estimates over the */
2451  /* subintervals, such that work(limit*3+iwork(1)),..., */
2452  /* work(limit*3+iwork(k)) form a decreasing */
2453  /* sequence, with k = last if last.le.(limit/2+2), and */
2454  /* k = limit+1-last otherwise */
2455  /* iwork(limit+1), ...,iwork(limit+last) contain the */
2456  /* subdivision levels of the subintervals, i.e. */
2457  /* if (aa,bb) is a subinterval of (p1,p2) */
2458  /* where p1 as well as p2 is a user-provided */
2459  /* break point or integration limit, then (aa,bb) has */
2460  /* level l if fabs(bb-aa) = fabs(p2-p1)*2**(-l), */
2461  /* iwork(limit*2+1), ..., iwork(limit*2+npts2) have */
2462  /* no significance for the user, */
2463  /* note that limit = (leniw-npts2)/2. */
2464 
2465  /* work - double precision */
2466  /* vector of dimension at least lenw */
2467  /* on return */
2468  /* work(1), ..., work(last) contain the left */
2469  /* end points of the subintervals in the */
2470  /* partition of (a,b), */
2471  /* work(limit+1), ..., work(limit+last) contain */
2472  /* the right end points, */
2473  /* work(limit*2+1), ..., work(limit*2+last) contain */
2474  /* the integral approximations over the subintervals, */
2475  /* work(limit*3+1), ..., work(limit*3+last) */
2476  /* contain the corresponding error estimates, */
2477  /* work(limit*4+1), ..., work(limit*4+npts2) */
2478  /* contain the integration limits and the */
2479  /* break points sorted in an ascending sequence. */
2480  /* note that limit = (leniw-npts2)/2. */
2481 
2482  /* ***references (none) */
2483  /* ***routines called dqagpe,xerror */
2484  /* ***end prologue dqagp */
2485 
2486 
2487 
2488 
2489  /* check validity of limit and lenw. */
2490 
2491  /* ***first executable statement dqagp */
2492  /* Parameter adjustments */
2493  --points;
2494  --iwork;
2495  --work;
2496 
2497  /* Function Body */
2498  *ier = 6;
2499  *neval = 0;
2500  *last = 0;
2501  *result = 0.;
2502  *abserr = 0.;
2503  if (*leniw < *npts2 * 3 - 2 || *lenw < (*leniw << 1) - *npts2 || *npts2 <
2504  2) {
2505  goto L10;
2506  }
2507 
2508  /* prepare call for dqagpe. */
2509 
2510  limit = (*leniw - *npts2) / 2;
2511  l1 = limit + 1;
2512  l2 = limit + l1;
2513  l3 = limit + l2;
2514  l4 = limit + l3;
2515 
2516  dqagpe_(f, a, b, npts2, &points[1], epsabs, epsrel, &limit, result,
2517  abserr, neval, ier, &work[1], &work[l1], &work[l2], &work[l3], &
2518  work[l4], &iwork[1], &iwork[l1], &iwork[l2], last);
2519 
2520  /* call error handler if necessary. */
2521 
2522  lvl = 0;
2523 L10:
2524  if (*ier == 6) {
2525  lvl = 1;
2526  }
2527  if (*ier != 0) {
2528  xerror_("abnormal return from dqagp", &c__26, ier, &lvl, 26);
2529  }
2530  return;
2531 } /* dqagp_ */
2532 
2533 void dqagse_(const D_fp& f, const double *a, const double *b,
2534  const double *epsabs, const double *epsrel, const long *limit,
2535  double *result,
2536  double *abserr, long *neval, long *ier, double *alist__,
2537  double *blist, double *rlist, double *elist, long *
2538  iord, long *last)
2539 {
2540  /* System generated locals */
2541  int i__1, i__2;
2542  double d__1, d__2;
2543 
2544  /* Local variables */
2545  long k;
2546  double a1, a2, b1, b2;
2547  long id;
2548  double area;
2549  double dres;
2550  long ksgn, nres;
2551  double area1, area2, area12;
2552  double small=0., erro12;
2553  long ierro;
2554  double defab1, defab2;
2555  long ktmin, nrmax;
2556  double oflow, uflow;
2557  bool noext;
2558  long iroff1, iroff2, iroff3;
2559  double res3la[3], error1, error2, rlist2[52];
2560  long numrl2;
2561  double defabs, epmach, erlarg=0., abseps, correc=0., errbnd, resabs;
2562  long jupbnd;
2563  double erlast, errmax;
2564  long maxerr;
2565  double reseps;
2566  bool extrap;
2567  double ertest=0., errsum;
2568 
2569  /* ***begin prologue dqagse */
2570  /* ***date written 800101 (yymmdd) */
2571  /* ***revision date 830518 (yymmdd) */
2572  /* ***category no. h2a1a1 */
2573  /* ***keywords automatic integrator, general-purpose, */
2574  /* (end point) singularities, extrapolation, */
2575  /* globally adaptive */
2576  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
2577  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
2578  /* ***purpose the routine calculates an approximation result to a given */
2579  /* definite integral i = integral of f over (a,b), */
2580  /* hopefully satisfying following claim for accuracy */
2581  /* fabs(i-result).le.max(epsabs,epsrel*fabs(i)). */
2582  /* ***description */
2583 
2584  /* computation of a definite integral */
2585  /* standard fortran subroutine */
2586  /* double precision version */
2587 
2588  /* parameters */
2589  /* on entry */
2590  /* f - double precision */
2591  /* function subprogram defining the integrand */
2592  /* function f(x). the actual name for f needs to be */
2593  /* declared e x t e r n a l in the driver program. */
2594 
2595  /* a - double precision */
2596  /* lower limit of integration */
2597 
2598  /* b - double precision */
2599  /* upper limit of integration */
2600 
2601  /* epsabs - double precision */
2602  /* absolute accuracy requested */
2603  /* epsrel - double precision */
2604  /* relative accuracy requested */
2605  /* if epsabs.le.0 */
2606  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
2607  /* the routine will end with ier = 6. */
2608 
2609  /* limit - long */
2610  /* gives an upperbound on the number of subintervals */
2611  /* in the partition of (a,b) */
2612 
2613  /* on return */
2614  /* result - double precision */
2615  /* approximation to the integral */
2616 
2617  /* abserr - double precision */
2618  /* estimate of the modulus of the absolute error, */
2619  /* which should equal or exceed fabs(i-result) */
2620 
2621  /* neval - long */
2622  /* number of integrand evaluations */
2623 
2624  /* ier - long */
2625  /* ier = 0 normal and reliable termination of the */
2626  /* routine. it is assumed that the requested */
2627  /* accuracy has been achieved. */
2628  /* ier.gt.0 abnormal termination of the routine */
2629  /* the estimates for integral and error are */
2630  /* less reliable. it is assumed that the */
2631  /* requested accuracy has not been achieved. */
2632  /* error messages */
2633  /* = 1 maximum number of subdivisions allowed */
2634  /* has been achieved. one can allow more sub- */
2635  /* divisions by increasing the value of limit */
2636  /* (and taking the according dimension */
2637  /* adjustments into account). however, if */
2638  /* this yields no improvement it is advised */
2639  /* to analyze the integrand in order to */
2640  /* determine the integration difficulties. if */
2641  /* the position of a local difficulty can be */
2642  /* determined (e.g. singularity, */
2643  /* discontinuity within the interval) one */
2644  /* will probably gain from splitting up the */
2645  /* interval at this point and calling the */
2646  /* integrator on the subranges. if possible, */
2647  /* an appropriate special-purpose integrator */
2648  /* should be used, which is designed for */
2649  /* handling the type of difficulty involved. */
2650  /* = 2 the occurrence of roundoff error is detec- */
2651  /* ted, which prevents the requested */
2652  /* tolerance from being achieved. */
2653  /* the error may be under-estimated. */
2654  /* = 3 extremely bad integrand behaviour */
2655  /* occurs at some points of the integration */
2656  /* interval. */
2657  /* = 4 the algorithm does not converge. */
2658  /* roundoff error is detected in the */
2659  /* extrapolation table. */
2660  /* it is presumed that the requested */
2661  /* tolerance cannot be achieved, and that the */
2662  /* returned result is the best which can be */
2663  /* obtained. */
2664  /* = 5 the integral is probably divergent, or */
2665  /* slowly convergent. it must be noted that */
2666  /* divergence can occur with any other value */
2667  /* of ier. */
2668  /* = 6 the input is invalid, because */
2669  /* epsabs.le.0 and */
2670  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28). */
2671  /* result, abserr, neval, last, rlist(1), */
2672  /* iord(1) and elist(1) are set to zero. */
2673  /* alist(1) and blist(1) are set to a and b */
2674  /* respectively. */
2675 
2676  /* alist - double precision */
2677  /* vector of dimension at least limit, the first */
2678  /* last elements of which are the left end points */
2679  /* of the subintervals in the partition of the */
2680  /* given integration range (a,b) */
2681 
2682  /* blist - double precision */
2683  /* vector of dimension at least limit, the first */
2684  /* last elements of which are the right end points */
2685  /* of the subintervals in the partition of the given */
2686  /* integration range (a,b) */
2687 
2688  /* rlist - double precision */
2689  /* vector of dimension at least limit, the first */
2690  /* last elements of which are the integral */
2691  /* approximations on the subintervals */
2692 
2693  /* elist - double precision */
2694  /* vector of dimension at least limit, the first */
2695  /* last elements of which are the moduli of the */
2696  /* absolute error estimates on the subintervals */
2697 
2698  /* iord - long */
2699  /* vector of dimension at least limit, the first k */
2700  /* elements of which are pointers to the */
2701  /* error estimates over the subintervals, */
2702  /* such that elist(iord(1)), ..., elist(iord(k)) */
2703  /* form a decreasing sequence, with k = last */
2704  /* if last.le.(limit/2+2), and k = limit+1-last */
2705  /* otherwise */
2706 
2707  /* last - long */
2708  /* number of subintervals actually produced in the */
2709  /* subdivision process */
2710 
2711  /* ***references (none) */
2712  /* ***routines called d1mach,dqelg,dqk21,dqpsrt */
2713  /* ***end prologue dqagse */
2714 
2715 
2716 
2717 
2718  /* the dimension of rlist2 is determined by the value of */
2719  /* limexp in subroutine dqelg (rlist2 should be of dimension */
2720  /* (limexp+2) at least). */
2721 
2722  /* list of major variables */
2723  /* ----------------------- */
2724 
2725  /* alist - list of left end points of all subintervals */
2726  /* considered up to now */
2727  /* blist - list of right end points of all subintervals */
2728  /* considered up to now */
2729  /* rlist(i) - approximation to the integral over */
2730  /* (alist(i),blist(i)) */
2731  /* rlist2 - array of dimension at least limexp+2 containing */
2732  /* the part of the epsilon table which is still */
2733  /* needed for further computations */
2734  /* elist(i) - error estimate applying to rlist(i) */
2735  /* maxerr - pointer to the interval with largest error */
2736  /* estimate */
2737  /* errmax - elist(maxerr) */
2738  /* erlast - error on the interval currently subdivided */
2739  /* (before that subdivision has taken place) */
2740  /* area - sum of the integrals over the subintervals */
2741  /* errsum - sum of the errors over the subintervals */
2742  /* errbnd - requested accuracy max(epsabs,epsrel* */
2743  /* fabs(result)) */
2744  /* *****1 - variable for the left interval */
2745  /* *****2 - variable for the right interval */
2746  /* last - index for subdivision */
2747  /* nres - number of calls to the extrapolation routine */
2748  /* numrl2 - number of elements currently in rlist2. if an */
2749  /* appropriate approximation to the compounded */
2750  /* integral has been obtained it is put in */
2751  /* rlist2(numrl2) after numrl2 has been increased */
2752  /* by one. */
2753  /* small - length of the smallest interval considered up */
2754  /* to now, multiplied by 1.5 */
2755  /* erlarg - sum of the errors over the intervals larger */
2756  /* than the smallest interval considered up to now */
2757  /* extrap - bool variable denoting that the routine is */
2758  /* attempting to perform extrapolation i.e. before */
2759  /* subdividing the smallest interval we try to */
2760  /* decrease the value of erlarg. */
2761  /* noext - bool variable denoting that extrapolation */
2762  /* is no longer allowed (true value) */
2763 
2764  /* machine dependent constants */
2765  /* --------------------------- */
2766 
2767  /* epmach is the largest relative spacing. */
2768  /* uflow is the smallest positive magnitude. */
2769  /* oflow is the largest positive magnitude. */
2770 
2771  /* ***first executable statement dqagse */
2772  /* Parameter adjustments */
2773  --iord;
2774  --elist;
2775  --rlist;
2776  --blist;
2777  --alist__;
2778 
2779  /* Function Body */
2780  epmach = d1mach(c__4);
2781 
2782  /* test on validity of parameters */
2783  /* ------------------------------ */
2784  *ier = 0;
2785  *neval = 0;
2786  *last = 0;
2787  *result = 0.;
2788  *abserr = 0.;
2789  alist__[1] = *a;
2790  blist[1] = *b;
2791  rlist[1] = 0.;
2792  elist[1] = 0.;
2793  /* Computing MAX */
2794  d__1 = epmach * 50.;
2795  if (*epsabs <= 0. && *epsrel < max(d__1,5e-29)) {
2796  *ier = 6;
2797  }
2798  if (*ier == 6) {
2799  goto L999;
2800  }
2801 
2802  /* first approximation to the integral */
2803  /* ----------------------------------- */
2804 
2805  uflow = d1mach(c__1);
2806  oflow = d1mach(c__2);
2807  ierro = 0;
2808  dqk21_(f, a, b, result, abserr, &defabs, &resabs);
2809 
2810  /* test on accuracy. */
2811 
2812  dres = fabs(*result);
2813  /* Computing MAX */
2814  d__1 = *epsabs, d__2 = *epsrel * dres;
2815  errbnd = max(d__1,d__2);
2816  *last = 1;
2817  rlist[1] = *result;
2818  elist[1] = *abserr;
2819  iord[1] = 1;
2820  if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) {
2821  *ier = 2;
2822  }
2823  if (*limit == 1) {
2824  *ier = 1;
2825  }
2826  if (*ier != 0 || ( *abserr <= errbnd && *abserr != resabs ) || *abserr == 0.)
2827  {
2828  goto L140;
2829  }
2830 
2831  /* initialization */
2832  /* -------------- */
2833 
2834  rlist2[0] = *result;
2835  errmax = *abserr;
2836  maxerr = 1;
2837  area = *result;
2838  errsum = *abserr;
2839  *abserr = oflow;
2840  nrmax = 1;
2841  nres = 0;
2842  numrl2 = 2;
2843  ktmin = 0;
2844  extrap = false;
2845  noext = false;
2846  iroff1 = 0;
2847  iroff2 = 0;
2848  iroff3 = 0;
2849  ksgn = -1;
2850  if (dres >= (1. - epmach * 50.) * defabs) {
2851  ksgn = 1;
2852  }
2853 
2854  /* main do-loop */
2855  /* ------------ */
2856 
2857  i__1 = *limit;
2858  for (*last = 2; *last <= i__1; ++(*last)) {
2859 
2860  /* bisect the subinterval with the nrmax-th largest error */
2861  /* estimate. */
2862 
2863  a1 = alist__[maxerr];
2864  b1 = (alist__[maxerr] + blist[maxerr]) * .5;
2865  a2 = b1;
2866  b2 = blist[maxerr];
2867  erlast = errmax;
2868  dqk21_(f, &a1, &b1, &area1, &error1, &resabs, &defab1);
2869  dqk21_(f, &a2, &b2, &area2, &error2, &resabs, &defab2);
2870 
2871  /* improve previous approximations to integral */
2872  /* and error and test for accuracy. */
2873 
2874  area12 = area1 + area2;
2875  erro12 = error1 + error2;
2876  errsum = errsum + erro12 - errmax;
2877  area = area + area12 - rlist[maxerr];
2878  if (defab1 == error1 || defab2 == error2) {
2879  goto L15;
2880  }
2881  if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
2882  erro12 < errmax * .99) {
2883  goto L10;
2884  }
2885  if (extrap) {
2886  ++iroff2;
2887  }
2888  if (! extrap) {
2889  ++iroff1;
2890  }
2891  L10:
2892  if (*last > 10 && erro12 > errmax) {
2893  ++iroff3;
2894  }
2895  L15:
2896  rlist[maxerr] = area1;
2897  rlist[*last] = area2;
2898  /* Computing MAX */
2899  d__1 = *epsabs, d__2 = *epsrel * fabs(area);
2900  errbnd = max(d__1,d__2);
2901 
2902  /* test for roundoff error and eventually set error flag. */
2903 
2904  if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
2905  *ier = 2;
2906  }
2907  if (iroff2 >= 5) {
2908  ierro = 3;
2909  }
2910 
2911  /* set error flag in the case that the number of subintervals */
2912  /* equals limit. */
2913 
2914  if (*last == *limit) {
2915  *ier = 1;
2916  }
2917 
2918  /* set error flag in the case of bad integrand behaviour */
2919  /* at a point of the integration range. */
2920 
2921  /* Computing MAX */
2922  d__1 = fabs(a1), d__2 = fabs(b2);
2923  if (max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
2924  {
2925  *ier = 4;
2926  }
2927 
2928  /* append the newly-created intervals to the list. */
2929 
2930  if (error2 > error1) {
2931  goto L20;
2932  }
2933  alist__[*last] = a2;
2934  blist[maxerr] = b1;
2935  blist[*last] = b2;
2936  elist[maxerr] = error1;
2937  elist[*last] = error2;
2938  goto L30;
2939  L20:
2940  alist__[maxerr] = a2;
2941  alist__[*last] = a1;
2942  blist[*last] = b1;
2943  rlist[maxerr] = area2;
2944  rlist[*last] = area1;
2945  elist[maxerr] = error2;
2946  elist[*last] = error1;
2947 
2948  /* call subroutine dqpsrt to maintain the descending ordering */
2949  /* in the list of error estimates and select the subinterval */
2950  /* with nrmax-th largest error estimate (to be bisected next). */
2951 
2952  L30:
2953  dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
2954  /* ***jump out of do-loop */
2955  if (errsum <= errbnd) {
2956  goto L115;
2957  }
2958  /* ***jump out of do-loop */
2959  if (*ier != 0) {
2960  goto L100;
2961  }
2962  if (*last == 2) {
2963  goto L80;
2964  }
2965  if (noext) {
2966  goto L90;
2967  }
2968  erlarg -= erlast;
2969  if ((d__1 = b1 - a1, fabs(d__1)) > small) {
2970  erlarg += erro12;
2971  }
2972  if (extrap) {
2973  goto L40;
2974  }
2975 
2976  /* test whether the interval to be bisected next is the */
2977  /* smallest interval. */
2978 
2979  if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
2980  goto L90;
2981  }
2982  extrap = true;
2983  nrmax = 2;
2984  L40:
2985  if (ierro == 3 || erlarg <= ertest) {
2986  goto L60;
2987  }
2988 
2989  /* the smallest interval has the largest error. */
2990  /* before bisecting decrease the sum of the errors over the */
2991  /* larger intervals (erlarg) and perform extrapolation. */
2992 
2993  id = nrmax;
2994  jupbnd = *last;
2995  if (*last > *limit / 2 + 2) {
2996  jupbnd = *limit + 3 - *last;
2997  }
2998  i__2 = jupbnd;
2999  for (k = id; k <= i__2; ++k) {
3000  maxerr = iord[nrmax];
3001  errmax = elist[maxerr];
3002  /* ***jump out of do-loop */
3003  if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
3004  goto L90;
3005  }
3006  ++nrmax;
3007  /* L50: */
3008  }
3009 
3010  /* perform extrapolation. */
3011 
3012  L60:
3013  ++numrl2;
3014  rlist2[numrl2 - 1] = area;
3015  dqelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
3016  ++ktmin;
3017  if (ktmin > 5 && *abserr < errsum * .001) {
3018  *ier = 5;
3019  }
3020  if (abseps >= *abserr) {
3021  goto L70;
3022  }
3023  ktmin = 0;
3024  *abserr = abseps;
3025  *result = reseps;
3026  correc = erlarg;
3027  /* Computing MAX */
3028  d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
3029  ertest = max(d__1,d__2);
3030  /* ***jump out of do-loop */
3031  if (*abserr <= ertest) {
3032  goto L100;
3033  }
3034 
3035  /* prepare bisection of the smallest interval. */
3036 
3037  L70:
3038  if (numrl2 == 1) {
3039  noext = true;
3040  }
3041  if (*ier == 5) {
3042  goto L100;
3043  }
3044  maxerr = iord[1];
3045  errmax = elist[maxerr];
3046  nrmax = 1;
3047  extrap = false;
3048  small *= .5;
3049  erlarg = errsum;
3050  goto L90;
3051  L80:
3052  small = (d__1 = *b - *a, fabs(d__1)) * .375;
3053  erlarg = errsum;
3054  ertest = errbnd;
3055  rlist2[1] = area;
3056  L90:
3057  ;
3058  }
3059 
3060  /* set final result and error estimate. */
3061  /* ------------------------------------ */
3062 
3063 L100:
3064  if (*abserr == oflow) {
3065  goto L115;
3066  }
3067  if (*ier + ierro == 0) {
3068  goto L110;
3069  }
3070  if (ierro == 3) {
3071  *abserr += correc;
3072  }
3073  if (*ier == 0) {
3074  *ier = 3;
3075  }
3076  if (*result != 0. && area != 0.) {
3077  goto L105;
3078  }
3079  if (*abserr > errsum) {
3080  goto L115;
3081  }
3082  if (area == 0.) {
3083  goto L130;
3084  }
3085  goto L110;
3086 L105:
3087  if (*abserr / fabs(*result) > errsum / fabs(area)) {
3088  goto L115;
3089  }
3090 
3091  /* test on divergence. */
3092 
3093 L110:
3094  /* Computing MAX */
3095  d__1 = fabs(*result), d__2 = fabs(area);
3096  if (ksgn == -1 && max(d__1,d__2) <= defabs * .01) {
3097  goto L130;
3098  }
3099  if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
3100  *ier = 6;
3101  }
3102  goto L130;
3103 
3104  /* compute global integral sum. */
3105 
3106 L115:
3107  *result = 0.;
3108  i__1 = *last;
3109  for (k = 1; k <= i__1; ++k) {
3110  *result += rlist[k];
3111  /* L120: */
3112  }
3113  *abserr = errsum;
3114 L130:
3115  if (*ier > 2) {
3116  --(*ier);
3117  }
3118 L140:
3119  *neval = *last * 42 - 21;
3120 L999:
3121  return;
3122 } /* dqagse_ */
3123 
3124 void dqags_(const D_fp& f, const double *a, const double *b, const double *epsabs,
3125  const double *epsrel, double *result, double *abserr,
3126  long *neval, long *ier, const long *limit, const long *lenw, long *
3127  last, long *iwork, double *work)
3128 {
3129  long l1, l2, l3, lvl;
3130 
3131  /* ***begin prologue dqags */
3132  /* ***date written 800101 (yymmdd) */
3133  /* ***revision date 830518 (yymmdd) */
3134  /* ***category no. h2a1a1 */
3135  /* ***keywords automatic integrator, general-purpose, */
3136  /* (end-point) singularities, extrapolation, */
3137  /* globally adaptive */
3138  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
3139  /* de doncker,elise,appl. math. & prog. div. - k.u.leuven */
3140  /* ***purpose the routine calculates an approximation result to a given */
3141  /* definite integral i = integral of f over (a,b), */
3142  /* hopefully satisfying following claim for accuracy */
3143  /* fabs(i-result).le.max(epsabs,epsrel*fabs(i)). */
3144  /* ***description */
3145 
3146  /* computation of a definite integral */
3147  /* standard fortran subroutine */
3148  /* double precision version */
3149 
3150 
3151  /* parameters */
3152  /* on entry */
3153  /* f - double precision */
3154  /* function subprogram defining the integrand */
3155  /* function f(x). the actual name for f needs to be */
3156  /* declared e x t e r n a l in the driver program. */
3157 
3158  /* a - double precision */
3159  /* lower limit of integration */
3160 
3161  /* b - double precision */
3162  /* upper limit of integration */
3163 
3164  /* epsabs - double precision */
3165  /* absolute accuracy requested */
3166  /* epsrel - double precision */
3167  /* relative accuracy requested */
3168  /* if epsabs.le.0 */
3169  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
3170  /* the routine will end with ier = 6. */
3171 
3172  /* on return */
3173  /* result - double precision */
3174  /* approximation to the integral */
3175 
3176  /* abserr - double precision */
3177  /* estimate of the modulus of the absolute error, */
3178  /* which should equal or exceed fabs(i-result) */
3179 
3180  /* neval - long */
3181  /* number of integrand evaluations */
3182 
3183  /* ier - long */
3184  /* ier = 0 normal and reliable termination of the */
3185  /* routine. it is assumed that the requested */
3186  /* accuracy has been achieved. */
3187  /* ier.gt.0 abnormal termination of the routine */
3188  /* the estimates for integral and error are */
3189  /* less reliable. it is assumed that the */
3190  /* requested accuracy has not been achieved. */
3191  /* error messages */
3192  /* ier = 1 maximum number of subdivisions allowed */
3193  /* has been achieved. one can allow more sub- */
3194  /* divisions by increasing the value of limit */
3195  /* (and taking the according dimension */
3196  /* adjustments into account. however, if */
3197  /* this yields no improvement it is advised */
3198  /* to analyze the integrand in order to */
3199  /* determine the integration difficulties. if */
3200  /* the position of a local difficulty can be */
3201  /* determined (e.g. singularity, */
3202  /* discontinuity within the interval) one */
3203  /* will probably gain from splitting up the */
3204  /* interval at this point and calling the */
3205  /* integrator on the subranges. if possible, */
3206  /* an appropriate special-purpose integrator */
3207  /* should be used, which is designed for */
3208  /* handling the type of difficulty involved. */
3209  /* = 2 the occurrence of roundoff error is detec- */
3210  /* ted, which prevents the requested */
3211  /* tolerance from being achieved. */
3212  /* the error may be under-estimated. */
3213  /* = 3 extremely bad integrand behaviour */
3214  /* occurs at some points of the integration */
3215  /* interval. */
3216  /* = 4 the algorithm does not converge. */
3217  /* roundoff error is detected in the */
3218  /* extrapolation table. it is presumed that */
3219  /* the requested tolerance cannot be */
3220  /* achieved, and that the returned result is */
3221  /* the best which can be obtained. */
3222  /* = 5 the integral is probably divergent, or */
3223  /* slowly convergent. it must be noted that */
3224  /* divergence can occur with any other value */
3225  /* of ier. */
3226  /* = 6 the input is invalid, because */
3227  /* (epsabs.le.0 and */
3228  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28) */
3229  /* or limit.lt.1 or lenw.lt.limit*4. */
3230  /* result, abserr, neval, last are set to */
3231  /* zero.except when limit or lenw is invalid, */
3232  /* iwork(1), work(limit*2+1) and */
3233  /* work(limit*3+1) are set to zero, work(1) */
3234  /* is set to a and work(limit+1) to b. */
3235 
3236  /* dimensioning parameters */
3237  /* limit - long */
3238  /* dimensioning parameter for iwork */
3239  /* limit determines the maximum number of subintervals */
3240  /* in the partition of the given integration interval */
3241  /* (a,b), limit.ge.1. */
3242  /* if limit.lt.1, the routine will end with ier = 6. */
3243 
3244  /* lenw - long */
3245  /* dimensioning parameter for work */
3246  /* lenw must be at least limit*4. */
3247  /* if lenw.lt.limit*4, the routine will end */
3248  /* with ier = 6. */
3249 
3250  /* last - long */
3251  /* on return, last equals the number of subintervals */
3252  /* produced in the subdivision process, detemines the */
3253  /* number of significant elements actually in the work */
3254  /* arrays. */
3255 
3256  /* work arrays */
3257  /* iwork - long */
3258  /* vector of dimension at least limit, the first k */
3259  /* elements of which contain pointers */
3260  /* to the error estimates over the subintervals */
3261  /* such that work(limit*3+iwork(1)),... , */
3262  /* work(limit*3+iwork(k)) form a decreasing */
3263  /* sequence, with k = last if last.le.(limit/2+2), */
3264  /* and k = limit+1-last otherwise */
3265 
3266  /* work - double precision */
3267  /* vector of dimension at least lenw */
3268  /* on return */
3269  /* work(1), ..., work(last) contain the left */
3270  /* end-points of the subintervals in the */
3271  /* partition of (a,b), */
3272  /* work(limit+1), ..., work(limit+last) contain */
3273  /* the right end-points, */
3274  /* work(limit*2+1), ..., work(limit*2+last) contain */
3275  /* the integral approximations over the subintervals, */
3276  /* work(limit*3+1), ..., work(limit*3+last) */
3277  /* contain the error estimates. */
3278 
3279  /* ***references (none) */
3280  /* ***routines called dqagse,xerror */
3281  /* ***end prologue dqags */
3282 
3283 
3284 
3285 
3286 
3287  /* check validity of limit and lenw. */
3288 
3289  /* ***first executable statement dqags */
3290  /* Parameter adjustments */
3291  --iwork;
3292  --work;
3293 
3294  /* Function Body */
3295  *ier = 6;
3296  *neval = 0;
3297  *last = 0;
3298  *result = 0.;
3299  *abserr = 0.;
3300  if (*limit < 1 || *lenw < *limit << 2) {
3301  goto L10;
3302  }
3303 
3304  /* prepare call for dqagse. */
3305 
3306  l1 = *limit + 1;
3307  l2 = *limit + l1;
3308  l3 = *limit + l2;
3309 
3310  dqagse_(f, a, b, epsabs, epsrel, limit, result, abserr, neval, ier,
3311  &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
3312 
3313  /* call error handler if necessary. */
3314 
3315  lvl = 0;
3316 L10:
3317  if (*ier == 6) {
3318  lvl = 1;
3319  }
3320  if (*ier != 0) {
3321  xerror_("abnormal return from dqags", &c__26, ier, &lvl, 26);
3322  }
3323  return;
3324 } /* dqags_ */
3325 
3326 void dqawce_(const D_fp& f, const double *a, const double *b, const double *c__,
3327  const double *epsabs, const double *epsrel, const long *limit,
3328  double *result, double *abserr, long *neval, long *ier,
3329  double *alist__, double *blist, double *rlist, double
3330  *elist, long *iord, long *last)
3331 {
3332  /* System generated locals */
3333  int i__1;
3334  double d__1, d__2;
3335 
3336  /* Local variables */
3337  long k;
3338  double a1, a2, b1, b2, aa, bb;
3339  long nev;
3340  double area, area1, area2, area12;
3341  double erro12;
3342  long krule, nrmax;
3343  double uflow;
3344  long iroff1, iroff2;
3345  double error1, error2, epmach, errbnd, errmax;
3346  long maxerr;
3347  double errsum;
3348 
3349  /* ***begin prologue dqawce */
3350  /* ***date written 800101 (yymmdd) */
3351  /* ***revision date 830518 (yymmdd) */
3352  /* ***category no. h2a2a1,j4 */
3353  /* ***keywords automatic integrator, special-purpose, */
3354  /* cauchy principal value, clenshaw-curtis method */
3355  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
3356  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
3357  /* *** purpose the routine calculates an approximation result to a */
3358  /* cauchy principal value i = integral of f*w over (a,b) */
3359  /* (w(x) = 1/(x-c), (c.ne.a, c.ne.b), hopefully satisfying */
3360  /* following claim for accuracy */
3361  /* fabs(i-result).le.max(epsabs,epsrel*fabs(i)) */
3362  /* ***description */
3363 
3364  /* computation of a cauchy principal value */
3365  /* standard fortran subroutine */
3366  /* double precision version */
3367 
3368  /* parameters */
3369  /* on entry */
3370  /* f - double precision */
3371  /* function subprogram defining the integrand */
3372  /* function f(x). the actual name for f needs to be */
3373  /* declared e x t e r n a l in the driver program. */
3374 
3375  /* a - double precision */
3376  /* lower limit of integration */
3377 
3378  /* b - double precision */
3379  /* upper limit of integration */
3380 
3381  /* c - double precision */
3382  /* parameter in the weight function, c.ne.a, c.ne.b */
3383  /* if c = a or c = b, the routine will end with */
3384  /* ier = 6. */
3385 
3386  /* epsabs - double precision */
3387  /* absolute accuracy requested */
3388  /* epsrel - double precision */
3389  /* relative accuracy requested */
3390  /* if epsabs.le.0 */
3391  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
3392  /* the routine will end with ier = 6. */
3393 
3394  /* limit - long */
3395  /* gives an upper bound on the number of subintervals */
3396  /* in the partition of (a,b), limit.ge.1 */
3397 
3398  /* on return */
3399  /* result - double precision */
3400  /* approximation to the integral */
3401 
3402  /* abserr - double precision */
3403  /* estimate of the modulus of the absolute error, */
3404  /* which should equal or exceed fabs(i-result) */
3405 
3406  /* neval - long */
3407  /* number of integrand evaluations */
3408 
3409  /* ier - long */
3410  /* ier = 0 normal and reliable termination of the */
3411  /* routine. it is assumed that the requested */
3412  /* accuracy has been achieved. */
3413  /* ier.gt.0 abnormal termination of the routine */
3414  /* the estimates for integral and error are */
3415  /* less reliable. it is assumed that the */
3416  /* requested accuracy has not been achieved. */
3417  /* error messages */
3418  /* ier = 1 maximum number of subdivisions allowed */
3419  /* has been achieved. one can allow more sub- */
3420  /* divisions by increasing the value of */
3421  /* limit. however, if this yields no */
3422  /* improvement it is advised to analyze the */
3423  /* the integrand, in order to determine the */
3424  /* the integration difficulties. if the */
3425  /* position of a local difficulty can be */
3426  /* determined (e.g. singularity, */
3427  /* discontinuity within the interval) one */
3428  /* will probably gain from splitting up the */
3429  /* interval at this point and calling */
3430  /* appropriate integrators on the subranges. */
3431  /* = 2 the occurrence of roundoff error is detec- */
3432  /* ted, which prevents the requested */
3433  /* tolerance from being achieved. */
3434  /* = 3 extremely bad integrand behaviour */
3435  /* occurs at some interior points of */
3436  /* the integration interval. */
3437  /* = 6 the input is invalid, because */
3438  /* c = a or c = b or */
3439  /* (epsabs.le.0 and */
3440  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
3441  /* or limit.lt.1. */
3442  /* result, abserr, neval, rlist(1), elist(1), */
3443  /* iord(1) and last are set to zero. alist(1) */
3444  /* and blist(1) are set to a and b */
3445  /* respectively. */
3446 
3447  /* alist - double precision */
3448  /* vector of dimension at least limit, the first */
3449  /* last elements of which are the left */
3450  /* end points of the subintervals in the partition */
3451  /* of the given integration range (a,b) */
3452 
3453  /* blist - double precision */
3454  /* vector of dimension at least limit, the first */
3455  /* last elements of which are the right */
3456  /* end points of the subintervals in the partition */
3457  /* of the given integration range (a,b) */
3458 
3459  /* rlist - double precision */
3460  /* vector of dimension at least limit, the first */
3461  /* last elements of which are the integral */
3462  /* approximations on the subintervals */
3463 
3464  /* elist - double precision */
3465  /* vector of dimension limit, the first last */
3466  /* elements of which are the moduli of the absolute */
3467  /* error estimates on the subintervals */
3468 
3469  /* iord - long */
3470  /* vector of dimension at least limit, the first k */
3471  /* elements of which are pointers to the error */
3472  /* estimates over the subintervals, so that */
3473  /* elist(iord(1)), ..., elist(iord(k)) with k = last */
3474  /* if last.le.(limit/2+2), and k = limit+1-last */
3475  /* otherwise, form a decreasing sequence */
3476 
3477  /* last - long */
3478  /* number of subintervals actually produced in */
3479  /* the subdivision process */
3480 
3481  /* ***references (none) */
3482  /* ***routines called d1mach,dqc25c,dqpsrt */
3483  /* ***end prologue dqawce */
3484 
3485 
3486 
3487 
3488  /* list of major variables */
3489  /* ----------------------- */
3490 
3491  /* alist - list of left end points of all subintervals */
3492  /* considered up to now */
3493  /* blist - list of right end points of all subintervals */
3494  /* considered up to now */
3495  /* rlist(i) - approximation to the integral over */
3496  /* (alist(i),blist(i)) */
3497  /* elist(i) - error estimate applying to rlist(i) */
3498  /* maxerr - pointer to the interval with largest */
3499  /* error estimate */
3500  /* errmax - elist(maxerr) */
3501  /* area - sum of the integrals over the subintervals */
3502  /* errsum - sum of the errors over the subintervals */
3503  /* errbnd - requested accuracy max(epsabs,epsrel* */
3504  /* fabs(result)) */
3505  /* *****1 - variable for the left subinterval */
3506  /* *****2 - variable for the right subinterval */
3507  /* last - index for subdivision */
3508 
3509 
3510  /* machine dependent constants */
3511  /* --------------------------- */
3512 
3513  /* epmach is the largest relative spacing. */
3514  /* uflow is the smallest positive magnitude. */
3515 
3516  /* ***first executable statement dqawce */
3517  /* Parameter adjustments */
3518  --iord;
3519  --elist;
3520  --rlist;
3521  --blist;
3522  --alist__;
3523 
3524  /* Function Body */
3525  epmach = d1mach(c__4);
3526  uflow = d1mach(c__1);
3527 
3528 
3529  /* test on validity of parameters */
3530  /* ------------------------------ */
3531 
3532  *ier = 6;
3533  *neval = 0;
3534  *last = 0;
3535  alist__[1] = *a;
3536  blist[1] = *b;
3537  rlist[1] = 0.;
3538  elist[1] = 0.;
3539  iord[1] = 0;
3540  *result = 0.;
3541  *abserr = 0.;
3542  /* Computing MAX */
3543  d__1 = epmach * 50.;
3544  if (*c__ == *a || *c__ == *b || ( *epsabs <= 0. && *epsrel < max(d__1,5e-29) )
3545  ) {
3546  goto L999;
3547  }
3548 
3549  /* first approximation to the integral */
3550  /* ----------------------------------- */
3551 
3552  aa = *a;
3553  bb = *b;
3554  if (*a <= *b) {
3555  goto L10;
3556  }
3557  aa = *b;
3558  bb = *a;
3559 L10:
3560  *ier = 0;
3561  krule = 1;
3562  dqc25c_(f, &aa, &bb, c__, result, abserr, &krule, neval);
3563  *last = 1;
3564  rlist[1] = *result;
3565  elist[1] = *abserr;
3566  iord[1] = 1;
3567  alist__[1] = *a;
3568  blist[1] = *b;
3569 
3570  /* test on accuracy */
3571 
3572  /* Computing MAX */
3573  d__1 = *epsabs, d__2 = *epsrel * fabs(*result);
3574  errbnd = max(d__1,d__2);
3575  if (*limit == 1) {
3576  *ier = 1;
3577  }
3578  /* Computing MIN */
3579  d__1 = fabs(*result) * .01;
3580  if (*abserr < min(d__1,errbnd) || *ier == 1) {
3581  goto L70;
3582  }
3583 
3584  /* initialization */
3585  /* -------------- */
3586 
3587  alist__[1] = aa;
3588  blist[1] = bb;
3589  rlist[1] = *result;
3590  errmax = *abserr;
3591  maxerr = 1;
3592  area = *result;
3593  errsum = *abserr;
3594  nrmax = 1;
3595  iroff1 = 0;
3596  iroff2 = 0;
3597 
3598  /* main do-loop */
3599  /* ------------ */
3600 
3601  i__1 = *limit;
3602  for (*last = 2; *last <= i__1; ++(*last)) {
3603 
3604  /* bisect the subinterval with nrmax-th largest */
3605  /* error estimate. */
3606 
3607  a1 = alist__[maxerr];
3608  b1 = (alist__[maxerr] + blist[maxerr]) * .5;
3609  b2 = blist[maxerr];
3610  if (*c__ <= b1 && *c__ > a1) {
3611  b1 = (*c__ + b2) * .5;
3612  }
3613  if (*c__ > b1 && *c__ < b2) {
3614  b1 = (a1 + *c__) * .5;
3615  }
3616  a2 = b1;
3617  krule = 2;
3618  dqc25c_(f, &a1, &b1, c__, &area1, &error1, &krule, &nev);
3619  *neval += nev;
3620  dqc25c_(f, &a2, &b2, c__, &area2, &error2, &krule, &nev);
3621  *neval += nev;
3622 
3623  /* improve previous approximations to integral */
3624  /* and error and test for accuracy. */
3625 
3626  area12 = area1 + area2;
3627  erro12 = error1 + error2;
3628  errsum = errsum + erro12 - errmax;
3629  area = area + area12 - rlist[maxerr];
3630  if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) < fabs(area12) * 1e-5 &&
3631  erro12 >= errmax * .99 && krule == 0) {
3632  ++iroff1;
3633  }
3634  if (*last > 10 && erro12 > errmax && krule == 0) {
3635  ++iroff2;
3636  }
3637  rlist[maxerr] = area1;
3638  rlist[*last] = area2;
3639  /* Computing MAX */
3640  d__1 = *epsabs, d__2 = *epsrel * fabs(area);
3641  errbnd = max(d__1,d__2);
3642  if (errsum <= errbnd) {
3643  goto L15;
3644  }
3645 
3646  /* test for roundoff error and eventually set error flag. */
3647 
3648  if (iroff1 >= 6 && iroff2 > 20) {
3649  *ier = 2;
3650  }
3651 
3652  /* set error flag in the case that number of interval */
3653  /* bisections exceeds limit. */
3654 
3655  if (*last == *limit) {
3656  *ier = 1;
3657  }
3658 
3659  /* set error flag in the case of bad integrand behaviour */
3660  /* at a point of the integration range. */
3661 
3662  /* Computing MAX */
3663  d__1 = fabs(a1), d__2 = fabs(b2);
3664  if (max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
3665  {
3666  *ier = 3;
3667  }
3668 
3669  /* append the newly-created intervals to the list. */
3670 
3671  L15:
3672  if (error2 > error1) {
3673  goto L20;
3674  }
3675  alist__[*last] = a2;
3676  blist[maxerr] = b1;
3677  blist[*last] = b2;
3678  elist[maxerr] = error1;
3679  elist[*last] = error2;
3680  goto L30;
3681  L20:
3682  alist__[maxerr] = a2;
3683  alist__[*last] = a1;
3684  blist[*last] = b1;
3685  rlist[maxerr] = area2;
3686  rlist[*last] = area1;
3687  elist[maxerr] = error2;
3688  elist[*last] = error1;
3689 
3690  /* call subroutine dqpsrt to maintain the descending ordering */
3691  /* in the list of error estimates and select the subinterval */
3692  /* with nrmax-th largest error estimate (to be bisected next). */
3693 
3694  L30:
3695  dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
3696  /* ***jump out of do-loop */
3697  if (*ier != 0 || errsum <= errbnd) {
3698  goto L50;
3699  }
3700  /* L40: */
3701  }
3702 
3703  /* compute final result. */
3704  /* --------------------- */
3705 
3706 L50:
3707  *result = 0.;
3708  i__1 = *last;
3709  for (k = 1; k <= i__1; ++k) {
3710  *result += rlist[k];
3711  /* L60: */
3712  }
3713  *abserr = errsum;
3714 L70:
3715  if (aa == *b) {
3716  *result = -(*result);
3717  }
3718 L999:
3719  return;
3720 } /* dqawce_ */
3721 
3722 void dqawc_(const D_fp& f, const double *a, const double *b, const double *c__,
3723  const double *epsabs, const double *epsrel, double *result,
3724  double *abserr, long *neval, long *ier, long *limit,
3725  const long *lenw, long *last, long *iwork, double *work)
3726 {
3727  long l1, l2, l3, lvl;
3728 
3729  /* ***begin prologue dqawc */
3730  /* ***date written 800101 (yymmdd) */
3731  /* ***revision date 830518 (yymmdd) */
3732  /* ***category no. h2a2a1,j4 */
3733  /* ***keywords automatic integrator, special-purpose, */
3734  /* cauchy principal value, */
3735  /* clenshaw-curtis, globally adaptive */
3736  /* ***author piessens,robert ,appl. math. & progr. div. - k.u.leuven */
3737  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
3738  /* ***purpose the routine calculates an approximation result to a */
3739  /* cauchy principal value i = integral of f*w over (a,b) */
3740  /* (w(x) = 1/((x-c), c.ne.a, c.ne.b), hopefully satisfying */
3741  /* following claim for accuracy */
3742  /* fabs(i-result).le.max(epsabe,epsrel*fabs(i)). */
3743  /* ***description */
3744 
3745  /* computation of a cauchy principal value */
3746  /* standard fortran subroutine */
3747  /* double precision version */
3748 
3749 
3750  /* parameters */
3751  /* on entry */
3752  /* f - double precision */
3753  /* function subprogram defining the integrand */
3754  /* function f(x). the actual name for f needs to be */
3755  /* declared e x t e r n a l in the driver program. */
3756 
3757  /* a - double precision */
3758  /* under limit of integration */
3759 
3760  /* b - double precision */
3761  /* upper limit of integration */
3762 
3763  /* c - parameter in the weight function, c.ne.a, c.ne.b. */
3764  /* if c = a or c = b, the routine will end with */
3765  /* ier = 6 . */
3766 
3767  /* epsabs - double precision */
3768  /* absolute accuracy requested */
3769  /* epsrel - double precision */
3770  /* relative accuracy requested */
3771  /* if epsabs.le.0 */
3772  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
3773  /* the routine will end with ier = 6. */
3774 
3775  /* on return */
3776  /* result - double precision */
3777  /* approximation to the integral */
3778 
3779  /* abserr - double precision */
3780  /* estimate or the modulus of the absolute error, */
3781  /* which should equal or exceed fabs(i-result) */
3782 
3783  /* neval - long */
3784  /* number of integrand evaluations */
3785 
3786  /* ier - long */
3787  /* ier = 0 normal and reliable termination of the */
3788  /* routine. it is assumed that the requested */
3789  /* accuracy has been achieved. */
3790  /* ier.gt.0 abnormal termination of the routine */
3791  /* the estimates for integral and error are */
3792  /* less reliable. it is assumed that the */
3793  /* requested accuracy has not been achieved. */
3794  /* error messages */
3795  /* ier = 1 maximum number of subdivisions allowed */
3796  /* has been achieved. one can allow more sub- */
3797  /* divisions by increasing the value of limit */
3798  /* (and taking the according dimension */
3799  /* adjustments into account). however, if */
3800  /* this yields no improvement it is advised */
3801  /* to analyze the integrand in order to */
3802  /* determine the integration difficulties. */
3803  /* if the position of a local difficulty */
3804  /* can be determined (e.g. singularity, */
3805  /* discontinuity within the interval) one */
3806  /* will probably gain from splitting up the */
3807  /* interval at this point and calling */
3808  /* appropriate integrators on the subranges. */
3809  /* = 2 the occurrence of roundoff error is detec- */
3810  /* ted, which prevents the requested */
3811  /* tolerance from being achieved. */
3812  /* = 3 extremely bad integrand behaviour occurs */
3813  /* at some points of the integration */
3814  /* interval. */
3815  /* = 6 the input is invalid, because */
3816  /* c = a or c = b or */
3817  /* (epsabs.le.0 and */
3818  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
3819  /* or limit.lt.1 or lenw.lt.limit*4. */
3820  /* result, abserr, neval, last are set to */
3821  /* zero. exept when lenw or limit is invalid, */
3822  /* iwork(1), work(limit*2+1) and */
3823  /* work(limit*3+1) are set to zero, work(1) */
3824  /* is set to a and work(limit+1) to b. */
3825 
3826  /* dimensioning parameters */
3827  /* limit - long */
3828  /* dimensioning parameter for iwork */
3829  /* limit determines the maximum number of subintervals */
3830  /* in the partition of the given integration interval */
3831  /* (a,b), limit.ge.1. */
3832  /* if limit.lt.1, the routine will end with ier = 6. */
3833 
3834  /* lenw - long */
3835  /* dimensioning parameter for work */
3836  /* lenw must be at least limit*4. */
3837  /* if lenw.lt.limit*4, the routine will end with */
3838  /* ier = 6. */
3839 
3840  /* last - long */
3841  /* on return, last equals the number of subintervals */
3842  /* produced in the subdivision process, which */
3843  /* determines the number of significant elements */
3844  /* actually in the work arrays. */
3845 
3846  /* work arrays */
3847  /* iwork - long */
3848  /* vector of dimension at least limit, the first k */
3849  /* elements of which contain pointers */
3850  /* to the error estimates over the subintervals, */
3851  /* such that work(limit*3+iwork(1)), ... , */
3852  /* work(limit*3+iwork(k)) form a decreasing */
3853  /* sequence, with k = last if last.le.(limit/2+2), */
3854  /* and k = limit+1-last otherwise */
3855 
3856  /* work - double precision */
3857  /* vector of dimension at least lenw */
3858  /* on return */
3859  /* work(1), ..., work(last) contain the left */
3860  /* end points of the subintervals in the */
3861  /* partition of (a,b), */
3862  /* work(limit+1), ..., work(limit+last) contain */
3863  /* the right end points, */
3864  /* work(limit*2+1), ..., work(limit*2+last) contain */
3865  /* the integral approximations over the subintervals, */
3866  /* work(limit*3+1), ..., work(limit*3+last) */
3867  /* contain the error estimates. */
3868 
3869  /* ***references (none) */
3870  /* ***routines called dqawce,xerror */
3871  /* ***end prologue dqawc */
3872 
3873 
3874 
3875 
3876  /* check validity of limit and lenw. */
3877 
3878  /* ***first executable statement dqawc */
3879  /* Parameter adjustments */
3880  --iwork;
3881  --work;
3882 
3883  /* Function Body */
3884  *ier = 6;
3885  *neval = 0;
3886  *last = 0;
3887  *result = 0.;
3888  *abserr = 0.;
3889  if (*limit < 1 || *lenw < *limit << 2) {
3890  goto L10;
3891  }
3892 
3893  /* prepare call for dqawce. */
3894 
3895  l1 = *limit + 1;
3896  l2 = *limit + l1;
3897  l3 = *limit + l2;
3898  dqawce_(f, a, b, c__, epsabs, epsrel, limit, result, abserr, neval,
3899  ier, &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
3900 
3901  /* call error handler if necessary. */
3902 
3903  lvl = 0;
3904 L10:
3905  if (*ier == 6) {
3906  lvl = 1;
3907  }
3908  if (*ier != 0) {
3909  xerror_("abnormal return from dqawc", &c__26, ier, &lvl, 26);
3910  }
3911  return;
3912 } /* dqawc_ */
3913 
3914 void dqawfe_(const D_fp& f, const double *a, const double *omega,
3915  const long *integr, const double *epsabs, const long *limlst,
3916  const long *limit, const long *maxp1,
3917  double *result, double *abserr, long *neval,
3918  long *ier, double *rslst, double *erlst, long *ierlst,
3919  long *lst, double *alist__, double *blist,
3920  double *rlist, double *elist, long *iord, long *nnlog,
3921  double *chebmo)
3922 {
3923  /* Initialized data */
3924 
3925  double p = .9;
3926  double pi = 3.1415926535897932384626433832795;
3927 
3928  /* System generated locals */
3929  int chebmo_dim1, chebmo_offset, i__1;
3930  double d__1, d__2;
3931 
3932  /* Local variables */
3933  long l;
3934  double c1, c2, p1, dl, ep;
3935  long ll=0;
3936  double drl=0., eps;
3937  long nev;
3938  double fact, epsa;
3939  long last, nres;
3940  double psum[52];
3941  double cycle;
3942  long ktmin;
3943  double uflow;
3944  double res3la[3];
3945  long numrl2;
3946  double abseps, correc;
3947  long momcom;
3948  double reseps, errsum;
3949 
3950  /* ***begin prologue dqawfe */
3951  /* ***date written 800101 (yymmdd) */
3952  /* ***revision date 830518 (yymmdd) */
3953  /* ***category no. h2a3a1 */
3954  /* ***keywords automatic integrator, special-purpose, */
3955  /* fourier integrals, */
3956  /* integration between zeros with dqawoe, */
3957  /* convergence acceleration with dqelg */
3958  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
3959  /* dedoncker,elise,appl. math. & progr. div. - k.u.leuven */
3960  /* ***purpose the routine calculates an approximation result to a */
3961  /* given fourier integal */
3962  /* i = integral of f(x)*w(x) over (a,infinity) */
3963  /* where w(x)=cos(omega*x) or w(x)=sin(omega*x), */
3964  /* hopefully satisfying following claim for accuracy */
3965  /* fabs(i-result).le.epsabs. */
3966  /* ***description */
3967 
3968  /* computation of fourier integrals */
3969  /* standard fortran subroutine */
3970  /* double precision version */
3971 
3972  /* parameters */
3973  /* on entry */
3974  /* f - double precision */
3975  /* function subprogram defining the integrand */
3976  /* function f(x). the actual name for f needs to */
3977  /* be declared e x t e r n a l in the driver program. */
3978 
3979  /* a - double precision */
3980  /* lower limit of integration */
3981 
3982  /* omega - double precision */
3983  /* parameter in the weight function */
3984 
3985  /* integr - long */
3986  /* indicates which weight function is used */
3987  /* integr = 1 w(x) = cos(omega*x) */
3988  /* integr = 2 w(x) = sin(omega*x) */
3989  /* if integr.ne.1.and.integr.ne.2, the routine will */
3990  /* end with ier = 6. */
3991 
3992  /* epsabs - double precision */
3993  /* absolute accuracy requested, epsabs.gt.0 */
3994  /* if epsabs.le.0, the routine will end with ier = 6. */
3995 
3996  /* limlst - long */
3997  /* limlst gives an upper bound on the number of */
3998  /* cycles, limlst.ge.1. */
3999  /* if limlst.lt.3, the routine will end with ier = 6. */
4000 
4001  /* limit - long */
4002  /* gives an upper bound on the number of subintervals */
4003  /* allowed in the partition of each cycle, limit.ge.1 */
4004  /* each cycle, limit.ge.1. */
4005 
4006  /* maxp1 - long */
4007  /* gives an upper bound on the number of */
4008  /* chebyshev moments which can be stored, i.e. */
4009  /* for the intervals of lengths fabs(b-a)*2**(-l), */
4010  /* l=0,1, ..., maxp1-2, maxp1.ge.1 */
4011 
4012  /* on return */
4013  /* result - double precision */
4014  /* approximation to the integral x */
4015 
4016  /* abserr - double precision */
4017  /* estimate of the modulus of the absolute error, */
4018  /* which should equal or exceed fabs(i-result) */
4019 
4020  /* neval - long */
4021  /* number of integrand evaluations */
4022 
4023  /* ier - ier = 0 normal and reliable termination of */
4024  /* the routine. it is assumed that the */
4025  /* requested accuracy has been achieved. */
4026  /* ier.gt.0 abnormal termination of the routine. the */
4027  /* estimates for integral and error are less */
4028  /* reliable. it is assumed that the requested */
4029  /* accuracy has not been achieved. */
4030  /* error messages */
4031  /* if omega.ne.0 */
4032  /* ier = 1 maximum number of cycles allowed */
4033  /* has been achieved., i.e. of subintervals */
4034  /* (a+(k-1)c,a+kc) where */
4035  /* c = (2*int(fabs(omega))+1)*pi/fabs(omega), */
4036  /* for k = 1, 2, ..., lst. */
4037  /* one can allow more cycles by increasing */
4038  /* the value of limlst (and taking the */
4039  /* according dimension adjustments into */
4040  /* account). */
4041  /* examine the array iwork which contains */
4042  /* the error flags on the cycles, in order to */
4043  /* look for eventual local integration */
4044  /* difficulties. if the position of a local */
4045  /* difficulty can be determined (e.g. */
4046  /* singularity, discontinuity within the */
4047  /* interval) one will probably gain from */
4048  /* splitting up the interval at this polong */
4049  /* and calling appropriate integrators on */
4050  /* the subranges. */
4051  /* = 4 the extrapolation table constructed for */
4052  /* convergence acceleration of the series */
4053  /* formed by the integral contributions over */
4054  /* the cycles, does not converge to within */
4055  /* the requested accuracy. as in the case of */
4056  /* ier = 1, it is advised to examine the */
4057  /* array iwork which contains the error */
4058  /* flags on the cycles. */
4059  /* = 6 the input is invalid because */
4060  /* (integr.ne.1 and integr.ne.2) or */
4061  /* epsabs.le.0 or limlst.lt.3. */
4062  /* result, abserr, neval, lst are set */
4063  /* to zero. */
4064  /* = 7 bad integrand behaviour occurs within one */
4065  /* or more of the cycles. location and type */
4066  /* of the difficulty involved can be */
4067  /* determined from the vector ierlst. here */
4068  /* lst is the number of cycles actually */
4069  /* needed (see below). */
4070  /* ierlst(k) = 1 the maximum number of */
4071  /* subdivisions (= limit) has */
4072  /* been achieved on the k th */
4073  /* cycle. */
4074  /* = 2 occurrence of roundoff error */
4075  /* is detected and prevents the */
4076  /* tolerance imposed on the */
4077  /* k th cycle, from being */
4078  /* achieved. */
4079  /* = 3 extremely bad integrand */
4080  /* behaviour occurs at some */
4081  /* points of the k th cycle. */
4082  /* = 4 the integration procedure */
4083  /* over the k th cycle does */
4084  /* not converge (to within the */
4085  /* required accuracy) due to */
4086  /* roundoff in the */
4087  /* extrapolation procedure */
4088  /* invoked on this cycle. it */
4089  /* is assumed that the result */
4090  /* on this interval is the */
4091  /* best which can be obtained. */
4092  /* = 5 the integral over the k th */
4093  /* cycle is probably divergent */
4094  /* or slowly convergent. it */
4095  /* must be noted that */
4096  /* divergence can occur with */
4097  /* any other value of */
4098  /* ierlst(k). */
4099  /* if omega = 0 and integr = 1, */
4100  /* the integral is calculated by means of dqagie */
4101  /* and ier = ierlst(1) (with meaning as described */
4102  /* for ierlst(k), k = 1). */
4103 
4104  /* rslst - double precision */
4105  /* vector of dimension at least limlst */
4106  /* rslst(k) contains the integral contribution */
4107  /* over the interval (a+(k-1)c,a+kc) where */
4108  /* c = (2*int(fabs(omega))+1)*pi/fabs(omega), */
4109  /* k = 1, 2, ..., lst. */
4110  /* note that, if omega = 0, rslst(1) contains */
4111  /* the value of the integral over (a,infinity). */
4112 
4113  /* erlst - double precision */
4114  /* vector of dimension at least limlst */
4115  /* erlst(k) contains the error estimate corresponding */
4116  /* with rslst(k). */
4117 
4118  /* ierlst - long */
4119  /* vector of dimension at least limlst */
4120  /* ierlst(k) contains the error flag corresponding */
4121  /* with rslst(k). for the meaning of the local error */
4122  /* flags see description of output parameter ier. */
4123 
4124  /* lst - long */
4125  /* number of subintervals needed for the integration */
4126  /* if omega = 0 then lst is set to 1. */
4127 
4128  /* alist, blist, rlist, elist - double precision */
4129  /* vector of dimension at least limit, */
4130 
4131  /* iord, nnlog - long */
4132  /* vector of dimension at least limit, providing */
4133  /* space for the quantities needed in the subdivision */
4134  /* process of each cycle */
4135 
4136  /* chebmo - double precision */
4137  /* array of dimension at least (maxp1,25), providing */
4138  /* space for the chebyshev moments needed within the */
4139  /* cycles */
4140 
4141  /* ***references (none) */
4142  /* ***routines called d1mach,dqagie,dqawoe,dqelg */
4143  /* ***end prologue dqawfe */
4144 
4145 
4146 
4147 
4148 
4149  /* the dimension of psum is determined by the value of */
4150  /* limexp in subroutine dqelg (psum must be of dimension */
4151  /* (limexp+2) at least). */
4152 
4153  /* list of major variables */
4154  /* ----------------------- */
4155 
4156  /* c1, c2 - end points of subinterval (of length cycle) */
4157  /* cycle - (2*int(fabs(omega))+1)*pi/fabs(omega) */
4158  /* psum - vector of dimension at least (limexp+2) */
4159  /* (see routine dqelg) */
4160  /* psum contains the part of the epsilon table */
4161  /* which is still needed for further computations. */
4162  /* each element of psum is a partial sum of the */
4163  /* series which should sum to the value of the */
4164  /* integral. */
4165  /* errsum - sum of error estimates over the subintervals, */
4166  /* calculated cumulatively */
4167  /* epsa - absolute tolerance requested over current */
4168  /* subinterval */
4169  /* chebmo - array containing the modified chebyshev */
4170  /* moments (see also routine dqc25f) */
4171 
4172  /* Parameter adjustments */
4173  --ierlst;
4174  --erlst;
4175  --rslst;
4176  --nnlog;
4177  --iord;
4178  --elist;
4179  --rlist;
4180  --blist;
4181  --alist__;
4182  chebmo_dim1 = *maxp1;
4183  chebmo_offset = 1 + chebmo_dim1;
4184  chebmo -= chebmo_offset;
4185 
4186  /* Function Body */
4187 
4188  /* test on validity of parameters */
4189  /* ------------------------------ */
4190 
4191  /* ***first executable statement dqawfe */
4192  *result = 0.;
4193  *abserr = 0.;
4194  *neval = 0;
4195  *lst = 0;
4196  *ier = 0;
4197  if ( (*integr != 1 && *integr != 2 ) || *epsabs <= 0. || *limlst < 3) {
4198  *ier = 6;
4199  }
4200  if (*ier == 6) {
4201  goto L999;
4202  }
4203  if (*omega != 0.) {
4204  goto L10;
4205  }
4206 
4207  /* integration by dqagie if omega is zero */
4208  /* -------------------------------------- */
4209 
4210  if (*integr == 1) {
4211  dqagie_(f, &c_b20, &c__1, epsabs, &c_b20, limit, result, abserr,
4212  neval, ier, &alist__[1], &blist[1], &rlist[1], &elist[1], &
4213  iord[1], &last);
4214  }
4215  rslst[1] = *result;
4216  erlst[1] = *abserr;
4217  ierlst[1] = *ier;
4218  *lst = 1;
4219  goto L999;
4220 
4221  /* initializations */
4222  /* --------------- */
4223 
4224 L10:
4225  l = (int) fabs(*omega);
4226  dl = (double) ((l << 1) + 1);
4227  cycle = dl * pi / fabs(*omega);
4228  *ier = 0;
4229  ktmin = 0;
4230  *neval = 0;
4231  numrl2 = 0;
4232  nres = 0;
4233  c1 = *a;
4234  c2 = cycle + *a;
4235  p1 = 1. - p;
4236  uflow = d1mach(c__1);
4237  eps = *epsabs;
4238  if (*epsabs > uflow / p1) {
4239  eps = *epsabs * p1;
4240  }
4241  ep = eps;
4242  fact = 1.;
4243  correc = 0.;
4244  *abserr = 0.;
4245  errsum = 0.;
4246 
4247  /* main do-loop */
4248  /* ------------ */
4249 
4250  i__1 = *limlst;
4251  for (*lst = 1; *lst <= i__1; ++(*lst)) {
4252 
4253  /* integrate over current subinterval. */
4254 
4255  epsa = eps * fact;
4256  dqawoe_(f, &c1, &c2, omega, integr, &epsa, &c_b20, limit, lst,
4257  maxp1, &rslst[*lst], &erlst[*lst], &nev, &ierlst[*lst], &last,
4258  &alist__[1], &blist[1], &rlist[1], &elist[1], &iord[1], &
4259  nnlog[1], &momcom, &chebmo[chebmo_offset]);
4260  *neval += nev;
4261  fact *= p;
4262  errsum += erlst[*lst];
4263  drl = (d__1 = rslst[*lst], fabs(d__1)) * 50.;
4264 
4265  /* test on accuracy with partial sum */
4266 
4267  if (errsum + drl <= *epsabs && *lst >= 6) {
4268  goto L80;
4269  }
4270  /* Computing MAX */
4271  d__1 = correc, d__2 = erlst[*lst];
4272  correc = max(d__1,d__2);
4273  if (ierlst[*lst] != 0) {
4274  /* Computing MAX */
4275  d__1 = ep, d__2 = correc * p1;
4276  eps = max(d__1,d__2);
4277  }
4278  if (ierlst[*lst] != 0) {
4279  *ier = 7;
4280  }
4281  if (*ier == 7 && errsum + drl <= correc * 10. && *lst > 5) {
4282  goto L80;
4283  }
4284  ++numrl2;
4285  if (*lst > 1) {
4286  goto L20;
4287  }
4288  psum[0] = rslst[1];
4289  goto L40;
4290  L20:
4291  psum[numrl2 - 1] = psum[ll - 1] + rslst[*lst];
4292  if (*lst == 2) {
4293  goto L40;
4294  }
4295 
4296  /* test on maximum number of subintervals */
4297 
4298  if (*lst == *limlst) {
4299  *ier = 1;
4300  }
4301 
4302  /* perform new extrapolation */
4303 
4304  dqelg_(&numrl2, psum, &reseps, &abseps, res3la, &nres);
4305 
4306  /* test whether extrapolated result is influenced by roundoff */
4307 
4308  ++ktmin;
4309  if (ktmin >= 15 && *abserr <= (errsum + drl) * .001) {
4310  *ier = 4;
4311  }
4312  if (abseps > *abserr && *lst != 3) {
4313  goto L30;
4314  }
4315  *abserr = abseps;
4316  *result = reseps;
4317  ktmin = 0;
4318 
4319  /* if ier is not 0, check whether direct result (partial sum) */
4320  /* or extrapolated result yields the best integral */
4321  /* approximation */
4322 
4323  if (*abserr + correc * 10. <= *epsabs ||
4324  ( *abserr <= *epsabs && correc * 10. >= *epsabs ) ) {
4325  goto L60;
4326  }
4327  L30:
4328  if (*ier != 0 && *ier != 7) {
4329  goto L60;
4330  }
4331  L40:
4332  ll = numrl2;
4333  c1 = c2;
4334  c2 += cycle;
4335  /* L50: */
4336  }
4337 
4338  /* set final result and error estimate */
4339  /* ----------------------------------- */
4340 
4341 L60:
4342  *abserr += correc * 10.;
4343  if (*ier == 0) {
4344  goto L999;
4345  }
4346  if (*result != 0. && psum[numrl2 - 1] != 0.) {
4347  goto L70;
4348  }
4349  if (*abserr > errsum) {
4350  goto L80;
4351  }
4352  if (psum[numrl2 - 1] == 0.) {
4353  goto L999;
4354  }
4355 L70:
4356  if (*abserr / fabs(*result) > (errsum + drl) / (d__1 = psum[numrl2 - 1],
4357  fabs(d__1))) {
4358  goto L80;
4359  }
4360  if (*ier >= 1 && *ier != 7) {
4361  *abserr += drl;
4362  }
4363  goto L999;
4364 L80:
4365  *result = psum[numrl2 - 1];
4366  *abserr = errsum + drl;
4367 L999:
4368  return;
4369 } /* dqawfe_ */
4370 
4371 void dqawf_(const D_fp& f, const double *a, const double *omega, const long *integr,
4372  const double *epsabs, double *result, double *abserr,
4373  long *neval, long *ier, long *limlst, long *lst, const long *
4374  leniw, const long *maxp1, const long *lenw, long *iwork, double *
4375  work)
4376 {
4377  long l1, l2, l3, l4, l5, l6, ll2, lvl, limit;
4378 
4379  /* ***begin prologue dqawf */
4380  /* ***date written 800101 (yymmdd) */
4381  /* ***revision date 830518 (yymmdd) */
4382  /* ***category no. h2a3a1 */
4383  /* ***keywords automatic integrator, special-purpose,fourier */
4384  /* integral, integration between zeros with dqawoe, */
4385  /* convergence acceleration with dqelg */
4386  /* ***author piessens,robert ,appl. math. & progr. div. - k.u.leuven */
4387  /* de doncker,elise,appl. math & progr. div. - k.u.leuven */
4388  /* ***purpose the routine calculates an approximation result to a given */
4389  /* fourier integral i=integral of f(x)*w(x) over (a,infinity) */
4390  /* where w(x) = cos(omega*x) or w(x) = sin(omega*x). */
4391  /* hopefully satisfying following claim for accuracy */
4392  /* fabs(i-result).le.epsabs. */
4393  /* ***description */
4394 
4395  /* computation of fourier integrals */
4396  /* standard fortran subroutine */
4397  /* double precision version */
4398 
4399 
4400  /* parameters */
4401  /* on entry */
4402  /* f - double precision */
4403  /* function subprogram defining the integrand */
4404  /* function f(x). the actual name for f needs to be */
4405  /* declared e x t e r n a l in the driver program. */
4406 
4407  /* a - double precision */
4408  /* lower limit of integration */
4409 
4410  /* omega - double precision */
4411  /* parameter in the integrand weight function */
4412 
4413  /* integr - long */
4414  /* indicates which of the weight functions is used */
4415  /* integr = 1 w(x) = cos(omega*x) */
4416  /* integr = 2 w(x) = sin(omega*x) */
4417  /* if integr.ne.1.and.integr.ne.2, the routine */
4418  /* will end with ier = 6. */
4419 
4420  /* epsabs - double precision */
4421  /* absolute accuracy requested, epsabs.gt.0. */
4422  /* if epsabs.le.0, the routine will end with ier = 6. */
4423 
4424  /* on return */
4425  /* result - double precision */
4426  /* approximation to the integral */
4427 
4428  /* abserr - double precision */
4429  /* estimate of the modulus of the absolute error, */
4430  /* which should equal or exceed fabs(i-result) */
4431 
4432  /* neval - long */
4433  /* number of integrand evaluations */
4434 
4435  /* ier - long */
4436  /* ier = 0 normal and reliable termination of the */
4437  /* routine. it is assumed that the requested */
4438  /* accuracy has been achieved. */
4439  /* ier.gt.0 abnormal termination of the routine. */
4440  /* the estimates for integral and error are */
4441  /* less reliable. it is assumed that the */
4442  /* requested accuracy has not been achieved. */
4443  /* error messages */
4444  /* if omega.ne.0 */
4445  /* ier = 1 maximum number of cycles allowed */
4446  /* has been achieved, i.e. of subintervals */
4447  /* (a+(k-1)c,a+kc) where */
4448  /* c = (2*int(fabs(omega))+1)*pi/fabs(omega), */
4449  /* for k = 1, 2, ..., lst. */
4450  /* one can allow more cycles by increasing */
4451  /* the value of limlst (and taking the */
4452  /* according dimension adjustments into */
4453  /* account). examine the array iwork which */
4454  /* contains the error flags on the cycles, in */
4455  /* order to look for eventual local */
4456  /* integration difficulties. */
4457  /* if the position of a local difficulty */
4458  /* can be determined (e.g. singularity, */
4459  /* discontinuity within the interval) one */
4460  /* will probably gain from splitting up the */
4461  /* interval at this point and calling */
4462  /* appropriate integrators on the subranges. */
4463  /* = 4 the extrapolation table constructed for */
4464  /* convergence accelaration of the series */
4465  /* formed by the integral contributions over */
4466  /* the cycles, does not converge to within */
4467  /* the requested accuracy. */
4468  /* as in the case of ier = 1, it is advised */
4469  /* to examine the array iwork which contains */
4470  /* the error flags on the cycles. */
4471  /* = 6 the input is invalid because */
4472  /* (integr.ne.1 and integr.ne.2) or */
4473  /* epsabs.le.0 or limlst.lt.1 or */
4474  /* leniw.lt.(limlst+2) or maxp1.lt.1 or */
4475  /* lenw.lt.(leniw*2+maxp1*25). */
4476  /* result, abserr, neval, lst are set to */
4477  /* zero. */
4478  /* = 7 bad integrand behaviour occurs within */
4479  /* one or more of the cycles. location and */
4480  /* type of the difficulty involved can be */
4481  /* determined from the first lst elements of */
4482  /* vector iwork. here lst is the number of */
4483  /* cycles actually needed (see below). */
4484  /* iwork(k) = 1 the maximum number of */
4485  /* subdivisions (=(leniw-limlst) */
4486  /* /2) has been achieved on the */
4487  /* k th cycle. */
4488  /* = 2 occurrence of roundoff error */
4489  /* is detected and prevents the */
4490  /* tolerance imposed on the k th */
4491  /* cycle, from being achieved */
4492  /* on this cycle. */
4493  /* = 3 extremely bad integrand */
4494  /* behaviour occurs at some */
4495  /* points of the k th cycle. */
4496  /* = 4 the integration procedure */
4497  /* over the k th cycle does */
4498  /* not converge (to within the */
4499  /* required accuracy) due to */
4500  /* roundoff in the extrapolation */
4501  /* procedure invoked on this */
4502  /* cycle. it is assumed that the */
4503  /* result on this interval is */
4504  /* the best which can be */
4505  /* obtained. */
4506  /* = 5 the integral over the k th */
4507  /* cycle is probably divergent */
4508  /* or slowly convergent. it must */
4509  /* be noted that divergence can */
4510  /* occur with any other value of */
4511  /* iwork(k). */
4512  /* if omega = 0 and integr = 1, */
4513  /* the integral is calculated by means of dqagie, */
4514  /* and ier = iwork(1) (with meaning as described */
4515  /* for iwork(k),k = 1). */
4516 
4517  /* dimensioning parameters */
4518  /* limlst - long */
4519  /* limlst gives an upper bound on the number of */
4520  /* cycles, limlst.ge.3. */
4521  /* if limlst.lt.3, the routine will end with ier = 6. */
4522 
4523  /* lst - long */
4524  /* on return, lst indicates the number of cycles */
4525  /* actually needed for the integration. */
4526  /* if omega = 0, then lst is set to 1. */
4527 
4528  /* leniw - long */
4529  /* dimensioning parameter for iwork. on entry, */
4530  /* (leniw-limlst)/2 equals the maximum number of */
4531  /* subintervals allowed in the partition of each */
4532  /* cycle, leniw.ge.(limlst+2). */
4533  /* if leniw.lt.(limlst+2), the routine will end with */
4534  /* ier = 6. */
4535 
4536  /* maxp1 - long */
4537  /* maxp1 gives an upper bound on the number of */
4538  /* chebyshev moments which can be stored, i.e. for */
4539  /* the intervals of lengths fabs(b-a)*2**(-l), */
4540  /* l = 0,1, ..., maxp1-2, maxp1.ge.1. */
4541  /* if maxp1.lt.1, the routine will end with ier = 6. */
4542  /* lenw - long */
4543  /* dimensioning parameter for work */
4544  /* lenw must be at least leniw*2+maxp1*25. */
4545  /* if lenw.lt.(leniw*2+maxp1*25), the routine will */
4546  /* end with ier = 6. */
4547 
4548  /* work arrays */
4549  /* iwork - long */
4550  /* vector of dimension at least leniw */
4551  /* on return, iwork(k) for k = 1, 2, ..., lst */
4552  /* contain the error flags on the cycles. */
4553 
4554  /* work - double precision */
4555  /* vector of dimension at least */
4556  /* on return, */
4557  /* work(1), ..., work(lst) contain the integral */
4558  /* approximations over the cycles, */
4559  /* work(limlst+1), ..., work(limlst+lst) contain */
4560  /* the error extimates over the cycles. */
4561  /* further elements of work have no specific */
4562  /* meaning for the user. */
4563 
4564  /* ***references (none) */
4565  /* ***routines called dqawfe,xerror */
4566  /* ***end prologue dqawf */
4567 
4568 
4569 
4570 
4571  /* check validity of limlst, leniw, maxp1 and lenw. */
4572 
4573  /* ***first executable statement dqawf */
4574  /* Parameter adjustments */
4575  --iwork;
4576  --work;
4577 
4578  /* Function Body */
4579  *ier = 6;
4580  *neval = 0;
4581  *result = 0.;
4582  *abserr = 0.;
4583  if (*limlst < 3 || *leniw < *limlst + 2 || *maxp1 < 1 || *lenw < (*leniw
4584  << 1) + *maxp1 * 25) {
4585  goto L10;
4586  }
4587 
4588  /* prepare call for dqawfe */
4589 
4590  limit = (*leniw - *limlst) / 2;
4591  l1 = *limlst + 1;
4592  l2 = *limlst + l1;
4593  l3 = limit + l2;
4594  l4 = limit + l3;
4595  l5 = limit + l4;
4596  l6 = limit + l5;
4597  ll2 = limit + l1;
4598  dqawfe_(f, a, omega, integr, epsabs, limlst, &limit, maxp1, result,
4599  abserr, neval, ier, &work[1], &work[l1], &iwork[1], lst, &work[l2]
4600  , &work[l3], &work[l4], &work[l5], &iwork[l1], &iwork[ll2], &work[
4601  l6]);
4602 
4603  /* call error handler if necessary */
4604 
4605  lvl = 0;
4606 L10:
4607  if (*ier == 6) {
4608  lvl = 1;
4609  }
4610  if (*ier != 0) {
4611  xerror_("abnormal return from dqawf", &c__26, ier, &lvl, 26);
4612  }
4613  return;
4614 } /* dqawf_ */
4615 
4616 void dqawoe_(const D_fp& f, const double *a, const double *b, const double
4617  *omega, const long *integr, const double *epsabs,
4618  const double *epsrel,
4619  const long *limit, const long *icall, const long *maxp1,
4620  double *result,
4621  double *abserr, long *neval, long *ier, long *last,
4622  double *alist__, double *blist, double *rlist, double
4623  *elist, long *iord, long *nnlog, long *momcom, double *
4624  chebmo)
4625 {
4626  /* System generated locals */
4627  int chebmo_dim1, chebmo_offset, i__1, i__2;
4628  double d__1, d__2;
4629 
4630  /* Local variables */
4631  long k;
4632  double a1, a2, b1, b2;
4633  long id, nev;
4634  double area, dres;
4635  long ksgn, nres;
4636  double area1, area2, area12;
4637  double small, erro12, width, defab1, defab2;
4638  long ierro, ktmin;
4639  double oflow;
4640  long nrmax, nrmom;
4641  double uflow;
4642  bool noext;
4643  long iroff1, iroff2, iroff3;
4644  double res3la[3], error1, error2, rlist2[52];
4645  long numrl2;
4646  double defabs, domega, epmach, erlarg=0., abseps, correc=0., errbnd,
4647  resabs;
4648  long jupbnd;
4649  bool extall;
4650  double erlast, errmax;
4651  long maxerr;
4652  double reseps;
4653  bool extrap;
4654  double ertest=0., errsum;
4655 
4656  /* ***begin prologue dqawoe */
4657  /* ***date written 800101 (yymmdd) */
4658  /* ***revision date 830518 (yymmdd) */
4659  /* ***category no. h2a2a1 */
4660  /* ***keywords automatic integrator, special-purpose, */
4661  /* integrand with oscillatory cos or sin factor, */
4662  /* clenshaw-curtis method, (end point) singularities, */
4663  /* extrapolation, globally adaptive */
4664  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
4665  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
4666  /* ***purpose the routine calculates an approximation result to a given */
4667  /* definite integral */
4668  /* i = integral of f(x)*w(x) over (a,b) */
4669  /* where w(x) = cos(omega*x) or w(x)=sin(omega*x), */
4670  /* hopefully satisfying following claim for accuracy */
4671  /* fabs(i-result).le.max(epsabs,epsrel*fabs(i)). */
4672  /* ***description */
4673 
4674  /* computation of oscillatory integrals */
4675  /* standard fortran subroutine */
4676  /* double precision version */
4677 
4678  /* parameters */
4679  /* on entry */
4680  /* f - double precision */
4681  /* function subprogram defining the integrand */
4682  /* function f(x). the actual name for f needs to be */
4683  /* declared e x t e r n a l in the driver program. */
4684 
4685  /* a - double precision */
4686  /* lower limit of integration */
4687 
4688  /* b - double precision */
4689  /* upper limit of integration */
4690 
4691  /* omega - double precision */
4692  /* parameter in the integrand weight function */
4693 
4694  /* integr - long */
4695  /* indicates which of the weight functions is to be */
4696  /* used */
4697  /* integr = 1 w(x) = cos(omega*x) */
4698  /* integr = 2 w(x) = sin(omega*x) */
4699  /* if integr.ne.1 and integr.ne.2, the routine */
4700  /* will end with ier = 6. */
4701 
4702  /* epsabs - double precision */
4703  /* absolute accuracy requested */
4704  /* epsrel - double precision */
4705  /* relative accuracy requested */
4706  /* if epsabs.le.0 */
4707  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
4708  /* the routine will end with ier = 6. */
4709 
4710  /* limit - long */
4711  /* gives an upper bound on the number of subdivisions */
4712  /* in the partition of (a,b), limit.ge.1. */
4713 
4714  /* icall - long */
4715  /* if dqawoe is to be used only once, icall must */
4716  /* be set to 1. assume that during this call, the */
4717  /* chebyshev moments (for clenshaw-curtis integration */
4718  /* of degree 24) have been computed for intervals of */
4719  /* lenghts (fabs(b-a))*2**(-l), l=0,1,2,...momcom-1. */
4720  /* if icall.gt.1 this means that dqawoe has been */
4721  /* called twice or more on intervals of the same */
4722  /* length fabs(b-a). the chebyshev moments already */
4723  /* computed are then re-used in subsequent calls. */
4724  /* if icall.lt.1, the routine will end with ier = 6. */
4725 
4726  /* maxp1 - long */
4727  /* gives an upper bound on the number of chebyshev */
4728  /* moments which can be stored, i.e. for the */
4729  /* intervals of lenghts fabs(b-a)*2**(-l), */
4730  /* l=0,1, ..., maxp1-2, maxp1.ge.1. */
4731  /* if maxp1.lt.1, the routine will end with ier = 6. */
4732 
4733  /* on return */
4734  /* result - double precision */
4735  /* approximation to the integral */
4736 
4737  /* abserr - double precision */
4738  /* estimate of the modulus of the absolute error, */
4739  /* which should equal or exceed fabs(i-result) */
4740 
4741  /* neval - long */
4742  /* number of integrand evaluations */
4743 
4744  /* ier - long */
4745  /* ier = 0 normal and reliable termination of the */
4746  /* routine. it is assumed that the */
4747  /* requested accuracy has been achieved. */
4748  /* - ier.gt.0 abnormal termination of the routine. */
4749  /* the estimates for integral and error are */
4750  /* less reliable. it is assumed that the */
4751  /* requested accuracy has not been achieved. */
4752  /* error messages */
4753  /* ier = 1 maximum number of subdivisions allowed */
4754  /* has been achieved. one can allow more */
4755  /* subdivisions by increasing the value of */
4756  /* limit (and taking according dimension */
4757  /* adjustments into account). however, if */
4758  /* this yields no improvement it is advised */
4759  /* to analyze the integrand, in order to */
4760  /* determine the integration difficulties. */
4761  /* if the position of a local difficulty can */
4762  /* be determined (e.g. singularity, */
4763  /* discontinuity within the interval) one */
4764  /* will probably gain from splitting up the */
4765  /* interval at this point and calling the */
4766  /* integrator on the subranges. if possible, */
4767  /* an appropriate special-purpose integrator */
4768  /* should be used which is designed for */
4769  /* handling the type of difficulty involved. */
4770  /* = 2 the occurrence of roundoff error is */
4771  /* detected, which prevents the requested */
4772  /* tolerance from being achieved. */
4773  /* the error may be under-estimated. */
4774  /* = 3 extremely bad integrand behaviour occurs */
4775  /* at some points of the integration */
4776  /* interval. */
4777  /* = 4 the algorithm does not converge. */
4778  /* roundoff error is detected in the */
4779  /* extrapolation table. */
4780  /* it is presumed that the requested */
4781  /* tolerance cannot be achieved due to */
4782  /* roundoff in the extrapolation table, */
4783  /* and that the returned result is the */
4784  /* best which can be obtained. */
4785  /* = 5 the integral is probably divergent, or */
4786  /* slowly convergent. it must be noted that */
4787  /* divergence can occur with any other value */
4788  /* of ier.gt.0. */
4789  /* = 6 the input is invalid, because */
4790  /* (epsabs.le.0 and */
4791  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
4792  /* or (integr.ne.1 and integr.ne.2) or */
4793  /* icall.lt.1 or maxp1.lt.1. */
4794  /* result, abserr, neval, last, rlist(1), */
4795  /* elist(1), iord(1) and nnlog(1) are set */
4796  /* to zero. alist(1) and blist(1) are set */
4797  /* to a and b respectively. */
4798 
4799  /* last - long */
4800  /* on return, last equals the number of */
4801  /* subintervals produces in the subdivision */
4802  /* process, which determines the number of */
4803  /* significant elements actually in the */
4804  /* work arrays. */
4805  /* alist - double precision */
4806  /* vector of dimension at least limit, the first */
4807  /* last elements of which are the left */
4808  /* end points of the subintervals in the partition */
4809  /* of the given integration range (a,b) */
4810 
4811  /* blist - double precision */
4812  /* vector of dimension at least limit, the first */
4813  /* last elements of which are the right */
4814  /* end points of the subintervals in the partition */
4815  /* of the given integration range (a,b) */
4816 
4817  /* rlist - double precision */
4818  /* vector of dimension at least limit, the first */
4819  /* last elements of which are the integral */
4820  /* approximations on the subintervals */
4821 
4822  /* elist - double precision */
4823  /* vector of dimension at least limit, the first */
4824  /* last elements of which are the moduli of the */
4825  /* absolute error estimates on the subintervals */
4826 
4827  /* iord - long */
4828  /* vector of dimension at least limit, the first k */
4829  /* elements of which are pointers to the error */
4830  /* estimates over the subintervals, */
4831  /* such that elist(iord(1)), ..., */
4832  /* elist(iord(k)) form a decreasing sequence, with */
4833  /* k = last if last.le.(limit/2+2), and */
4834  /* k = limit+1-last otherwise. */
4835 
4836  /* nnlog - long */
4837  /* vector of dimension at least limit, containing the */
4838  /* subdivision levels of the subintervals, i.e. */
4839  /* iwork(i) = l means that the subinterval */
4840  /* numbered i is of length fabs(b-a)*2**(1-l) */
4841 
4842  /* on entry and return */
4843  /* momcom - long */
4844  /* indicating that the chebyshev moments */
4845  /* have been computed for intervals of lengths */
4846  /* (fabs(b-a))*2**(-l), l=0,1,2, ..., momcom-1, */
4847  /* momcom.lt.maxp1 */
4848 
4849  /* chebmo - double precision */
4850  /* array of dimension (maxp1,25) containing the */
4851  /* chebyshev moments */
4852 
4853  /* ***references (none) */
4854  /* ***routines called d1mach,dqc25f,dqelg,dqpsrt */
4855  /* ***end prologue dqawoe */
4856 
4857 
4858 
4859 
4860  /* the dimension of rlist2 is determined by the value of */
4861  /* limexp in subroutine dqelg (rlist2 should be of */
4862  /* dimension (limexp+2) at least). */
4863 
4864  /* list of major variables */
4865  /* ----------------------- */
4866 
4867  /* alist - list of left end points of all subintervals */
4868  /* considered up to now */
4869  /* blist - list of right end points of all subintervals */
4870  /* considered up to now */
4871  /* rlist(i) - approximation to the integral over */
4872  /* (alist(i),blist(i)) */
4873  /* rlist2 - array of dimension at least limexp+2 */
4874  /* containing the part of the epsilon table */
4875  /* which is still needed for further computations */
4876  /* elist(i) - error estimate applying to rlist(i) */
4877  /* maxerr - pointer to the interval with largest */
4878  /* error estimate */
4879  /* errmax - elist(maxerr) */
4880  /* erlast - error on the interval currently subdivided */
4881  /* area - sum of the integrals over the subintervals */
4882  /* errsum - sum of the errors over the subintervals */
4883  /* errbnd - requested accuracy max(epsabs,epsrel* */
4884  /* fabs(result)) */
4885  /* *****1 - variable for the left subinterval */
4886  /* *****2 - variable for the right subinterval */
4887  /* last - index for subdivision */
4888  /* nres - number of calls to the extrapolation routine */
4889  /* numrl2 - number of elements in rlist2. if an appropriate */
4890  /* approximation to the compounded integral has */
4891  /* been obtained it is put in rlist2(numrl2) after */
4892  /* numrl2 has been increased by one */
4893  /* small - length of the smallest interval considered */
4894  /* up to now, multiplied by 1.5 */
4895  /* erlarg - sum of the errors over the intervals larger */
4896  /* than the smallest interval considered up to now */
4897  /* extrap - bool variable denoting that the routine is */
4898  /* attempting to perform extrapolation, i.e. before */
4899  /* subdividing the smallest interval we try to */
4900  /* decrease the value of erlarg */
4901  /* noext - bool variable denoting that extrapolation */
4902  /* is no longer allowed (true value) */
4903 
4904  /* machine dependent constants */
4905  /* --------------------------- */
4906 
4907  /* epmach is the largest relative spacing. */
4908  /* uflow is the smallest positive magnitude. */
4909  /* oflow is the largest positive magnitude. */
4910 
4911  /* ***first executable statement dqawoe */
4912  /* Parameter adjustments */
4913  --nnlog;
4914  --iord;
4915  --elist;
4916  --rlist;
4917  --blist;
4918  --alist__;
4919  chebmo_dim1 = *maxp1;
4920  chebmo_offset = 1 + chebmo_dim1;
4921  chebmo -= chebmo_offset;
4922 
4923  /* Function Body */
4924  epmach = d1mach(c__4);
4925 
4926  /* test on validity of parameters */
4927  /* ------------------------------ */
4928 
4929  *ier = 0;
4930  *neval = 0;
4931  *last = 0;
4932  *result = 0.;
4933  *abserr = 0.;
4934  alist__[1] = *a;
4935  blist[1] = *b;
4936  rlist[1] = 0.;
4937  elist[1] = 0.;
4938  iord[1] = 0;
4939  nnlog[1] = 0;
4940  /* Computing MAX */
4941  d__1 = epmach * 50.;
4942  if ( ( *integr != 1 && *integr != 2 ) ||
4943  ( *epsabs <= 0. && *epsrel < max(d__1, 5e-29) ) ||
4944  *icall < 1 || *maxp1 < 1) {
4945  *ier = 6;
4946  }
4947  if (*ier == 6) {
4948  goto L999;
4949  }
4950 
4951  /* first approximation to the integral */
4952  /* ----------------------------------- */
4953 
4954  domega = fabs(*omega);
4955  nrmom = 0;
4956  if (*icall > 1) {
4957  goto L5;
4958  }
4959  *momcom = 0;
4960 L5:
4961  dqc25f_(f, a, b, &domega, integr, &nrmom, maxp1, &c__0, result,
4962  abserr, neval, &defabs, &resabs, momcom, &chebmo[chebmo_offset]);
4963 
4964  /* test on accuracy. */
4965 
4966  dres = fabs(*result);
4967  /* Computing MAX */
4968  d__1 = *epsabs, d__2 = *epsrel * dres;
4969  errbnd = max(d__1,d__2);
4970  rlist[1] = *result;
4971  elist[1] = *abserr;
4972  iord[1] = 1;
4973  if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) {
4974  *ier = 2;
4975  }
4976  if (*limit == 1) {
4977  *ier = 1;
4978  }
4979  if (*ier != 0 || *abserr <= errbnd) {
4980  goto L200;
4981  }
4982 
4983  /* initializations */
4984  /* --------------- */
4985 
4986  uflow = d1mach(c__1);
4987  oflow = d1mach(c__2);
4988  errmax = *abserr;
4989  maxerr = 1;
4990  area = *result;
4991  errsum = *abserr;
4992  *abserr = oflow;
4993  nrmax = 1;
4994  extrap = false;
4995  noext = false;
4996  ierro = 0;
4997  iroff1 = 0;
4998  iroff2 = 0;
4999  iroff3 = 0;
5000  ktmin = 0;
5001  small = (d__1 = *b - *a, fabs(d__1)) * .75;
5002  nres = 0;
5003  numrl2 = 0;
5004  extall = false;
5005  if ((d__1 = *b - *a, fabs(d__1)) * .5 * domega > 2.) {
5006  goto L10;
5007  }
5008  numrl2 = 1;
5009  extall = true;
5010  rlist2[0] = *result;
5011 L10:
5012  if ((d__1 = *b - *a, fabs(d__1)) * .25 * domega <= 2.) {
5013  extall = true;
5014  }
5015  ksgn = -1;
5016  if (dres >= (1. - epmach * 50.) * defabs) {
5017  ksgn = 1;
5018  }
5019 
5020  /* main do-loop */
5021  /* ------------ */
5022 
5023  i__1 = *limit;
5024  for (*last = 2; *last <= i__1; ++(*last)) {
5025 
5026  /* bisect the subinterval with the nrmax-th largest */
5027  /* error estimate. */
5028 
5029  nrmom = nnlog[maxerr] + 1;
5030  a1 = alist__[maxerr];
5031  b1 = (alist__[maxerr] + blist[maxerr]) * .5;
5032  a2 = b1;
5033  b2 = blist[maxerr];
5034  erlast = errmax;
5035  dqc25f_(f, &a1, &b1, &domega, integr, &nrmom, maxp1, &c__0, &
5036  area1, &error1, &nev, &resabs, &defab1, momcom, &chebmo[
5037  chebmo_offset]);
5038  *neval += nev;
5039  dqc25f_(f, &a2, &b2, &domega, integr, &nrmom, maxp1, &c__1, &
5040  area2, &error2, &nev, &resabs, &defab2, momcom, &chebmo[
5041  chebmo_offset]);
5042  *neval += nev;
5043 
5044  /* improve previous approximations to integral */
5045  /* and error and test for accuracy. */
5046 
5047  area12 = area1 + area2;
5048  erro12 = error1 + error2;
5049  errsum = errsum + erro12 - errmax;
5050  area = area + area12 - rlist[maxerr];
5051  if (defab1 == error1 || defab2 == error2) {
5052  goto L25;
5053  }
5054  if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
5055  erro12 < errmax * .99) {
5056  goto L20;
5057  }
5058  if (extrap) {
5059  ++iroff2;
5060  }
5061  if (! extrap) {
5062  ++iroff1;
5063  }
5064  L20:
5065  if (*last > 10 && erro12 > errmax) {
5066  ++iroff3;
5067  }
5068  L25:
5069  rlist[maxerr] = area1;
5070  rlist[*last] = area2;
5071  nnlog[maxerr] = nrmom;
5072  nnlog[*last] = nrmom;
5073  /* Computing MAX */
5074  d__1 = *epsabs, d__2 = *epsrel * fabs(area);
5075  errbnd = max(d__1,d__2);
5076 
5077  /* test for roundoff error and eventually set error flag. */
5078 
5079  if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
5080  *ier = 2;
5081  }
5082  if (iroff2 >= 5) {
5083  ierro = 3;
5084  }
5085 
5086  /* set error flag in the case that the number of */
5087  /* subintervals equals limit. */
5088 
5089  if (*last == *limit) {
5090  *ier = 1;
5091  }
5092 
5093  /* set error flag in the case of bad integrand behaviour */
5094  /* at a point of the integration range. */
5095 
5096  /* Computing MAX */
5097  d__1 = fabs(a1), d__2 = fabs(b2);
5098  if (max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
5099  {
5100  *ier = 4;
5101  }
5102 
5103  /* append the newly-created intervals to the list. */
5104 
5105  if (error2 > error1) {
5106  goto L30;
5107  }
5108  alist__[*last] = a2;
5109  blist[maxerr] = b1;
5110  blist[*last] = b2;
5111  elist[maxerr] = error1;
5112  elist[*last] = error2;
5113  goto L40;
5114  L30:
5115  alist__[maxerr] = a2;
5116  alist__[*last] = a1;
5117  blist[*last] = b1;
5118  rlist[maxerr] = area2;
5119  rlist[*last] = area1;
5120  elist[maxerr] = error2;
5121  elist[*last] = error1;
5122 
5123  /* call subroutine dqpsrt to maintain the descending ordering */
5124  /* in the list of error estimates and select the subinterval */
5125  /* with nrmax-th largest error estimate (to bisected next). */
5126 
5127  L40:
5128  dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
5129  /* ***jump out of do-loop */
5130  if (errsum <= errbnd) {
5131  goto L170;
5132  }
5133  if (*ier != 0) {
5134  goto L150;
5135  }
5136  if (*last == 2 && extall) {
5137  goto L120;
5138  }
5139  if (noext) {
5140  goto L140;
5141  }
5142  if (! extall) {
5143  goto L50;
5144  }
5145  erlarg -= erlast;
5146  if ((d__1 = b1 - a1, fabs(d__1)) > small) {
5147  erlarg += erro12;
5148  }
5149  if (extrap) {
5150  goto L70;
5151  }
5152 
5153  /* test whether the interval to be bisected next is the */
5154  /* smallest interval. */
5155 
5156  L50:
5157  width = (d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1));
5158  if (width > small) {
5159  goto L140;
5160  }
5161  if (extall) {
5162  goto L60;
5163  }
5164 
5165  /* test whether we can start with the extrapolation procedure */
5166  /* (we do this if we integrate over the next interval with */
5167  /* use of a gauss-kronrod rule - see subroutine dqc25f). */
5168 
5169  small *= .5;
5170  if (width * .25 * domega > 2.) {
5171  goto L140;
5172  }
5173  extall = true;
5174  goto L130;
5175  L60:
5176  extrap = true;
5177  nrmax = 2;
5178  L70:
5179  if (ierro == 3 || erlarg <= ertest) {
5180  goto L90;
5181  }
5182 
5183  /* the smallest interval has the largest error. */
5184  /* before bisecting decrease the sum of the errors over */
5185  /* the larger intervals (erlarg) and perform extrapolation. */
5186 
5187  jupbnd = *last;
5188  if (*last > *limit / 2 + 2) {
5189  jupbnd = *limit + 3 - *last;
5190  }
5191  id = nrmax;
5192  i__2 = jupbnd;
5193  for (k = id; k <= i__2; ++k) {
5194  maxerr = iord[nrmax];
5195  errmax = elist[maxerr];
5196  if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
5197  goto L140;
5198  }
5199  ++nrmax;
5200  /* L80: */
5201  }
5202 
5203  /* perform extrapolation. */
5204 
5205  L90:
5206  ++numrl2;
5207  rlist2[numrl2 - 1] = area;
5208  if (numrl2 < 3) {
5209  goto L110;
5210  }
5211  dqelg_(&numrl2, rlist2, &reseps, &abseps, res3la, &nres);
5212  ++ktmin;
5213  if (ktmin > 5 && *abserr < errsum * .001) {
5214  *ier = 5;
5215  }
5216  if (abseps >= *abserr) {
5217  goto L100;
5218  }
5219  ktmin = 0;
5220  *abserr = abseps;
5221  *result = reseps;
5222  correc = erlarg;
5223  /* Computing MAX */
5224  d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
5225  ertest = max(d__1,d__2);
5226  /* ***jump out of do-loop */
5227  if (*abserr <= ertest) {
5228  goto L150;
5229  }
5230 
5231  /* prepare bisection of the smallest interval. */
5232 
5233  L100:
5234  if (numrl2 == 1) {
5235  noext = true;
5236  }
5237  if (*ier == 5) {
5238  goto L150;
5239  }
5240  L110:
5241  maxerr = iord[1];
5242  errmax = elist[maxerr];
5243  nrmax = 1;
5244  extrap = false;
5245  small *= .5;
5246  erlarg = errsum;
5247  goto L140;
5248  L120:
5249  small *= .5;
5250  ++numrl2;
5251  rlist2[numrl2 - 1] = area;
5252  L130:
5253  ertest = errbnd;
5254  erlarg = errsum;
5255  L140:
5256  ;
5257  }
5258 
5259  /* set the final result. */
5260  /* --------------------- */
5261 
5262 L150:
5263  if (*abserr == oflow || nres == 0) {
5264  goto L170;
5265  }
5266  if (*ier + ierro == 0) {
5267  goto L165;
5268  }
5269  if (ierro == 3) {
5270  *abserr += correc;
5271  }
5272  if (*ier == 0) {
5273  *ier = 3;
5274  }
5275  if (*result != 0. && area != 0.) {
5276  goto L160;
5277  }
5278  if (*abserr > errsum) {
5279  goto L170;
5280  }
5281  if (area == 0.) {
5282  goto L190;
5283  }
5284  goto L165;
5285 L160:
5286  if (*abserr / fabs(*result) > errsum / fabs(area)) {
5287  goto L170;
5288  }
5289 
5290  /* test on divergence. */
5291 
5292 L165:
5293  /* Computing MAX */
5294  d__1 = fabs(*result), d__2 = fabs(area);
5295  if (ksgn == -1 && max(d__1,d__2) <= defabs * .01) {
5296  goto L190;
5297  }
5298  if (.01 > *result / area || *result / area > 100. || errsum >= fabs(area))
5299  {
5300  *ier = 6;
5301  }
5302  goto L190;
5303 
5304  /* compute global integral sum. */
5305 
5306 L170:
5307  *result = 0.;
5308  i__1 = *last;
5309  for (k = 1; k <= i__1; ++k) {
5310  *result += rlist[k];
5311  /* L180: */
5312  }
5313  *abserr = errsum;
5314 L190:
5315  if (*ier > 2) {
5316  --(*ier);
5317  }
5318 L200:
5319  if (*integr == 2 && *omega < 0.) {
5320  *result = -(*result);
5321  }
5322 L999:
5323  return;
5324 } /* dqawoe_ */
5325 
5326 void dqawo_(const D_fp& f, const double *a, const double *b, const double *omega,
5327  const long *integr, const double *epsabs, const double *epsrel,
5328  double *result, double *abserr, long *neval, long *ier,
5329  const long *leniw, long *maxp1, const long *lenw, long *last,
5330  long *iwork, double *work)
5331 {
5332  long l1, l2, l3, l4, lvl, limit;
5333  long momcom;
5334 
5335  /* ***begin prologue dqawo */
5336  /* ***date written 800101 (yymmdd) */
5337  /* ***revision date 830518 (yymmdd) */
5338  /* ***category no. h2a2a1 */
5339  /* ***keywords automatic integrator, special-purpose, */
5340  /* integrand with oscillatory cos or sin factor, */
5341  /* clenshaw-curtis method, (end point) singularities, */
5342  /* extrapolation, globally adaptive */
5343  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
5344  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
5345  /* ***purpose the routine calculates an approximation result to a given */
5346  /* definite integral i=integral of f(x)*w(x) over (a,b) */
5347  /* where w(x) = cos(omega*x) */
5348  /* or w(x) = sin(omega*x), */
5349  /* hopefully satisfying following claim for accuracy */
5350  /* fabs(i-result).le.max(epsabs,epsrel*fabs(i)). */
5351  /* ***description */
5352 
5353  /* computation of oscillatory integrals */
5354  /* standard fortran subroutine */
5355  /* double precision version */
5356 
5357  /* parameters */
5358  /* on entry */
5359  /* f - double precision */
5360  /* function subprogram defining the function */
5361  /* f(x). the actual name for f needs to be */
5362  /* declared e x t e r n a l in the driver program. */
5363 
5364  /* a - double precision */
5365  /* lower limit of integration */
5366 
5367  /* b - double precision */
5368  /* upper limit of integration */
5369 
5370  /* omega - double precision */
5371  /* parameter in the integrand weight function */
5372 
5373  /* integr - long */
5374  /* indicates which of the weight functions is used */
5375  /* integr = 1 w(x) = cos(omega*x) */
5376  /* integr = 2 w(x) = sin(omega*x) */
5377  /* if integr.ne.1.and.integr.ne.2, the routine will */
5378  /* end with ier = 6. */
5379 
5380  /* epsabs - double precision */
5381  /* absolute accuracy requested */
5382  /* epsrel - double precision */
5383  /* relative accuracy requested */
5384  /* if epsabs.le.0 and */
5385  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
5386  /* the routine will end with ier = 6. */
5387 
5388  /* on return */
5389  /* result - double precision */
5390  /* approximation to the integral */
5391 
5392  /* abserr - double precision */
5393  /* estimate of the modulus of the absolute error, */
5394  /* which should equal or exceed fabs(i-result) */
5395 
5396  /* neval - long */
5397  /* number of integrand evaluations */
5398 
5399  /* ier - long */
5400  /* ier = 0 normal and reliable termination of the */
5401  /* routine. it is assumed that the requested */
5402  /* accuracy has been achieved. */
5403  /* - ier.gt.0 abnormal termination of the routine. */
5404  /* the estimates for integral and error are */
5405  /* less reliable. it is assumed that the */
5406  /* requested accuracy has not been achieved. */
5407  /* error messages */
5408  /* ier = 1 maximum number of subdivisions allowed */
5409  /* (= leniw/2) has been achieved. one can */
5410  /* allow more subdivisions by increasing the */
5411  /* value of leniw (and taking the according */
5412  /* dimension adjustments into account). */
5413  /* however, if this yields no improvement it */
5414  /* is advised to analyze the integrand in */
5415  /* order to determine the integration */
5416  /* difficulties. if the position of a local */
5417  /* difficulty can be determined (e.g. */
5418  /* singularity, discontinuity within the */
5419  /* interval) one will probably gain from */
5420  /* splitting up the interval at this polong */
5421  /* and calling the integrator on the */
5422  /* subranges. if possible, an appropriate */
5423  /* special-purpose integrator should be used */
5424  /* which is designed for handling the type of */
5425  /* difficulty involved. */
5426  /* = 2 the occurrence of roundoff error is */
5427  /* detected, which prevents the requested */
5428  /* tolerance from being achieved. */
5429  /* the error may be under-estimated. */
5430  /* = 3 extremely bad integrand behaviour occurs */
5431  /* at some interior points of the */
5432  /* integration interval. */
5433  /* = 4 the algorithm does not converge. */
5434  /* roundoff error is detected in the */
5435  /* extrapolation table. it is presumed that */
5436  /* the requested tolerance cannot be achieved */
5437  /* due to roundoff in the extrapolation */
5438  /* table, and that the returned result is */
5439  /* the best which can be obtained. */
5440  /* = 5 the integral is probably divergent, or */
5441  /* slowly convergent. it must be noted that */
5442  /* divergence can occur with any other value */
5443  /* of ier. */
5444  /* = 6 the input is invalid, because */
5445  /* (epsabs.le.0 and */
5446  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
5447  /* or (integr.ne.1 and integr.ne.2), */
5448  /* or leniw.lt.2 or maxp1.lt.1 or */
5449  /* lenw.lt.leniw*2+maxp1*25. */
5450  /* result, abserr, neval, last are set to */
5451  /* zero. except when leniw, maxp1 or lenw are */
5452  /* invalid, work(limit*2+1), work(limit*3+1), */
5453  /* iwork(1), iwork(limit+1) are set to zero, */
5454  /* work(1) is set to a and work(limit+1) to */
5455  /* b. */
5456 
5457  /* dimensioning parameters */
5458  /* leniw - long */
5459  /* dimensioning parameter for iwork. */
5460  /* leniw/2 equals the maximum number of subintervals */
5461  /* allowed in the partition of the given integration */
5462  /* interval (a,b), leniw.ge.2. */
5463  /* if leniw.lt.2, the routine will end with ier = 6. */
5464 
5465  /* maxp1 - long */
5466  /* gives an upper bound on the number of chebyshev */
5467  /* moments which can be stored, i.e. for the */
5468  /* intervals of lengths fabs(b-a)*2**(-l), */
5469  /* l=0,1, ..., maxp1-2, maxp1.ge.1 */
5470  /* if maxp1.lt.1, the routine will end with ier = 6. */
5471 
5472  /* lenw - long */
5473  /* dimensioning parameter for work */
5474  /* lenw must be at least leniw*2+maxp1*25. */
5475  /* if lenw.lt.(leniw*2+maxp1*25), the routine will */
5476  /* end with ier = 6. */
5477 
5478  /* last - long */
5479  /* on return, last equals the number of subintervals */
5480  /* produced in the subdivision process, which */
5481  /* determines the number of significant elements */
5482  /* actually in the work arrays. */
5483 
5484  /* work arrays */
5485  /* iwork - long */
5486  /* vector of dimension at least leniw */
5487  /* on return, the first k elements of which contain */
5488  /* pointers to the error estimates over the */
5489  /* subintervals, such that work(limit*3+iwork(1)), .. */
5490  /* work(limit*3+iwork(k)) form a decreasing */
5491  /* sequence, with limit = lenw/2 , and k = last */
5492  /* if last.le.(limit/2+2), and k = limit+1-last */
5493  /* otherwise. */
5494  /* furthermore, iwork(limit+1), ..., iwork(limit+ */
5495  /* last) indicate the subdivision levels of the */
5496  /* subintervals, such that iwork(limit+i) = l means */
5497  /* that the subinterval numbered i is of length */
5498  /* fabs(b-a)*2**(1-l). */
5499 
5500  /* work - double precision */
5501  /* vector of dimension at least lenw */
5502  /* on return */
5503  /* work(1), ..., work(last) contain the left */
5504  /* end points of the subintervals in the */
5505  /* partition of (a,b), */
5506  /* work(limit+1), ..., work(limit+last) contain */
5507  /* the right end points, */
5508  /* work(limit*2+1), ..., work(limit*2+last) contain */
5509  /* the integral approximations over the */
5510  /* subintervals, */
5511  /* work(limit*3+1), ..., work(limit*3+last) */
5512  /* contain the error estimates. */
5513  /* work(limit*4+1), ..., work(limit*4+maxp1*25) */
5514  /* provide space for storing the chebyshev moments. */
5515  /* note that limit = lenw/2. */
5516 
5517  /* ***references (none) */
5518  /* ***routines called dqawoe,xerror */
5519  /* ***end prologue dqawo */
5520 
5521 
5522 
5523 
5524  /* check validity of leniw, maxp1 and lenw. */
5525 
5526  /* ***first executable statement dqawo */
5527  /* Parameter adjustments */
5528  --iwork;
5529  --work;
5530 
5531  /* Function Body */
5532  *ier = 6;
5533  *neval = 0;
5534  *last = 0;
5535  *result = 0.;
5536  *abserr = 0.;
5537  if (*leniw < 2 || *maxp1 < 1 || *lenw < (*leniw << 1) + *maxp1 * 25) {
5538  goto L10;
5539  }
5540 
5541  /* prepare call for dqawoe */
5542 
5543  limit = *leniw / 2;
5544  l1 = limit + 1;
5545  l2 = limit + l1;
5546  l3 = limit + l2;
5547  l4 = limit + l3;
5548  dqawoe_(f, a, b, omega, integr, epsabs, epsrel, &limit, &c__1,
5549  maxp1, result, abserr, neval, ier, last, &work[1], &work[l1], &
5550  work[l2], &work[l3], &iwork[1], &iwork[l1], &momcom, &work[l4]);
5551 
5552  /* call error handler if necessary */
5553 
5554  lvl = 0;
5555 L10:
5556  if (*ier == 6) {
5557  lvl = 0;
5558  }
5559  if (*ier != 0) {
5560  xerror_("abnormal return from dqawo", &c__26, ier, &lvl, 26);
5561  }
5562  return;
5563 } /* dqawo_ */
5564 
5565 void dqawse_(const D_fp& f, const double *a, const double *b, const double *alfa,
5566  const double *beta, const long *integr, const double *epsabs,
5567  const double *epsrel, const long *limit, double *result,
5568  double *abserr, long *neval, long *ier, double *alist__,
5569  double *blist, double *rlist, double *elist, long *iord,
5570  long *last)
5571 {
5572  /* System generated locals */
5573  int i__1;
5574  double d__1, d__2;
5575 
5576  /* Local variables */
5577  long k;
5578  double a1, a2, b1, b2, rg[25], rh[25], ri[25], rj[25];
5579  long nev;
5580  double area, area1, area2, area12;
5581  double erro12;
5582  long nrmax;
5583  double uflow;
5584  long iroff1, iroff2;
5585  double resas1, resas2, error1, error2, epmach, errbnd, centre;
5586  double errmax;
5587  long maxerr;
5588  double errsum;
5589 
5590  /* ***begin prologue dqawse */
5591  /* ***date written 800101 (yymmdd) */
5592  /* ***revision date 830518 (yymmdd) */
5593  /* ***category no. h2a2a1 */
5594  /* ***keywords automatic integrator, special-purpose, */
5595  /* algebraico-logarithmic end point singularities, */
5596  /* clenshaw-curtis method */
5597  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
5598  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
5599  /* ***purpose the routine calculates an approximation result to a given */
5600  /* definite integral i = integral of f*w over (a,b), */
5601  /* (where w shows a singular behaviour at the end points, */
5602  /* see parameter integr). */
5603  /* hopefully satisfying following claim for accuracy */
5604  /* fabs(i-result).le.max(epsabs,epsrel*fabs(i)). */
5605  /* ***description */
5606 
5607  /* integration of functions having algebraico-logarithmic */
5608  /* end point singularities */
5609  /* standard fortran subroutine */
5610  /* double precision version */
5611 
5612  /* parameters */
5613  /* on entry */
5614  /* f - double precision */
5615  /* function subprogram defining the integrand */
5616  /* function f(x). the actual name for f needs to be */
5617  /* declared e x t e r n a l in the driver program. */
5618 
5619  /* a - double precision */
5620  /* lower limit of integration */
5621 
5622  /* b - double precision */
5623  /* upper limit of integration, b.gt.a */
5624  /* if b.le.a, the routine will end with ier = 6. */
5625 
5626  /* alfa - double precision */
5627  /* parameter in the weight function, alfa.gt.(-1) */
5628  /* if alfa.le.(-1), the routine will end with */
5629  /* ier = 6. */
5630 
5631  /* beta - double precision */
5632  /* parameter in the weight function, beta.gt.(-1) */
5633  /* if beta.le.(-1), the routine will end with */
5634  /* ier = 6. */
5635 
5636  /* integr - long */
5637  /* indicates which weight function is to be used */
5638  /* = 1 (x-a)**alfa*(b-x)**beta */
5639  /* = 2 (x-a)**alfa*(b-x)**beta*log(x-a) */
5640  /* = 3 (x-a)**alfa*(b-x)**beta*log(b-x) */
5641  /* = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x) */
5642  /* if integr.lt.1 or integr.gt.4, the routine */
5643  /* will end with ier = 6. */
5644 
5645  /* epsabs - double precision */
5646  /* absolute accuracy requested */
5647  /* epsrel - double precision */
5648  /* relative accuracy requested */
5649  /* if epsabs.le.0 */
5650  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
5651  /* the routine will end with ier = 6. */
5652 
5653  /* limit - long */
5654  /* gives an upper bound on the number of subintervals */
5655  /* in the partition of (a,b), limit.ge.2 */
5656  /* if limit.lt.2, the routine will end with ier = 6. */
5657 
5658  /* on return */
5659  /* result - double precision */
5660  /* approximation to the integral */
5661 
5662  /* abserr - double precision */
5663  /* estimate of the modulus of the absolute error, */
5664  /* which should equal or exceed fabs(i-result) */
5665 
5666  /* neval - long */
5667  /* number of integrand evaluations */
5668 
5669  /* ier - long */
5670  /* ier = 0 normal and reliable termination of the */
5671  /* routine. it is assumed that the requested */
5672  /* accuracy has been achieved. */
5673  /* ier.gt.0 abnormal termination of the routine */
5674  /* the estimates for the integral and error */
5675  /* are less reliable. it is assumed that the */
5676  /* requested accuracy has not been achieved. */
5677  /* error messages */
5678  /* = 1 maximum number of subdivisions allowed */
5679  /* has been achieved. one can allow more */
5680  /* subdivisions by increasing the value of */
5681  /* limit. however, if this yields no */
5682  /* improvement, it is advised to analyze the */
5683  /* integrand in order to determine the */
5684  /* integration difficulties which prevent the */
5685  /* requested tolerance from being achieved. */
5686  /* in case of a jump discontinuity or a local */
5687  /* singularity of algebraico-logarithmic type */
5688  /* at one or more interior points of the */
5689  /* integration range, one should proceed by */
5690  /* splitting up the interval at these */
5691  /* points and calling the integrator on the */
5692  /* subranges. */
5693  /* = 2 the occurrence of roundoff error is */
5694  /* detected, which prevents the requested */
5695  /* tolerance from being achieved. */
5696  /* = 3 extremely bad integrand behaviour occurs */
5697  /* at some points of the integration */
5698  /* interval. */
5699  /* = 6 the input is invalid, because */
5700  /* b.le.a or alfa.le.(-1) or beta.le.(-1), or */
5701  /* integr.lt.1 or integr.gt.4, or */
5702  /* (epsabs.le.0 and */
5703  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
5704  /* or limit.lt.2. */
5705  /* result, abserr, neval, rlist(1), elist(1), */
5706  /* iord(1) and last are set to zero. alist(1) */
5707  /* and blist(1) are set to a and b */
5708  /* respectively. */
5709 
5710  /* alist - double precision */
5711  /* vector of dimension at least limit, the first */
5712  /* last elements of which are the left */
5713  /* end points of the subintervals in the partition */
5714  /* of the given integration range (a,b) */
5715 
5716  /* blist - double precision */
5717  /* vector of dimension at least limit, the first */
5718  /* last elements of which are the right */
5719  /* end points of the subintervals in the partition */
5720  /* of the given integration range (a,b) */
5721 
5722  /* rlist - double precision */
5723  /* vector of dimension at least limit,the first */
5724  /* last elements of which are the integral */
5725  /* approximations on the subintervals */
5726 
5727  /* elist - double precision */
5728  /* vector of dimension at least limit, the first */
5729  /* last elements of which are the moduli of the */
5730  /* absolute error estimates on the subintervals */
5731 
5732  /* iord - long */
5733  /* vector of dimension at least limit, the first k */
5734  /* of which are pointers to the error */
5735  /* estimates over the subintervals, so that */
5736  /* elist(iord(1)), ..., elist(iord(k)) with k = last */
5737  /* if last.le.(limit/2+2), and k = limit+1-last */
5738  /* otherwise form a decreasing sequence */
5739 
5740  /* last - long */
5741  /* number of subintervals actually produced in */
5742  /* the subdivision process */
5743 
5744  /* ***references (none) */
5745  /* ***routines called d1mach,dqc25s,dqmomo,dqpsrt */
5746  /* ***end prologue dqawse */
5747 
5748 
5749 
5750 
5751  /* list of major variables */
5752  /* ----------------------- */
5753 
5754  /* alist - list of left end points of all subintervals */
5755  /* considered up to now */
5756  /* blist - list of right end points of all subintervals */
5757  /* considered up to now */
5758  /* rlist(i) - approximation to the integral over */
5759  /* (alist(i),blist(i)) */
5760  /* elist(i) - error estimate applying to rlist(i) */
5761  /* maxerr - pointer to the interval with largest */
5762  /* error estimate */
5763  /* errmax - elist(maxerr) */
5764  /* area - sum of the integrals over the subintervals */
5765  /* errsum - sum of the errors over the subintervals */
5766  /* errbnd - requested accuracy max(epsabs,epsrel* */
5767  /* fabs(result)) */
5768  /* *****1 - variable for the left subinterval */
5769  /* *****2 - variable for the right subinterval */
5770  /* last - index for subdivision */
5771 
5772 
5773  /* machine dependent constants */
5774  /* --------------------------- */
5775 
5776  /* epmach is the largest relative spacing. */
5777  /* uflow is the smallest positive magnitude. */
5778 
5779  /* ***first executable statement dqawse */
5780  /* Parameter adjustments */
5781  --iord;
5782  --elist;
5783  --rlist;
5784  --blist;
5785  --alist__;
5786 
5787  /* Function Body */
5788  epmach = d1mach(c__4);
5789  uflow = d1mach(c__1);
5790 
5791  /* test on validity of parameters */
5792  /* ------------------------------ */
5793 
5794  *ier = 6;
5795  *neval = 0;
5796  *last = 0;
5797  rlist[1] = 0.;
5798  elist[1] = 0.;
5799  iord[1] = 0;
5800  *result = 0.;
5801  *abserr = 0.;
5802  /* Computing MAX */
5803  d__1 = epmach * 50.;
5804  if (*b <= *a || ( *epsabs == 0. && *epsrel < max(d__1,5e-29) ) || *alfa <=
5805  -1. || *beta <= -1. || *integr < 1 || *integr > 4 || *limit < 2) {
5806  goto L999;
5807  }
5808  *ier = 0;
5809 
5810  /* compute the modified chebyshev moments. */
5811 
5812  dqmomo_(alfa, beta, ri, rj, rg, rh, integr);
5813 
5814  /* integrate over the intervals (a,(a+b)/2) and ((a+b)/2,b). */
5815 
5816  centre = (*b + *a) * .5;
5817  dqc25s_(f, a, b, a, &centre, alfa, beta, ri, rj, rg, rh, &area1, &
5818  error1, &resas1, integr, &nev);
5819  *neval = nev;
5820  dqc25s_(f, a, b, &centre, b, alfa, beta, ri, rj, rg, rh, &area2, &
5821  error2, &resas2, integr, &nev);
5822  *last = 2;
5823  *neval += nev;
5824  *result = area1 + area2;
5825  *abserr = error1 + error2;
5826 
5827  /* test on accuracy. */
5828 
5829  /* Computing MAX */
5830  d__1 = *epsabs, d__2 = *epsrel * fabs(*result);
5831  errbnd = max(d__1,d__2);
5832 
5833  /* initialization */
5834  /* -------------- */
5835 
5836  if (error2 > error1) {
5837  goto L10;
5838  }
5839  alist__[1] = *a;
5840  alist__[2] = centre;
5841  blist[1] = centre;
5842  blist[2] = *b;
5843  rlist[1] = area1;
5844  rlist[2] = area2;
5845  elist[1] = error1;
5846  elist[2] = error2;
5847  goto L20;
5848 L10:
5849  alist__[1] = centre;
5850  alist__[2] = *a;
5851  blist[1] = *b;
5852  blist[2] = centre;
5853  rlist[1] = area2;
5854  rlist[2] = area1;
5855  elist[1] = error2;
5856  elist[2] = error1;
5857 L20:
5858  iord[1] = 1;
5859  iord[2] = 2;
5860  if (*limit == 2) {
5861  *ier = 1;
5862  }
5863  if (*abserr <= errbnd || *ier == 1) {
5864  goto L999;
5865  }
5866  errmax = elist[1];
5867  maxerr = 1;
5868  nrmax = 1;
5869  area = *result;
5870  errsum = *abserr;
5871  iroff1 = 0;
5872  iroff2 = 0;
5873 
5874  /* main do-loop */
5875  /* ------------ */
5876 
5877  i__1 = *limit;
5878  for (*last = 3; *last <= i__1; ++(*last)) {
5879 
5880  /* bisect the subinterval with largest error estimate. */
5881 
5882  a1 = alist__[maxerr];
5883  b1 = (alist__[maxerr] + blist[maxerr]) * .5;
5884  a2 = b1;
5885  b2 = blist[maxerr];
5886 
5887  dqc25s_(f, a, b, &a1, &b1, alfa, beta, ri, rj, rg, rh, &area1, &
5888  error1, &resas1, integr, &nev);
5889  *neval += nev;
5890  dqc25s_(f, a, b, &a2, &b2, alfa, beta, ri, rj, rg, rh, &area2, &
5891  error2, &resas2, integr, &nev);
5892  *neval += nev;
5893 
5894  /* improve previous approximations integral and error */
5895  /* and test for accuracy. */
5896 
5897  area12 = area1 + area2;
5898  erro12 = error1 + error2;
5899  errsum = errsum + erro12 - errmax;
5900  area = area + area12 - rlist[maxerr];
5901  if (*a == a1 || *b == b2) {
5902  goto L30;
5903  }
5904  if (resas1 == error1 || resas2 == error2) {
5905  goto L30;
5906  }
5907 
5908  /* test for roundoff error. */
5909 
5910  if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) < fabs(area12) * 1e-5 &&
5911  erro12 >= errmax * .99) {
5912  ++iroff1;
5913  }
5914  if (*last > 10 && erro12 > errmax) {
5915  ++iroff2;
5916  }
5917  L30:
5918  rlist[maxerr] = area1;
5919  rlist[*last] = area2;
5920 
5921  /* test on accuracy. */
5922 
5923  /* Computing MAX */
5924  d__1 = *epsabs, d__2 = *epsrel * fabs(area);
5925  errbnd = max(d__1,d__2);
5926  if (errsum <= errbnd) {
5927  goto L35;
5928  }
5929 
5930  /* set error flag in the case that the number of interval */
5931  /* bisections exceeds limit. */
5932 
5933  if (*last == *limit) {
5934  *ier = 1;
5935  }
5936 
5937 
5938  /* set error flag in the case of roundoff error. */
5939 
5940  if (iroff1 >= 6 || iroff2 >= 20) {
5941  *ier = 2;
5942  }
5943 
5944  /* set error flag in the case of bad integrand behaviour */
5945  /* at interior points of integration range. */
5946 
5947  /* Computing MAX */
5948  d__1 = fabs(a1), d__2 = fabs(b2);
5949  if (max(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
5950  {
5951  *ier = 3;
5952  }
5953 
5954  /* append the newly-created intervals to the list. */
5955 
5956  L35:
5957  if (error2 > error1) {
5958  goto L40;
5959  }
5960  alist__[*last] = a2;
5961  blist[maxerr] = b1;
5962  blist[*last] = b2;
5963  elist[maxerr] = error1;
5964  elist[*last] = error2;
5965  goto L50;
5966  L40:
5967  alist__[maxerr] = a2;
5968  alist__[*last] = a1;
5969  blist[*last] = b1;
5970  rlist[maxerr] = area2;
5971  rlist[*last] = area1;
5972  elist[maxerr] = error2;
5973  elist[*last] = error1;
5974 
5975  /* call subroutine dqpsrt to maintain the descending ordering */
5976  /* in the list of error estimates and select the subinterval */
5977  /* with largest error estimate (to be bisected next). */
5978 
5979  L50:
5980  dqpsrt_(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
5981  /* ***jump out of do-loop */
5982  if (*ier != 0 || errsum <= errbnd) {
5983  goto L70;
5984  }
5985  /* L60: */
5986  }
5987 
5988  /* compute final result. */
5989  /* --------------------- */
5990 
5991 L70:
5992  *result = 0.;
5993  i__1 = *last;
5994  for (k = 1; k <= i__1; ++k) {
5995  *result += rlist[k];
5996  /* L80: */
5997  }
5998  *abserr = errsum;
5999 L999:
6000  return;
6001 } /* dqawse_ */
6002 
6003 void dqaws_(const D_fp& f, const double *a, const double *b, const double *alfa,
6004  const double *beta, const long *integr, const double *epsabs,
6005  const double *epsrel, double *result, double *abserr,
6006  long *neval, long *ier, long *limit, const long *lenw, long *last,
6007  long *iwork, double *work)
6008 {
6009  long l1, l2, l3, lvl;
6010 
6011  /* ***begin prologue dqaws */
6012  /* ***date written 800101 (yymmdd) */
6013  /* ***revision date 830518 (yymmdd) */
6014  /* ***category no. h2a2a1 */
6015  /* ***keywords automatic integrator, special-purpose, */
6016  /* algebraico-logarithmic end-point singularities, */
6017  /* clenshaw-curtis, globally adaptive */
6018  /* ***author piessens,robert,appl. math. & progr. div. -k.u.leuven */
6019  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
6020  /* ***purpose the routine calculates an approximation result to a given */
6021  /* definite integral i = integral of f*w over (a,b), */
6022  /* (where w shows a singular behaviour at the end points */
6023  /* see parameter integr). */
6024  /* hopefully satisfying following claim for accuracy */
6025  /* fabs(i-result).le.max(epsabs,epsrel*fabs(i)). */
6026  /* ***description */
6027 
6028  /* integration of functions having algebraico-logarithmic */
6029  /* end point singularities */
6030  /* standard fortran subroutine */
6031  /* double precision version */
6032 
6033  /* parameters */
6034  /* on entry */
6035  /* f - double precision */
6036  /* function subprogram defining the integrand */
6037  /* function f(x). the actual name for f needs to be */
6038  /* declared e x t e r n a l in the driver program. */
6039 
6040  /* a - double precision */
6041  /* lower limit of integration */
6042 
6043  /* b - double precision */
6044  /* upper limit of integration, b.gt.a */
6045  /* if b.le.a, the routine will end with ier = 6. */
6046 
6047  /* alfa - double precision */
6048  /* parameter in the integrand function, alfa.gt.(-1) */
6049  /* if alfa.le.(-1), the routine will end with */
6050  /* ier = 6. */
6051 
6052  /* beta - double precision */
6053  /* parameter in the integrand function, beta.gt.(-1) */
6054  /* if beta.le.(-1), the routine will end with */
6055  /* ier = 6. */
6056 
6057  /* integr - long */
6058  /* indicates which weight function is to be used */
6059  /* = 1 (x-a)**alfa*(b-x)**beta */
6060  /* = 2 (x-a)**alfa*(b-x)**beta*log(x-a) */
6061  /* = 3 (x-a)**alfa*(b-x)**beta*log(b-x) */
6062  /* = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x) */
6063  /* if integr.lt.1 or integr.gt.4, the routine */
6064  /* will end with ier = 6. */
6065 
6066  /* epsabs - double precision */
6067  /* absolute accuracy requested */
6068  /* epsrel - double precision */
6069  /* relative accuracy requested */
6070  /* if epsabs.le.0 */
6071  /* and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), */
6072  /* the routine will end with ier = 6. */
6073 
6074  /* on return */
6075  /* result - double precision */
6076  /* approximation to the integral */
6077 
6078  /* abserr - double precision */
6079  /* estimate of the modulus of the absolute error, */
6080  /* which should equal or exceed fabs(i-result) */
6081 
6082  /* neval - long */
6083  /* number of integrand evaluations */
6084 
6085  /* ier - long */
6086  /* ier = 0 normal and reliable termination of the */
6087  /* routine. it is assumed that the requested */
6088  /* accuracy has been achieved. */
6089  /* ier.gt.0 abnormal termination of the routine */
6090  /* the estimates for the integral and error */
6091  /* are less reliable. it is assumed that the */
6092  /* requested accuracy has not been achieved. */
6093  /* error messages */
6094  /* ier = 1 maximum number of subdivisions allowed */
6095  /* has been achieved. one can allow more */
6096  /* subdivisions by increasing the value of */
6097  /* limit (and taking the according dimension */
6098  /* adjustments into account). however, if */
6099  /* this yields no improvement it is advised */
6100  /* to analyze the integrand, in order to */
6101  /* determine the integration difficulties */
6102  /* which prevent the requested tolerance from */
6103  /* being achieved. in case of a jump */
6104  /* discontinuity or a local singularity */
6105  /* of algebraico-logarithmic type at one or */
6106  /* more interior points of the integration */
6107  /* range, one should proceed by splitting up */
6108  /* the interval at these points and calling */
6109  /* the integrator on the subranges. */
6110  /* = 2 the occurrence of roundoff error is */
6111  /* detected, which prevents the requested */
6112  /* tolerance from being achieved. */
6113  /* = 3 extremely bad integrand behaviour occurs */
6114  /* at some points of the integration */
6115  /* interval. */
6116  /* = 6 the input is invalid, because */
6117  /* b.le.a or alfa.le.(-1) or beta.le.(-1) or */
6118  /* or integr.lt.1 or integr.gt.4 or */
6119  /* (epsabs.le.0 and */
6120  /* epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) */
6121  /* or limit.lt.2 or lenw.lt.limit*4. */
6122  /* result, abserr, neval, last are set to */
6123  /* zero. except when lenw or limit is invalid */
6124  /* iwork(1), work(limit*2+1) and */
6125  /* work(limit*3+1) are set to zero, work(1) */
6126  /* is set to a and work(limit+1) to b. */
6127 
6128  /* dimensioning parameters */
6129  /* limit - long */
6130  /* dimensioning parameter for iwork */
6131  /* limit determines the maximum number of */
6132  /* subintervals in the partition of the given */
6133  /* integration interval (a,b), limit.ge.2. */
6134  /* if limit.lt.2, the routine will end with ier = 6. */
6135 
6136  /* lenw - long */
6137  /* dimensioning parameter for work */
6138  /* lenw must be at least limit*4. */
6139  /* if lenw.lt.limit*4, the routine will end */
6140  /* with ier = 6. */
6141 
6142  /* last - long */
6143  /* on return, last equals the number of */
6144  /* subintervals produced in the subdivision process, */
6145  /* which determines the significant number of */
6146  /* elements actually in the work arrays. */
6147 
6148  /* work arrays */
6149  /* iwork - long */
6150  /* vector of dimension limit, the first k */
6151  /* elements of which contain pointers */
6152  /* to the error estimates over the subintervals, */
6153  /* such that work(limit*3+iwork(1)), ..., */
6154  /* work(limit*3+iwork(k)) form a decreasing */
6155  /* sequence with k = last if last.le.(limit/2+2), */
6156  /* and k = limit+1-last otherwise */
6157 
6158  /* work - double precision */
6159  /* vector of dimension lenw */
6160  /* on return */
6161  /* work(1), ..., work(last) contain the left */
6162  /* end points of the subintervals in the */
6163  /* partition of (a,b), */
6164  /* work(limit+1), ..., work(limit+last) contain */
6165  /* the right end points, */
6166  /* work(limit*2+1), ..., work(limit*2+last) */
6167  /* contain the integral approximations over */
6168  /* the subintervals, */
6169  /* work(limit*3+1), ..., work(limit*3+last) */
6170  /* contain the error estimates. */
6171 
6172  /* ***references (none) */
6173  /* ***routines called dqawse,xerror */
6174  /* ***end prologue dqaws */
6175 
6176 
6177 
6178 
6179  /* check validity of limit and lenw. */
6180 
6181  /* ***first executable statement dqaws */
6182  /* Parameter adjustments */
6183  --iwork;
6184  --work;
6185 
6186  /* Function Body */
6187  *ier = 6;
6188  *neval = 0;
6189  *last = 0;
6190  *result = 0.;
6191  *abserr = 0.;
6192  if (*limit < 2 || *lenw < *limit << 2) {
6193  goto L10;
6194  }
6195 
6196  /* prepare call for dqawse. */
6197 
6198  l1 = *limit + 1;
6199  l2 = *limit + l1;
6200  l3 = *limit + l2;
6201 
6202  dqawse_(f, a, b, alfa, beta, integr, epsabs, epsrel, limit, result,
6203  abserr, neval, ier, &work[1], &work[l1], &work[l2], &work[l3], &
6204  iwork[1], last);
6205 
6206  /* call error handler if necessary. */
6207 
6208  lvl = 0;
6209 L10:
6210  if (*ier == 6) {
6211  lvl = 1;
6212  }
6213  if (*ier != 0) {
6214  xerror_("abnormal return from dqaws", &c__26, ier, &lvl, 26);
6215  }
6216  return;
6217 } /* dqaws_ */
6218 
6219 void dqc25c_(const D_fp& f, const double *a, const double *b, const double *c__,
6220  double *result, double *abserr, long *krul, long *neval)
6221 {
6222  /* Initialized data */
6223 
6224  double x[11] = { .991444861373810411144557526928563,
6225  .965925826289068286749743199728897,
6226  .923879532511286756128183189396788,
6227  .866025403784438646763723170752936,
6228  .793353340291235164579776961501299,
6229  .707106781186547524400844362104849,
6230  .608761429008720639416097542898164,.5,
6231  .382683432365089771728459984030399,
6232  .258819045102520762348898837624048,
6233  .130526192220051591548406227895489 };
6234 
6235  /* System generated locals */
6236  double d__1;
6237 
6238  /* Local variables */
6239  long i__, k;
6240  double u, p2, p3, p4, cc;
6241  long kp;
6242  double ak22, fval[25], res12, res24;
6243  long isym;
6244  double amom0, amom1, amom2, cheb12[13], cheb24[25], hlgth,
6245  centr;
6246  double resabs, resasc;
6247 
6248  /* ***begin prologue dqc25c */
6249  /* ***date written 810101 (yymmdd) */
6250  /* ***revision date 830518 (yymmdd) */
6251  /* ***category no. h2a2a2,j4 */
6252  /* ***keywords 25-point clenshaw-curtis integration */
6253  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
6254  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
6255  /* ***purpose to compute i = integral of f*w over (a,b) with */
6256  /* error estimate, where w(x) = 1/(x-c) */
6257  /* ***description */
6258 
6259  /* integration rules for the computation of cauchy */
6260  /* principal value integrals */
6261  /* standard fortran subroutine */
6262  /* double precision version */
6263 
6264  /* parameters */
6265  /* f - double precision */
6266  /* function subprogram defining the integrand function */
6267  /* f(x). the actual name for f needs to be declared */
6268  /* e x t e r n a l in the driver program. */
6269 
6270  /* a - double precision */
6271  /* left end point of the integration interval */
6272 
6273  /* b - double precision */
6274  /* right end point of the integration interval, b.gt.a */
6275 
6276  /* c - double precision */
6277  /* parameter in the weight function */
6278 
6279  /* result - double precision */
6280  /* approximation to the integral */
6281  /* result is computed by using a generalized */
6282  /* clenshaw-curtis method if c lies within ten percent */
6283  /* of the integration interval. in the other case the */
6284  /* 15-point kronrod rule obtained by optimal addition */
6285  /* of abscissae to the 7-point gauss rule, is applied. */
6286 
6287  /* abserr - double precision */
6288  /* estimate of the modulus of the absolute error, */
6289  /* which should equal or exceed fabs(i-result) */
6290 
6291  /* krul - long */
6292  /* key which is decreased by 1 if the 15-polong */
6293  /* gauss-kronrod scheme has been used */
6294 
6295  /* neval - long */
6296  /* number of integrand evaluations */
6297 
6298  /* ....................................................................... */
6299  /* ***references (none) */
6300  /* ***routines called dqcheb,dqk15w,dqwgtc */
6301  /* ***end prologue dqc25c */
6302 
6303 
6304 
6305 
6306  /* the vector x contains the values cos(k*pi/24), */
6307  /* k = 1, ..., 11, to be used for the chebyshev series */
6308  /* expansion of f */
6309 
6310 
6311  /* list of major variables */
6312  /* ---------------------- */
6313  /* fval - value of the function f at the points */
6314  /* cos(k*pi/24), k = 0, ..., 24 */
6315  /* cheb12 - chebyshev series expansion coefficients, */
6316  /* for the function f, of degree 12 */
6317  /* cheb24 - chebyshev series expansion coefficients, */
6318  /* for the function f, of degree 24 */
6319  /* res12 - approximation to the integral corresponding */
6320  /* to the use of cheb12 */
6321  /* res24 - approximation to the integral corresponding */
6322  /* to the use of cheb24 */
6323  /* dqwgtc - external function subprogram defining */
6324  /* the weight function */
6325  /* hlgth - half-length of the interval */
6326  /* centr - mid point of the interval */
6327 
6328 
6329  /* check the position of c. */
6330 
6331  /* ***first executable statement dqc25c */
6332  cc = (2. * *c__ - *b - *a) / (*b - *a);
6333  if (fabs(cc) < 1.1) {
6334  goto L10;
6335  }
6336 
6337  /* apply the 15-point gauss-kronrod scheme. */
6338 
6339  --(*krul);
6340  dqk15w_(f, dqwgtc_, c__, &p2, &p3, &p4, &kp, a, b, result,
6341  abserr, &resabs, &resasc);
6342  *neval = 15;
6343  if (resasc == *abserr) {
6344  ++(*krul);
6345  }
6346  goto L50;
6347 
6348  /* use the generalized clenshaw-curtis method. */
6349 
6350 L10:
6351  hlgth = (*b - *a) * .5;
6352  centr = (*b + *a) * .5;
6353  *neval = 25;
6354  d__1 = hlgth + centr;
6355  fval[0] = f(d__1) * .5;
6356  fval[12] = f(centr);
6357  d__1 = centr - hlgth;
6358  fval[24] = f(d__1) * .5;
6359  for (i__ = 2; i__ <= 12; ++i__) {
6360  u = hlgth * x[i__ - 2];
6361  isym = 26 - i__;
6362  d__1 = u + centr;
6363  fval[i__ - 1] = f(d__1);
6364  d__1 = centr - u;
6365  fval[isym - 1] = f(d__1);
6366  /* L20: */
6367  }
6368 
6369  /* compute the chebyshev series expansion. */
6370 
6371  dqcheb_(x, fval, cheb12, cheb24);
6372 
6373  /* the modified chebyshev moments are computed by forward */
6374  /* recursion, using amom0 and amom1 as starting values. */
6375 
6376  amom0 = log((d__1 = (1. - cc) / (cc + 1.), fabs(d__1)));
6377  amom1 = cc * amom0 + 2.;
6378  res12 = cheb12[0] * amom0 + cheb12[1] * amom1;
6379  res24 = cheb24[0] * amom0 + cheb24[1] * amom1;
6380  for (k = 3; k <= 13; ++k) {
6381  amom2 = cc * 2. * amom1 - amom0;
6382  ak22 = (double) ((k - 2) * (k - 2));
6383  if (k / 2 << 1 == k) {
6384  amom2 -= 4. / (ak22 - 1.);
6385  }
6386  res12 += cheb12[k - 1] * amom2;
6387  res24 += cheb24[k - 1] * amom2;
6388  amom0 = amom1;
6389  amom1 = amom2;
6390  /* L30: */
6391  }
6392  for (k = 14; k <= 25; ++k) {
6393  amom2 = cc * 2. * amom1 - amom0;
6394  ak22 = (double) ((k - 2) * (k - 2));
6395  if (k / 2 << 1 == k) {
6396  amom2 -= 4. / (ak22 - 1.);
6397  }
6398  res24 += cheb24[k - 1] * amom2;
6399  amom0 = amom1;
6400  amom1 = amom2;
6401  /* L40: */
6402  }
6403  *result = res24;
6404  *abserr = (d__1 = res24 - res12, fabs(d__1));
6405 L50:
6406  return;
6407 } /* dqc25c_ */
6408 
6409 void dqc25f_(const D_fp& f, const double *a, const double *b, const double
6410  *omega, const long *integr, const long *nrmom, const long *maxp1,
6411  const long *ksave, double *result, double *abserr, long *neval,
6412  double *resabs, double *resasc, long *momcom, double *
6413  chebmo)
6414 {
6415  /* Initialized data */
6416 
6417  double x[11] = { .991444861373810411144557526928563,
6418  .965925826289068286749743199728897,
6419  .923879532511286756128183189396788,
6420  .866025403784438646763723170752936,
6421  .793353340291235164579776961501299,
6422  .707106781186547524400844362104849,
6423  .608761429008720639416097542898164,.5,
6424  .382683432365089771728459984030399,
6425  .258819045102520762348898837624048,
6426  .130526192220051591548406227895489 };
6427 
6428  /* System generated locals */
6429  int chebmo_dim1, chebmo_offset, i__1;
6430  double d__1, d__2;
6431 
6432  /* Local variables */
6433  double d__[25];
6434  long i__, j, k, m=0;
6435  double v[28], d1[25], d2[25], p2, p3, p4, ac, an, as, an2, ass,
6436  par2, conc, asap, par22, fval[25], estc, cons;
6437  long iers;
6438  double ests;
6439  long isym, noeq1;
6440  double cheb12[13], cheb24[25], resc12, resc24, hlgth, centr;
6441  double ress12, ress24, oflow;
6442  long noequ;
6443  double cospar;
6444  double parint, sinpar;
6445 
6446  /* ***begin prologue dqc25f */
6447  /* ***date written 810101 (yymmdd) */
6448  /* ***revision date 830518 (yymmdd) */
6449  /* ***category no. h2a2a2 */
6450  /* ***keywords integration rules for functions with cos or sin */
6451  /* factor, clenshaw-curtis, gauss-kronrod */
6452  /* ***author piessens,robert,appl. math. & progr. div. - k.u.leuven */
6453  /* de doncker,elise,appl. math. & progr. div. - k.u.leuven */
6454  /* ***purpose to compute the integral i=integral of f(x) over (a,b) */
6455  /* where w(x) = cos(omega*x) or w(x)=sin(omega*x) and to */
6456  /* compute j = integral of fabs(f) over (a,b). for small value */
6457  /* of omega or small intervals (a,b) the 15-point gauss-kronro */
6458  /* rule is used. otherwise a generalized clenshaw-curtis */
6459  /* method is used. */
6460  /* ***description */
6461 
6462  /* integration rules for functions with cos or sin factor */
6463  /* standard fortran subroutine */
6464  /* double precision version */
6465 
6466  /* parameters */
6467  /* on entry */
6468  /* f - double precision */
6469  /* function subprogram defining the integrand */
6470  /* function f(x). the actual name for f needs to */
6471  /* be declared e x t e r n a l in the calling program. */
6472 
6473  /* a - double precision */
6474  /* lower limit of integration */
6475 
6476  /* b - double precision */
6477  /* upper limit of integration */
6478 
6479  /* omega - double precision */
6480  /* parameter in the weight function */
6481 
6482  /* integr - long */
6483  /* indicates which weight function is to be used */
6484  /* integr = 1 w(x) = cos(omega*x) */
6485  /* integr = 2 w(x) = sin(omega*x) */
6486 
6487  /* nrmom - long */
6488  /* the length of interval (a,b) is equal to the length */
6489  /* of the original integration interval divided by */
6490  /* 2**nrmom (we suppose that the routine is used in an */
6491  /* adaptive integration process, otherwise set */
6492  /* nrmom = 0). nrmom must be zero at the first call. */
6493 
6494  /* maxp1 - long */
6495  /* gives an upper bound on the number of chebyshev */
6496  /* moments which can be stored, i.e. for the */
6497  /* intervals of lengths fabs(bb-aa)*2**(-l), */
6498  /* l = 0,1,2, ..., maxp1-2. */
6499 
6500  /* ksave - long */
6501  /* key which is one when the moments for the */
6502  /* current interval have been computed */
6503 
6504  /* on return */
6505  /* result - double precision */
6506  /* approximation to the integral i */
6507 
6508  /* abserr - double precision */
6509  /* estimate of the modulus of the absolute */
6510  /* error, which should equal or exceed fabs(i-result) */
6511 
6512  /* neval - long */
6513  /* number of integrand evaluations */
6514 
6515  /* resabs - double precision */
6516  /* approximation to the integral j */
6517 
6518  /* resasc - double precision */
6519  /* approximation to the integral of fabs(f-i/(b-a)) */
6520 
6521  /* on entry and return */
6522  /* momcom - long */
6523  /* for each interval length we need to compute the */
6524  /* chebyshev moments. momcom counts the number of */
6525  /* intervals for which these moments have already been */
6526  /* computed. if nrmom.lt.momcom or ksave = 1, the */
6527  /* chebyshev moments for the interval (a,b) have */
6528  /* already been computed and stored, otherwise we */
6529  /* compute them and we increase momcom. */
6530 
6531  /* chebmo - double precision */
6532  /* array of dimension at least (maxp1,25) containing */
6533  /* the modified chebyshev moments for the first momcom */
6534  /* momcom interval lengths */
6535 
6536  /* ...................................................................... */
6537  /* ***references (none) */
6538  /* ***routines called d1mach,dgtsl,dqcheb,dqk15w,dqwgtf */
6539  /* ***end prologue dqc25f */
6540 
6541 
6542 
6543 
6544  /* the vector x contains the values cos(k*pi/24) */
6545  /* k = 1, ...,11, to be used for the chebyshev expansion of f */
6546 
6547  /* Parameter adjustments */
6548  chebmo_dim1 = *maxp1;
6549  chebmo_offset = 1 + chebmo_dim1;
6550  chebmo -= chebmo_offset;
6551 
6552  /* Function Body */
6553 
6554  /* list of major variables */
6555  /* ----------------------- */
6556 
6557  /* centr - mid point of the integration interval */
6558  /* hlgth - half-length of the integration interval */
6559  /* fval - value of the function f at the points */
6560  /* (b-a)*0.5*cos(k*pi/12) + (b+a)*0.5, k = 0, ..., 24 */
6561  /* cheb12 - coefficients of the chebyshev series expansion */
6562  /* of degree 12, for the function f, in the */
6563  /* interval (a,b) */
6564  /* cheb24 - coefficients of the chebyshev series expansion */
6565  /* of degree 24, for the function f, in the */
6566  /* interval (a,b) */
6567  /* resc12 - approximation to the integral of */
6568  /* cos(0.5*(b-a)*omega*x)*f(0.5*(b-a)*x+0.5*(b+a)) */
6569  /* over (-1,+1), using the chebyshev series */
6570  /* expansion of degree 12 */
6571  /* resc24 - approximation to the same integral, using the */
6572  /* chebyshev series expansion of degree 24 */
6573  /* ress12 - the analogue of resc12 for the sine */
6574  /* ress24 - the analogue of resc24 for the sine */
6575 
6576 
6577  /* machine dependent constant */
6578  /* -------------------------- */
6579 
6580  /* oflow is the largest positive magnitude. */
6581 
6582  /* ***first executable statement dqc25f */
6583  oflow = d1mach(c__2);
6584 
6585  centr = (*b + *a) * .5;
6586  hlgth = (*b - *a) * .5;
6587  parint = *omega * hlgth;
6588 
6589  /* compute the integral using the 15-point gauss-kronrod */
6590  /* formula if the value of the parameter in the integrand */
6591  /* is small. */
6592 
6593  if (fabs(parint) > 2.) {
6594  goto L10;
6595  }
6596  dqk15w_(f, dqwgtf_, omega, &p2, &p3, &p4, integr, a, b,
6597  result, abserr, resabs, resasc);
6598  *neval = 15;
6599  goto L170;
6600 
6601  /* compute the integral using the generalized clenshaw- */
6602  /* curtis method. */
6603 
6604 L10:
6605  conc = hlgth * cos(centr * *omega);
6606  cons = hlgth * sin(centr * *omega);
6607  *resasc = oflow;
6608  *neval = 25;
6609 
6610  /* check whether the chebyshev moments for this interval */
6611  /* have already been computed. */
6612 
6613  if (*nrmom < *momcom || *ksave == 1) {
6614  goto L120;
6615  }
6616 
6617  /* compute a new set of chebyshev moments. */
6618 
6619  m = *momcom + 1;
6620  par2 = parint * parint;
6621  par22 = par2 + 2.;
6622  sinpar = sin(parint);
6623  cospar = cos(parint);
6624 
6625  /* compute the chebyshev moments with respect to cosine. */
6626 
6627  v[0] = sinpar * 2. / parint;
6628  v[1] = (cospar * 8. + (par2 + par2 - 8.) * sinpar / parint) / par2;
6629  v[2] = ((par2 - 12.) * 32. * cospar + ((par2 - 80.) * par2 + 192.) * 2. *
6630  sinpar / parint) / (par2 * par2);
6631  ac = cospar * 8.;
6632  as = parint * 24. * sinpar;
6633  if (fabs(parint) > 24.) {
6634  goto L30;
6635  }
6636 
6637  /* compute the chebyshev moments as the solutions of a */
6638  /* boundary value problem with 1 initial value (v(3)) and 1 */
6639  /* end value (computed using an asymptotic formula). */
6640 
6641  noequ = 25;
6642  noeq1 = noequ - 1;
6643  an = 6.;
6644  i__1 = noeq1;
6645  for (k = 1; k <= i__1; ++k) {
6646  an2 = an * an;
6647  d__[k - 1] = (an2 - 4.) * -2. * (par22 - an2 - an2);
6648  d2[k - 1] = (an - 1.) * (an - 2.) * par2;
6649  d1[k] = (an + 3.) * (an + 4.) * par2;
6650  v[k + 2] = as - (an2 - 4.) * ac;
6651  an += 2.;
6652  /* L20: */
6653  }
6654  an2 = an * an;
6655  d__[noequ - 1] = (an2 - 4.) * -2. * (par22 - an2 - an2);
6656  v[noequ + 2] = as - (an2 - 4.) * ac;
6657  v[3] -= par2 * 56. * v[2];
6658  ass = parint * sinpar;
6659  asap = (((((par2 * 210. - 1.) * cospar - (par2 * 105. - 63.) * ass) / an2
6660  - (1. - par2 * 15.) * cospar + ass * 15.) / an2 - cospar + ass *
6661  3.) / an2 - cospar) / an2;
6662  v[noequ + 2] -= asap * 2. * par2 * (an - 1.) * (an - 2.);
6663 
6664  /* solve the tridiagonal system by means of gaussian */
6665  /* elimination with partial pivoting. */
6666 
6667  /* *** call to dgtsl must be replaced by call to */
6668  /* *** double precision version of linpack routine sgtsl */
6669 
6670  dgtsl_(&noequ, d1, d__, d2, &v[3], &iers);
6671  goto L50;
6672 
6673  /* compute the chebyshev moments by means of forward */
6674  /* recursion. */
6675 
6676 L30:
6677  an = 4.;
6678  for (i__ = 4; i__ <= 13; ++i__) {
6679  an2 = an * an;
6680  v[i__ - 1] = ((an2 - 4.) * ((par22 - an2 - an2) * 2. * v[i__ - 2] -
6681  ac) + as - par2 * (an + 1.) * (an + 2.) * v[i__ - 3]) / (par2
6682  * (an - 1.) * (an - 2.));
6683  an += 2.;
6684  /* L40: */
6685  }
6686 L50:
6687  for (j = 1; j <= 13; ++j) {
6688  chebmo[m + ((j << 1) - 1) * chebmo_dim1] = v[j - 1];
6689  /* L60: */
6690  }
6691 
6692  /* compute the chebyshev moments with respect to sine. */
6693 
6694  v[0] = (sinpar - parint * cospar) * 2. / par2;
6695  v[1] = (18. - 48. / par2) * sinpar / par2 + (48. / par2 - 2.) * cospar /
6696  parint;
6697  ac = parint * -24. * cospar;
6698  as = sinpar * -8.;
6699  if (fabs(parint) > 24.) {
6700  goto L80;
6701  }
6702 
6703  /* compute the chebyshev moments as the solutions of a boundary */
6704  /* value problem with 1 initial value (v(2)) and 1 end value */
6705  /* (computed using an asymptotic formula). */
6706 
6707  an = 5.;
6708  i__1 = noeq1;
6709  for (k = 1; k <= i__1; ++k) {
6710  an2 = an * an;
6711  d__[k - 1] = (an2 - 4.) * -2. * (par22 - an2 - an2);
6712  d2[k - 1] = (an - 1.) * (an - 2.) * par2;
6713  d1[k] = (an + 3.) * (an + 4.) * par2;
6714  v[k + 1] = ac + (an2 - 4.) * as;
6715  an += 2.;
6716  /* L70: */
6717  }
6718  an2 = an * an;
6719  d__[noequ - 1] = (an2 - 4.) * -2. * (par22 - an2 - an2);
6720  v[noequ + 1] = ac + (an2 - 4.) * as;
6721  v[2] -= par2 * 42. * v[1];
6722  ass = parint * cospar;
6723  asap = (((((par2 * 105. - 63.) * ass + (par2 * 210. - 1.) * sinpar) / an2
6724  + (par2 * 15. - 1.) * sinpar - ass * 15.) / an2 - ass * 3. -
6725  sinpar) / an2 - sinpar) / an2;
6726  v[noequ + 1] -= asap * 2. * par2 * (an - 1.) * (an - 2.);
6727 
6728  /* solve the tridiagonal system by means of gaussian */
6729  /* elimination with partial pivoting. */
6730 
6731  /* *** call to dgtsl must be replaced by call to */
6732  /* *** double precision version of linpack routine sgtsl */
6733 
6734  dgtsl_(&noequ, d1, d__, d2, &v[2], &iers);
6735  goto L100;
6736 
6737  /* compute the chebyshev moments by means of forward recursion. */
6738 
6739 L80:
6740  an = 3.;
6741  for (i__ = 3; i__ <= 12; ++i__) {
6742  an2 = an * an;
6743  v[i__ - 1] = ((an2 - 4.) * ((par22 - an2 - an2) * 2. * v[i__ - 2] +
6744  as) + ac - par2 * (an + 1.) * (an + 2.) * v[i__ - 3]) / (par2