The R Project SVN – R – Diff – /trunk/src/appl/integrate.c

 /*
  *  R : A Computer Language for Statistical Data Analysis
- *  Copyright (C) 2001  the R Development Core Team
+ *  Copyright (C) 2001-2002  the R Development Core Team
+ *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation; either version 2 of the License, or
  *  (at your option) any later version.
  *  GNU General Public License for more details.
+ *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ * C backend of R's integrate() --- via
+ *	 .External("call_dqags", ...) -> Rdqags()  -- for finite     interval
+ *	 .External("call_dqagi", ...) -> Rdqagi()  -- for indefinite interval
- */
+*/
 #include <Rinternals.h>
-#include <R_ext/Boolean.h>
-#ifdef Macintosh
-#include <fp.h>
-#else
-#include <math.h>
-#endif
 #include <Rmath.h>
+#include <R_ext/Applic.h> /* exporting the API , particularly */
+/*--- typedef void integr_fn(double *x, int n, void *ex) ---
+ * vectorizing function   f(x[1:n], ...) -> x[]  {overwriting x[]}.
+ * Vectorization can be used to speed up the integrand
+ * instead of calling it  n  times.
+*/
+/* Only these two are called via .External(.) :*/
+SEXP call_dqags(SEXP args);
+SEXP call_dqagi(SEXP args);
 typedef struct int_struct
+{
     SEXP f;    /* function */
     SEXP env;  /* where to evaluate the calls */
 } int_struct, *IntStruct;
-static void Rdqags(IntStruct IS, double *lower, double *upper,
-		   double *epsabs, double *epsrel, double *result,
-		   double *abserr, int *neval, int *ier, int *limit,
-		   int *lenw, int *last, int *iwork, double *work);
-static void Rdqagi(IntStruct IS, double *bound, int *inf,
-		   double *epsabs, double *epsrel, double *result,
-		   double *abserr, int *neval, int *ier, int *limit,
-		   int *lenw, int *last, int *iwork, double *work);
+/* This is *the* ``integr_fn f'' used when called from R : */
-static void Rintfn(double *x, int n, IntStruct IS)
+static void Rintfn(double *x, int n, void *ex)
+{
     SEXP args, resultsxp;
     int i;
+    IntStruct IS = (IntStruct) ex;
     PROTECT(args = allocVector(REALSXP, n));
     for(i = 0; i < n; i++) REAL(args)[i] = x[i];
     PROTECT(resultsxp = eval(lang2(IS->f , args), IS->env));
     if(length(resultsxp) != n)
 	error("evaluation of function gave a result of wrong length");
     for(i = 0; i < n; i++) {
 	x[i] = REAL(resultsxp)[i];
 	if(!R_FINITE(x[i]))
 	    error("non-finite function value");
+    }
     UNPROTECT(2);
     return;
+}
 SEXP call_dqags(SEXP args)
+{
     int_struct is;
     SEXP ans, ansnames;
     double lower, upper, epsabs, epsrel, result, abserr, *work;
     limit = asInteger(CAR(args)); args = CDR(args);
     lenw = 4 * limit;
     iwork = (int *) R_alloc(limit, sizeof(int));
     work = (double *) R_alloc(lenw, sizeof(double));
+    Rdqags(Rintfn, (void*)&is,
-    Rdqags(&is, &lower, &upper, &epsabs, &epsrel, &result,
+	   &lower, &upper, &epsabs, &epsrel, &result,
 	   &abserr, &neval, &ier, &limit, &lenw, &last, iwork, work);
     PROTECT(ans = allocVector(VECSXP, 4));
     PROTECT(ansnames = allocVector(STRSXP, 4));
     SET_STRING_ELT(ansnames, 0, mkChar("value"));
     limit = asInteger(CAR(args)); args = CDR(args);
     lenw = 4 * limit;
     iwork = (int *) R_alloc(limit, sizeof(int));
     work = (double *) R_alloc(lenw, sizeof(double));
-    Rdqagi(&is, &bound,&inf,&epsabs,&epsrel,&result,
+    Rdqagi(Rintfn, (void*)&is, &bound,&inf,&epsabs,&epsrel,&result,
 	   &abserr,&neval,&ier,&limit,&lenw,&last,iwork,work);
     PROTECT(ans = allocVector(VECSXP, 4));
     PROTECT(ansnames = allocVector(STRSXP, 4));
     SET_STRING_ELT(ansnames, 0, mkChar("value"));
     setAttrib(ans, R_NamesSymbol, ansnames);
     UNPROTECT(2);
     return ans;
+}
 /* f2c-ed translations + modifications of QUADPACK functions from here down */
+static void rdqagie(integr_fn f, void *ex,
-static void rdqagie(IntStruct IS, double *, int *, double * , double *, int *,
+		    double *, int *, double * , double *, int *,
 		    double *, double *, int *,
 		    int *, double *, double *, double *, double *,
 		    int *, int *);
+static void rdqk15i(integr_fn f, void *ex,
-static void  rdqk15i(IntStruct IS, double *, int *, double * , double *,
+		    double *, int *, double * , double *,
-		     double *, double *, double *, double *);
+		    double *, double *, double *, double *);
+static void rdqagse(integr_fn f, void *ex, double *, double *,
+		    double *, double *, int *, double *, double *,
+		    int *, int *, double *, double *, double *,
+		    double *, int *, int *);
+static void rdqk21(integr_fn f, void *ex,
-static void  rdqelg(int *, double *, double *, double *, double *, int *);
+		   double *, double *, double *, double *, double *, double *);
 static void rdqpsrt(int *, int *, int *, double *, double *, int *, int *);
-static void  rdqagse(IntStruct IS, double *, double *,
-		     double *, double *, int *, double *, double *,
-		     int *, int *, double *, double *, double *,
-		     double *, int *, int *);
-static void rdqk21(IntStruct IS, double *, double *, double *, double *, double *, double *);
 static void rdqelg(int *, double *, double *, double *, double *, int *);
 /* Table of constant values */
 static double c_b6 = 0.;
 static double c_b7 = 1.;
-static
-void Rdqagi(IntStruct IS, double *bound, int *inf, double *
+void Rdqagi(integr_fn f, void *ex, double *bound, int *inf,
+	    double *epsabs, double *epsrel,
-	    epsabs, double *epsrel, double *result, double *abserr,
+	    double *result, double *abserr, int *neval, int *ier,
-	    int *neval, int *ier, int *limit, int *lenw, int *
+	    int *limit, int *lenw, int *last,
-	    last, int *iwork, double *work)
+	    int *iwork, double *work)
+{
     int l1, l2, l3;
 /*
 ***begin prologue  dqagi
 ***purpose  the routine calculates an approximation result to a given
            integral   i = integral of f over (bound,+infinity)
            or i = integral of f over (-infinity,bound)
            or i = integral of f over (-infinity,+infinity)
            hopefully satisfying following claim for accuracy
-           abs(i-result).le.max(epsabs,epsrel*abs(i)).
+           abs(i-result) <= max(epsabs,epsrel*abs(i)).
 ***description
        integration over infinite intervals
        standard fortran subroutine
            epsabs - double precision
                     absolute accuracy requested
            epsrel - double precision
                     relative accuracy requested
-                    if  epsabs.le.0
+                    if  epsabs <= 0
-                    and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+                    and epsrel < max(50*rel.mach.acc.,0.5d-28),
                     the routine will end with ier = 6.
         on return
            result - double precision
            ier    - int
                     ier = 0 normal and reliable termination of the
                             routine. it is assumed that the requested
                             accuracy has been achieved.
-                  - ier.gt.0 abnormal termination of the routine. the
+                  - ier > 0 abnormal termination of the routine. the
                             estimates for result and error are less
                             reliable. it is assumed that the requested
                             accuracy has not been achieved.
            error messages
                     ier = 1 maximum number of subdivisions allowed
                         = 5 the integral is probably divergent, or
                             slowly convergent. it must be noted that
                             divergence can occur with any other value
                             of ier.
                         = 6 the input is invalid, because
-                            (epsabs.le.0 and
+                            (epsabs <= 0 and
-                             epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
+                             epsrel < max(50*rel.mach.acc.,0.5d-28))
-                             or limit.lt.1 or leniw.lt.limit*4.
+                             or limit < 1 or leniw < limit*4.
                             result, abserr, neval, last are set to
                             zero. exept when limit or leniw is
                             invalid, iwork(1), work(limit*2+1) and
                             work(limit*3+1) are set to zero, work(1)
                             is set to a and work(limit+1) to b.
         dimensioning parameters
            limit - int
                    dimensioning parameter for iwork
                    limit determines the maximum number of subintervals
                    in the partition of the given integration interval
-                   (a,b), limit.ge.1.
+                   (a,b), limit >= 1.
-                   if limit.lt.1, the routine will end with ier = 6.
+                   if limit < 1, the routine will end with ier = 6.
            lenw  - int
                    dimensioning parameter for work
                    lenw must be at least limit*4.
-                   if lenw.lt.limit*4, the routine will end
+                   if lenw < limit*4, the routine will end
                    with ier = 6.
            last  - int
                    on return, last equals the number of subintervals
                    produced in the subdivision process, which
                    vector of dimension at least limit, the first
                    k elements of which contain pointers
                    to the error estimates over the subintervals,
                    such that work(limit*3+iwork(1)),... ,
                    work(limit*3+iwork(k)) form a decreasing
-                   sequence, with k = last if last.le.(limit/2+2), and
+                   sequence, with k = last if last <= (limit/2+2), and
                    k = limit+1-last otherwise
            work  - double precision
                    vector of dimension at least lenw
                    on return
     l1 = *limit + 1;
     l2 = *limit + l1;
     l3 = *limit + l2;
-    rdqagie(IS, bound, inf, epsabs, epsrel, limit, result, abserr, neval, ier,
+    rdqagie(f, ex, bound, inf, epsabs, epsrel, limit, result, abserr, neval, ier,
 	    &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
     return;
 } /* Rdqagi */
 static
-void rdqagie(IntStruct IS, double *bound, int *inf, double *
+void rdqagie(integr_fn f, void *ex, double *bound, int *inf, double *
 	     epsabs, double *epsrel, int *limit, double *result,
 	     double *abserr, int *neval, int *ier, double *alist__,
 	     double *blist, double *rlist, double *elist, int *
 	     iord, int *last)
+{
     /* System generated locals */
-    int i__1, i__2;
     double d__1, d__2;
     /* Local variables */
     double area, dres;
     int ksgn;
 ***purpose  the routine calculates an approximation result to a given
            integral   i = integral of f over (bound,+infinity)
            or i = integral of f over (-infinity,bound)
            or i = integral of f over (-infinity,+infinity),
            hopefully satisfying following claim for accuracy
-           abs(i-result).le.max(epsabs,epsrel*abs(i))
+           abs(i-result) <= max(epsabs,epsrel*abs(i))
 ***description
 integration over infinite intervals
 standard fortran subroutine
            epsabs - double precision
                     absolute accuracy requested
            epsrel - double precision
                     relative accuracy requested
-                    if  epsabs.le.0
+                    if  epsabs <= 0
-                    and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+                    and epsrel < max(50*rel.mach.acc.,0.5d-28),
                     the routine will end with ier = 6.
            limit  - int
                     gives an upper bound on the number of subintervals
-                    in the partition of (a,b), limit.ge.1
+                    in the partition of (a,b), limit >= 1
         on return
            result - double precision
                     approximation to the integral
            ier    - int
                     ier = 0 normal and reliable termination of the
                             routine. it is assumed that the requested
                             accuracy has been achieved.
-                  - ier.gt.0 abnormal termination of the routine. the
+                  - ier > 0 abnormal termination of the routine. the
                             estimates for result and error are less
                             reliable. it is assumed that the requested
                             accuracy has not been achieved.
            error messages
                     ier = 1 maximum number of subdivisions allowed
                         = 5 the integral is probably divergent, or
                             slowly convergent. it must be noted that
                             divergence can occur with any other value
                             of ier.
                         = 6 the input is invalid, because
-                            (epsabs.le.0 and
+                            (epsabs <= 0 and
-                             epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+                             epsrel < max(50*rel.mach.acc.,0.5d-28),
                             result, abserr, neval, last, rlist(1),
                             elist(1) and iord(1) are set to zero.
                             alist(1) and blist(1) are set to 0
                             and 1 respectively.
                     vector of dimension limit, the first k
                     elements of which are pointers to the
                     error estimates over the subintervals,
                     such that elist(iord(1)), ..., elist(iord(k))
                     form a decreasing sequence, with k = last
-                    if last.le.(limit/2+2), and k = limit+1-last
+                    if last <= (limit/2+2), and k = limit+1-last
                     otherwise
            last   - int
                     number of subintervals actually produced
                     in the subdivision process
     alist__[1] = 0.;
     blist[1] = 1.;
     rlist[1] = 0.;
     elist[1] = 0.;
     iord[1] = 0;
-/* Computing MAX */
-    d__1 = epmach * 50.;
     if (*epsabs <= 0. && (*epsrel < fmax2(epmach * 50., 5e-29))) *ier = 6;
     if (*ier == 6) return;
 /*           first approximation to the integral */
 /*           ----------------------------------- */
     boun = *bound;
     if (*inf == 2) {
 	boun = 0.;
+    }
-    rdqk15i(IS, &boun, inf, &c_b6, &c_b7, result, abserr, &defabs, &resabs);
+    rdqk15i(f, ex, &boun, inf, &c_b6, &c_b7, result, abserr, &defabs, &resabs);
 /*           test on accuracy */
     *last = 1;
     rlist[1] = *result;
     elist[1] = *abserr;
     iord[1] = 1;
     dres = fabs(*result);
-/* Computing MAX */
-    d__1 = *epsabs, d__2 = *epsrel * dres;
+    errbnd = fmax2(*epsabs, *epsrel * dres);
-    errbnd = fmax2(d__1,d__2);
     if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) *ier = 2;
     if (*limit == 1) *ier = 1;
     if (*ier != 0 || (*abserr <= errbnd && *abserr != resabs)
 	|| *abserr == 0.) goto L130;
+    }
 /*           main do-loop */
 /*           ------------ */
-    i__1 = *limit;
-    for (*last = 2; *last <= i__1; ++(*last)) {
+    for (*last = 2; *last <= *limit; ++(*last)) {
 /*           bisect the subinterval with nrmax-th largest error estimate. */
 	a1 = alist__[maxerr];
 	b1 = (alist__[maxerr] + blist[maxerr]) * .5;
 	a2 = b1;
 	b2 = blist[maxerr];
 	erlast = errmax;
-	rdqk15i(IS, &boun, inf, &a1, &b1, &area1, &error1, &resabs, &defab1);
+	rdqk15i(f, ex, &boun, inf, &a1, &b1, &area1, &error1, &resabs, &defab1);
-	rdqk15i(IS, &boun, inf, &a2, &b2, &area2, &error2, &resabs, &defab2);
+	rdqk15i(f, ex, &boun, inf, &a2, &b2, &area2, &error2, &resabs, &defab2);
-/*           improve previous approximations to integral */
+/*           improve previous approximations to integral
-/*           and error and test for accuracy. */
+	     and error and test for accuracy. */
 	area12 = area1 + area2;
 	erro12 = error1 + error2;
 	errsum = errsum + erro12 - errmax;
 	area = area + area12 - rlist[maxerr];
 	if (defab1 == error1 || defab2 == error2) {
 	    goto L15;
+	}
-	if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
+	if (fabs(rlist[maxerr] - area12) > fabs(area12) * 1e-5 ||
 		 erro12 < errmax * .99) {
 	    goto L10;
+	}
 	if (extrap) {
 	    ++iroff2;
 	    ++iroff3;
+	}
 L15:
 	rlist[maxerr] = area1;
 	rlist[*last] = area2;
-/* Computing MAX */
-	d__1 = *epsabs, d__2 = *epsrel * fabs(area);
+	errbnd = fmax2(*epsabs, *epsrel * fabs(area));
-	errbnd = fmax2(d__1,d__2);
 /*           test for roundoff error and eventually set error flag. */
-	if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
+	if (iroff1 + iroff2 >= 10 || iroff3 >= 20)
 	    *ier = 2;
-	}
-	if (iroff2 >= 5) {
+	if (iroff2 >= 5)
 	    ierro = 3;
-	}
-/*           set error flag in the case that the number of */
+/*           set error flag in the case that the number of
-/*           subintervals equals limit. */
+	     subintervals equals limit. */
-	if (*last == *limit) {
+	if (*last == *limit)
 	    *ier = 1;
-	}
-/*           set error flag in the case of bad integrand behaviour */
+/*           set error flag in the case of bad integrand behaviour
-/*           at some points of the integration range. */
+	     at some points of the integration range. */
-/* Computing MAX */
-	d__1 = fabs(a1), d__2 = fabs(b2);
+	if (fmax2(fabs(a1), fabs(b2)) <=
-	if (fmax2(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
+	    (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
+		{
 	    *ier = 4;
+	}
 /*           append the newly-created intervals to the list. */
 	rlist[maxerr] = area2;
 	rlist[*last] = area1;
 	elist[maxerr] = error2;
 	elist[*last] = error1;
-/*           call subroutine dqpsrt to maintain the descending ordering */
+/*           call subroutine dqpsrt to maintain the descending ordering
-/*           in the list of error estimates and select the subinterval */
+	     in the list of error estimates and select the subinterval
-/*           with nrmax-th largest error estimate (to be bisected next). */
+	     with nrmax-th largest error estimate (to be bisected next). */
 L30:
 	rdqpsrt(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
 	if (errsum <= errbnd) {
 	    goto L115;
+	}
-	if (*ier != 0) {
+	if (*ier != 0)	    goto L100;
-	    goto L100;
-	}
-	if (*last == 2) {
+	if (*last == 2)     goto L80;
-	    goto L80;
-	}
-	if (noext) {
-	    goto L90;
+	if (noext) 	    goto L90;
-	}
 	erlarg -= erlast;
-	if ((d__1 = b1 - a1, fabs(d__1)) > small) {
+	if (fabs(b1 - a1) > small) {
 	    erlarg += erro12;
+	}
 	if (extrap) {
 	    goto L40;
+	}
-/*           test whether the interval to be bisected next is the */
+/*           test whether the interval to be bisected next is the
-/*           smallest interval. */
+	     smallest interval. */
-	if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
+	if (fabs(blist[maxerr] - alist__[maxerr]) > small) {
 	    goto L90;
+	}
 	extrap = TRUE;
 	nrmax = 2;
 L40:
 	if (ierro == 3 || erlarg <= ertest) {
 	    goto L60;
+	}
-/*           the smallest interval has the largest error. */
+/*           the smallest interval has the largest error.
-/*           before bisecting decrease the sum of the errors over the */
+	     before bisecting decrease the sum of the errors over the
-/*           larger intervals (erlarg) and perform extrapolation. */
+	     larger intervals (erlarg) and perform extrapolation. */
 	id = nrmax;
 	jupbnd = *last;
 	if (*last > *limit / 2 + 2) {
 	    jupbnd = *limit + 3 - *last;
+	}
-	i__2 = jupbnd;
-	for (k = id; k <= i__2; ++k) {
+	for (k = id; k <= jupbnd; ++k) {
 	    maxerr = iord[nrmax];
 	    errmax = elist[maxerr];
-	    if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
+	    if (fabs(blist[maxerr] - alist__[maxerr]) > small) {
 		goto L90;
+	    }
 	    ++nrmax;
 /* L50: */
+	}
 /*           compute global integral sum. */
 L115:
     *result = 0.;
-    i__1 = *last;
-    for (k = 1; k <= i__1; ++k) {
+    for (k = 1; k <= *last; ++k)
 	*result += rlist[k];
-/* L120: */
-    }
     *abserr = errsum;
 L130:
     *neval = *last * 30 - 15;
     if (*inf == 2) {
 	*neval <<= 1;
 	--(*ier);
+    }
     return;
 } /* rdqagie_ */
-static
-void Rdqags(IntStruct IS, double *a, double *b, double *epsabs,
+void Rdqags(integr_fn f, void *ex, double *a, double *b,
-	    double *epsrel, double *result, double *abserr, int *
+	    double *epsabs, double *epsrel,
-	    neval, int *ier, int *limit, int *lenw, int *last,
+	    double *result, double *abserr, int *neval, int *ier,
-	    int *iwork, double *work)
+	    int *limit, int *lenw, int *last, int *iwork, double *work)
+{
     int l1, l2, l3;
 /*
 ***begin prologue  dqags
 ***date written   800101   (yymmdd)
 ***revision date  830518   (yymmdd)
 ***category no.  h2a1a1
 ***keywords  automatic integrator, general-purpose,
 ***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
           de doncker,elise,appl. math. & prog. div. - k.u.leuven
 ***purpose  the routine calculates an approximation result to a given
            definite integral  i = integral of f over (a,b),
            hopefully satisfying following claim for accuracy
-           abs(i-result).le.max(epsabs,epsrel*abs(i)).
+           abs(i-result) <= max(epsabs,epsrel*abs(i)).
 ***description
        computation of a definite integral
        standard fortran subroutine
        double precision version
            epsabs - double precision
                     absolute accuracy requested
            epsrel - double precision
                     relative accuracy requested
-                    if  epsabs.le.0
+                    if  epsabs <= 0
-                    and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+                    and epsrel < max(50*rel.mach.acc.,0.5d-28),
                     the routine will end with ier = 6.
         on return
            result - double precision
                     approximation to the integral
            ier    - int
                     ier = 0 normal and reliable termination of the
                             routine. it is assumed that the requested
                             accuracy has been achieved.
-                    ier.gt.0 abnormal termination of the routine
+                    ier > 0 abnormal termination of the routine
                             the estimates for integral and error are
                             less reliable. it is assumed that the
                             requested accuracy has not been achieved.
            error messages
                     ier = 1 maximum number of subdivisions allowed
                         = 5 the integral is probably divergent, or
                             slowly convergent. it must be noted that
                             divergence can occur with any other value
                             of ier.
                         = 6 the input is invalid, because
-                            (epsabs.le.0 and
+                            (epsabs <= 0 and
-                             epsrel.lt.max(50*rel.mach.acc.,0.5d-28)
+                             epsrel < max(50*rel.mach.acc.,0.5d-28)
-                            or limit.lt.1 or lenw.lt.limit*4.
+                            or limit < 1 or lenw < limit*4.
                             result, abserr, neval, last are set to
                             zero.except when limit or lenw is invalid,
                             iwork(1), work(limit*2+1) and
                             work(limit*3+1) are set to zero, work(1)
                             is set to a and work(limit+1) to b.
         dimensioning parameters
            limit - int
                    dimensioning parameter for iwork
                    limit determines the maximum number of subintervals
                    in the partition of the given integration interval
-                   (a,b), limit.ge.1.
+                   (a,b), limit >= 1.
-                   if limit.lt.1, the routine will end with ier = 6.
+                   if limit < 1, the routine will end with ier = 6.
            lenw  - int
                    dimensioning parameter for work
                    lenw must be at least limit*4.
-                   if lenw.lt.limit*4, the routine will end
+                   if lenw < limit*4, the routine will end
                    with ier = 6.
            last  - int
                    on return, last equals the number of subintervals
                    produced in the subdivision process, detemines the
                    vector of dimension at least limit, the first k
                    elements of which contain pointers
                    to the error estimates over the subintervals
                    such that work(limit*3+iwork(1)),... ,
                    work(limit*3+iwork(k)) form a decreasing
-                   sequence, with k = last if last.le.(limit/2+2),
+                   sequence, with k = last if last <= (limit/2+2),
                    and k = limit+1-last otherwise
            work  - double precision
                    vector of dimension at least lenw
                    on return
     l1 = *limit + 1;
     l2 = *limit + l1;
     l3 = *limit + l2;
-    rdqagse(IS, a, b, epsabs, epsrel, limit, result, abserr, neval, ier,
+    rdqagse(f, ex, a, b, epsabs, epsrel, limit, result, abserr, neval, ier,
 	    &work[1], &work[l1], &work[l2], &work[l3], &iwork[1], last);
     return;
 } /* rdqags_ */
 static
-void rdqagse(IntStruct IS, double *a, double *b, double *
+void rdqagse(integr_fn f, void *ex, double *a, double *b, double *
 	     epsabs, double *epsrel, int *limit, double *result,
 	     double *abserr, int *neval, int *ier, double *alist__,
 	     double *blist, double *rlist, double *elist, int *
 	     iord, int *last)
+{
-    /* System generated locals */
-    int i__1, i__2;
-    double d__1, d__2;
     /* Local variables */
-    double area, dres;
+    Rboolean noext, extrap;
-    int ksgn, nres;
+    int k,ksgn, nres;
-    double area1, area2, area12;
-    int k;
-    double small = 0.0, erro12;
     int ierro;
-    double a1, a2, b1, b2, defab1, defab2, oflow;
     int ktmin, nrmax;
-    double uflow;
-    Rboolean noext;
     int iroff1, iroff2, iroff3;
-    double res3la[3], error1, error2;
     int id;
-    double rlist2[52];
     int numrl2;
-    double defabs, epmach, erlarg = 0.0, abseps, correc = 0.0, errbnd, resabs;
     int jupbnd;
-    double erlast, errmax;
     int maxerr;
-    double reseps;
+    double res3la[3];
-    Rboolean extrap;
+    double rlist2[52];
-    double ertest = 0.0, errsum;
+    double abseps, area, area1, area2, area12, dres, epmach;
+    double a1, a2, b1, b2, defabs, defab1, defab2, oflow, uflow, resabs, reseps;
+    double error1, error2, erro12, errbnd, erlast, errmax, errsum;
+    double correc = 0.0, erlarg = 0.0, ertest = 0.0, small = 0.0;
 /*
 ***begin prologue  dqagse
 ***date written   800101   (yymmdd)
 ***revision date  830518   (yymmdd)
 ***category no.  h2a1a1
 ***keywords  automatic integrator, general-purpose,
 ***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
           de doncker,elise,appl. math. & progr. div. - k.u.leuven
 ***purpose  the routine calculates an approximation result to a given
            definite integral i = integral of f over (a,b),
            hopefully satisfying following claim for accuracy
-           abs(i-result).le.max(epsabs,epsrel*abs(i)).
+           abs(i-result) <= max(epsabs,epsrel*abs(i)).
 ***description
        computation of a definite integral
        standard fortran subroutine
        double precision version
            epsabs - double precision
                     absolute accuracy requested
            epsrel - double precision
                     relative accuracy requested
-                    if  epsabs.le.0
+                    if  epsabs <= 0
-                    and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+                    and epsrel < max(50*rel.mach.acc.,0.5d-28),
                     the routine will end with ier = 6.
            limit  - int
                     gives an upperbound on the number of subintervals
                     in the partition of (a,b)
            ier    - int
                     ier = 0 normal and reliable termination of the
                             routine. it is assumed that the requested
                             accuracy has been achieved.
-                    ier.gt.0 abnormal termination of the routine
+                    ier > 0 abnormal termination of the routine
                             the estimates for integral and error are
                             less reliable. it is assumed that the
                             requested accuracy has not been achieved.
            error messages
                         = 1 maximum number of subdivisions allowed
                         = 5 the integral is probably divergent, or
                             slowly convergent. it must be noted that
                             divergence can occur with any other value
                             of ier.
                         = 6 the input is invalid, because
-                            epsabs.le.0 and
+                            epsabs <= 0 and
-                            epsrel.lt.max(50*rel.mach.acc.,0.5d-28).
+                            epsrel < max(50*rel.mach.acc.,0.5d-28).
                             result, abserr, neval, last, rlist(1),
                             iord(1) and elist(1) are set to zero.
                             alist(1) and blist(1) are set to a and b
                             respectively.
                     vector of dimension at least limit, the first k
                     elements of which are pointers to the
                     error estimates over the subintervals,
                     such that elist(iord(1)), ..., elist(iord(k))
                     form a decreasing sequence, with k = last
-                    if last.le.(limit/2+2), and k = limit+1-last
+                    if last <= (limit/2+2), and k = limit+1-last
                     otherwise
            last   - int
                     number of subintervals actually produced in the
                     subdivision process
     *abserr = 0.;
     alist__[1] = *a;
     blist[1] = *b;
     rlist[1] = 0.;
     elist[1] = 0.;
-/* Computing MAX */
-    d__1 = epmach * 50.;
-    if (*epsabs <= 0. && (*epsrel < fmax2(d__1,5e-29))) {
+    if (*epsabs <= 0. && *epsrel < fmax2(epmach * 50., 5e-29)) {
 	*ier = 6;
+	return;
+    }
-    if (*ier == 6) return;
 /*           first approximation to the integral */
 /*           ----------------------------------- */
     uflow = DBL_MIN;
     oflow = DBL_MAX;
     ierro = 0;
-    rdqk21(IS, a, b, result, abserr, &defabs, &resabs);
+    rdqk21(f, ex, a, b, result, abserr, &defabs, &resabs);
 /*           test on accuracy. */
     dres = fabs(*result);
-/* Computing MAX */
-    d__1 = *epsabs, d__2 = *epsrel * dres;
+    errbnd = fmax2(*epsabs, *epsrel * dres);
-    errbnd = fmax2(d__1,d__2);
     *last = 1;
     rlist[1] = *result;
     elist[1] = *abserr;
     iord[1] = 1;
-    if (*abserr <= epmach * 100. * defabs && *abserr > errbnd) {
+    if (*abserr <= epmach * 100. * defabs && *abserr > errbnd)
 	*ier = 2;
-    }
-    if (*limit == 1) {
+    if (*limit == 1)
 	*ier = 1;
-    }
     if (*ier != 0 || (*abserr <= errbnd && *abserr != resabs)
 	|| *abserr == 0.) goto L140;
 /*           initialization */
 /*           -------------- */
+    }
 /*           main do-loop */
 /*           ------------ */
-    i__1 = *limit;
-    for (*last = 2; *last <= i__1; ++(*last)) {
+    for (*last = 2; *last <= *limit; ++(*last)) {
-/*           bisect the subinterval with the nrmax-th largest error */
+/*           bisect the subinterval with the nrmax-th largest error estimate. */
-/*           estimate. */
 	a1 = alist__[maxerr];
 	b1 = (alist__[maxerr] + blist[maxerr]) * .5;
 	a2 = b1;
 	b2 = blist[maxerr];
 	erlast = errmax;
-	rdqk21(IS, &a1, &b1, &area1, &error1, &resabs, &defab1);
+	rdqk21(f, ex, &a1, &b1, &area1, &error1, &resabs, &defab1);
-	rdqk21(IS, &a2, &b2, &area2, &error2, &resabs, &defab2);
+	rdqk21(f, ex, &a2, &b2, &area2, &error2, &resabs, &defab2);
-/*           improve previous approximations to integral */
+/*           improve previous approximations to integral
-/*           and error and test for accuracy. */
+	     and error and test for accuracy. */
 	area12 = area1 + area2;
 	erro12 = error1 + error2;
 	errsum = errsum + erro12 - errmax;
 	area = area + area12 - rlist[maxerr];
 	if (defab1 == error1 || defab2 == error2) {
 	    goto L15;
+	}
-	if ((d__1 = rlist[maxerr] - area12, fabs(d__1)) > fabs(area12) * 1e-5 ||
+	if (fabs(rlist[maxerr] - area12) > fabs(area12) * 1e-5 ||
 		 erro12 < errmax * .99) {
 	    goto L10;
+	}
 	if (extrap) {
 	    ++iroff2;
 	    ++iroff3;
+	}
 L15:
 	rlist[maxerr] = area1;
 	rlist[*last] = area2;
-/* Computing MAX */
-	d__1 = *epsabs, d__2 = *epsrel * fabs(area);
+	errbnd = fmax2(*epsabs, *epsrel * fabs(area));
-	errbnd = fmax2(d__1,d__2);
 /*           test for roundoff error and eventually set error flag. */
-	if (iroff1 + iroff2 >= 10 || iroff3 >= 20) {
+	if (iroff1 + iroff2 >= 10 || iroff3 >= 20)
 	    *ier = 2;
-	}
-	if (iroff2 >= 5) {
+	if (iroff2 >= 5)
 	    ierro = 3;
-	}
-/*           set error flag in the case that the number of subintervals */
-/*           equals limit. */
+/* set error flag in the case that the number of subintervals equals limit. */
-	if (*last == *limit) {
+	if (*last == *limit)
 	    *ier = 1;
-	}
-/*           set error flag in the case of bad integrand behaviour */
+/*           set error flag in the case of bad integrand behaviour
-/*           at a point of the integration range. */
+	     at a point of the integration range. */
-/* Computing MAX */
-	d__1 = fabs(a1), d__2 = fabs(b2);
+	if (fmax2(fabs(a1), fabs(b2)) <=
-	if (fmax2(d__1,d__2) <= (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3))
+	    (epmach * 100. + 1.) * (fabs(a2) + uflow * 1e3)) {
-		{
 	    *ier = 4;
+	}
 /*           append the newly-created intervals to the list. */
 	if (error2 > error1) {
-	    goto L20;
+	    alist__[maxerr] = a2;
-	}
-	alist__[*last] = a2;
+	    alist__[*last] = a1;
+	    blist[*last] = b1;
-	blist[maxerr] = b1;
+	    rlist[maxerr] = area2;
-	blist[*last] = b2;
+	    rlist[*last] = area1;
-	elist[maxerr] = error1;
+	    elist[maxerr] = error2;
-	elist[*last] = error2;
+	    elist[*last] = error1;
-	goto L30;
+	} else {
-L20:
-	alist__[maxerr] = a2;
+	    alist__[*last] = a2;
-	alist__[*last] = a1;
+	    blist[maxerr] = b1;
-	blist[*last] = b1;
+	    blist[*last] = b2;
-	rlist[maxerr] = area2;
+	    elist[maxerr] = error1;
-	rlist[*last] = area1;
+	    elist[*last] = error2;
+	}
+/*           call subroutine dqpsrt to maintain the descending ordering
-	elist[maxerr] = error2;
+	     in the list of error estimates and select the subinterval
-	elist[*last] = error1;
+	     with nrmax-th largest error estimate (to be bisected next). */
-/*           call subroutine dqpsrt to maintain the descending ordering */
-/*           in the list of error estimates and select the subinterval */
-/*           with nrmax-th largest error estimate (to be bisected next). */
-L30:
+/*L30:*/
 	rdqpsrt(limit, last, &maxerr, &errmax, &elist[1], &iord[1], &nrmax);
-/* ***jump out of do-loop */
-	if (errsum <= errbnd) {
-	    goto L115;
-	}
-/* ***jump out of do-loop */
+	if (errsum <= errbnd)   goto L115;/* ***jump out of do-loop */
-	if (*ier != 0) {
+	if (*ier != 0)		goto L100;/* ***jump out of do-loop */
-	    goto L100;
-	}
-	if (*last == 2) {
+	if (*last == 2)		goto L80;
-	    goto L80;
-	}
-	if (noext) {
+	if (noext)		goto L90;
-	    goto L90;
-	}
 	erlarg -= erlast;
-	if ((d__1 = b1 - a1, fabs(d__1)) > small) {
+	if (fabs(b1 - a1) > small) {
 	    erlarg += erro12;
+	}
 	if (extrap) {
 	    goto L40;
+	}
-/*           test whether the interval to be bisected next is the */
+/*           test whether the interval to be bisected next is the
-/*           smallest interval. */
+	     smallest interval. */
-	if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
+	if (fabs(blist[maxerr] - alist__[maxerr]) > small) {
 	    goto L90;
+	}
 	extrap = TRUE;
 	nrmax = 2;
 L40:
 	if (ierro == 3 || erlarg <= ertest) {
 	    goto L60;
+	}
-/*           the smallest interval has the largest error. */
+/*           the smallest interval has the largest error.
-/*           before bisecting decrease the sum of the errors over the */
+	     before bisecting decrease the sum of the errors over the
-/*           larger intervals (erlarg) and perform extrapolation. */
+	     larger intervals (erlarg) and perform extrapolation. */
 	id = nrmax;
 	jupbnd = *last;
 	if (*last > *limit / 2 + 2) {
 	    jupbnd = *limit + 3 - *last;
+	}
-	i__2 = jupbnd;
-	for (k = id; k <= i__2; ++k) {
+	for (k = id; k <= jupbnd; ++k) {
 	    maxerr = iord[nrmax];
 	    errmax = elist[maxerr];
-/* ***jump out of do-loop */
-	    if ((d__1 = blist[maxerr] - alist__[maxerr], fabs(d__1)) > small) {
+	    if (fabs(blist[maxerr] - alist__[maxerr]) > small) {
-		goto L90;
+		goto L90;/* ***jump out of do-loop */
+	    }
 	    ++nrmax;
 /* L50: */
+	}
+	}
 	ktmin = 0;
 	*abserr = abseps;
 	*result = reseps;
 	correc = erlarg;
-/* Computing MAX */
-	d__1 = *epsabs, d__2 = *epsrel * fabs(reseps);
+	ertest = fmax2(*epsabs, *epsrel * fabs(reseps));
-	ertest = fmax2(d__1,d__2);
-/* ***jump out of do-loop */
 	if (*abserr <= ertest) {
-	    goto L100;
+	    goto L100;/* ***jump out of do-loop */
+	}
 /*           prepare bisection of the smallest interval. */
 L70:
 	extrap = FALSE;
 	small *= .5;
 	erlarg = errsum;
 	goto L90;
 L80:
-	small = (d__1 = *b - *a, fabs(d__1)) * .375;
+	small = fabs(*b - *a) * .375;
 	erlarg = errsum;
 	ertest = errbnd;
 	rlist2[1] = area;
 L90:
+	;
+    }
-/*           set final result and error estimate. */
-/*           ------------------------------------ */
-L100:
+L100:/*		set final result and error estimate. */
+/*		------------------------------------ */
-    if (*abserr == oflow) {
+    if (*abserr == oflow) 	goto L115;
-	goto L115;
-    }
-    if (*ier + ierro == 0) {
+    if (*ier + ierro == 0) 	goto L110;
-	goto L110;
-    }
-    if (ierro == 3) {
+    if (ierro == 3)
 	*abserr += correc;
-    }
-    if (*ier == 0) {
+    if (*ier == 0)
 	*ier = 3;
-    }
-    if (*result != 0. && area != 0.) {
+    if (*result != 0. && area != 0.) goto L105;
-	goto L105;
-    }
-    if (*abserr > errsum) {
+    if (*abserr > errsum) 	goto L115;
-	goto L115;
-    }
-    if (area == 0.) {
+    if (area == 0.) 		goto L130;
-	goto L130;
-    }
     goto L110;
 L105:
     if (*abserr / fabs(*result) > errsum / fabs(area)) {
 	goto L115;
+    }
-/*           test on divergence. */
+L110:/*		test on divergence. */
-L110:
-/* Computing MAX */
-    d__1 = fabs(*result), d__2 = fabs(area);
-    if (ksgn == -1 && fmax2(d__1,d__2) <= defabs * .01) {
+    if (ksgn == -1 && fmax2(fabs(*result), fabs(area)) <= defabs * .01) {
 	goto L130;
+    }
     if (.01 > *result / area || *result / area > 100. || errsum > fabs(area)) {
 	*ier = 5;
+    }
     goto L130;
-/*           compute global integral sum. */
+L115:/*		compute global integral sum. */
-L115:
     *result = 0.;
-    i__1 = *last;
-    for (k = 1; k <= i__1; ++k) {
+    for (k = 1; k <= *last; ++k)
 	*result += rlist[k];
-/* L120: */
-    }
     *abserr = errsum;
 L130:
     if (*ier > 2)
 L140:
     *neval = *last * 42 - 21;
     return;
 } /* rdqagse_ */
-static void rdqk15i(IntStruct IS, double *boun, int *inf, double *a,
+static void rdqk15i(integr_fn f, void *ex,
-		    double *b, double *result, double *abserr, double *
+		    double *boun, int *inf, double *a, double *b,
-		    resabs, double *resasc)
+		    double *result,
+		    double *abserr, double *resabs, double *resasc)
+{
     /* Initialized data */
+    static double wg[8] = {
-    static double wg[8] = { 0.,.129484966168869693270611432679082,0.,
+., .129484966168869693270611432679082,
-	    .27970539148927666790146777142378,0.,
+., .27970539148927666790146777142378,
-	    .381830050505118944950369775488975,0.,
+., .381830050505118944950369775488975,
-	    .417959183673469387755102040816327 };
+., .417959183673469387755102040816327 };
+    static double xgk[8] = {
-    static double xgk[8] = { .991455371120812639206854697526329,
+	    .991455371120812639206854697526329,
 	    .949107912342758524526189684047851,
 	    .864864423359769072789712788640926,
 	    .741531185599394439863864773280788,
 	    .58608723546769113029414483825873,
 	    .405845151377397166906606412076961,
-	    .207784955007898467600689403773245,0. };
+	    .207784955007898467600689403773245, 0. };
+    static double wgk[8] = {
-    static double wgk[8] = { .02293532201052922496373200805897,
+	    .02293532201052922496373200805897,
 	    .063092092629978553290700663189204,
 	    .104790010322250183839876322541518,
 	    .140653259715525918745189590510238,
 	    .16900472663926790282658342659855,
 	    .190350578064785409913256402421014,
 	    .204432940075298892414161999234649,
 	    .209482141084727828012999174891714 };
-    /* System generated locals */
-    double d__1, d__2, d__3;
     /* Local variables */
     double absc, dinf, resg, resk, fsum, absc1, absc2, fval1, fval2;
     int j;
     double hlgth, centr, reskh, uflow;
     double tabsc1, tabsc2, fc, epmach;
     double fv1[7], fv2[7], vec[15], vec2[15];
 /*
 ***begin prologue  dqk15i
 ***date written   800101   (yymmdd)
 ***revision date  830518   (yymmdd)
 ***category no.  h2a3a2,h2a4a2
 ***keywords  15-point transformed gauss-kronrod rules
                       approximation to the integral of
                       abs((transformed integrand)-i/(b-a)) over (a,b)
 ***references  (none)
 ***end prologue  dqk15i
           the abscissae and weights are supplied for the interval
           (-1,1).  because of symmetry only the positive abscissae and
           their corresponding weights are given.
           machine dependent constants
           ---------------------------
           epmach is the largest relative spacing.
-          uflow is the smallest positive magnitude. */
+          uflow is the smallest positive magnitude.
+*/
 /* ***first executable statement  dqk15i */
     epmach = DBL_EPSILON;
     uflow = DBL_MIN;
     dinf = (double) imin2(1, *inf);
 	    vec2[(j << 1) - 1] = -tabsc1;
 	    vec2[j * 2] = -tabsc2;
+	}
 /* L5: */
+    }
-    Rintfn(vec, 15, IS);
+    f(vec, 15, ex); /* -> new vec[] overwriting old vec[] */
-    if (*inf == 2) Rintfn(vec2, 15, IS);
+    if (*inf == 2) f(vec2, 15, ex);
     fval1 = vec[0];
     if (*inf == 2) fval1 += vec2[0];
     fc = fval1 / centr / centr;
-/*           compute the 15-point kronrod approximation to */
+/*           compute the 15-point kronrod approximation to
-/*           the integral, and estimate the error. */
+	     the integral, and estimate the error. */
     resg = wg[7] * fc;
     resk = wgk[7] * fc;
     *resabs = fabs(resk);
     for (j = 1; j <= 7; ++j) {
 	resk += wgk[j - 1] * fsum;
 	*resabs += wgk[j - 1] * (fabs(fval1) + fabs(fval2));
 /* L10: */
+    }
     reskh = resk * .5;
-    *resasc = wgk[7] * (d__1 = fc - reskh, fabs(d__1));
+    *resasc = wgk[7] * fabs(fc - reskh);
     for (j = 1; j <= 7; ++j) {
-	*resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
+	*resasc += wgk[j - 1] * (fabs(fv1[j - 1] - reskh) +
-		d__2 = fv2[j - 1] - reskh, fabs(d__2)));
+				 fabs(fv2[j - 1] - reskh));
 /* L20: */
+    }
     *result = resk * hlgth;
     *resasc *= hlgth;
     *resabs *= hlgth;
-    *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
+    *abserr = fabs((resk - resg) * hlgth);
     if (*resasc != 0. && *abserr != 0.) {
-/* Computing MIN */
-	d__3 = *abserr * 200. / *resasc;
-	d__1 = 1., d__2 = pow(d__3, 1.5);
-	*abserr = *resasc * fmin2(d__1,d__2);
+	*abserr = *resasc * fmin2(1., pow(*abserr * 200. / *resasc, 1.5));
+    }
     if (*resabs > uflow / (epmach * 50.)) {
-/* Computing MAX */
-	d__1 = epmach * 50. * *resabs;
-	*abserr = fmax2(d__1,*abserr);
+	*abserr = fmax2(epmach * 50. * *resabs, *abserr);
+    }
     return;
 } /* rdqk15i_ */
 static void rdqelg(int *n, double *epstab, double *
 		   result, double *abserr, double *res3la, int *nres)
+{
-    /* System generated locals */
-    int i__1;
-    double d__1, d__2, d__3;
     /* Local variables */
-    int indx;
-    double e1abs;
-    int i__;
-    double e0, e1, e2, e3, error, oflow;
+    int i__, indx, ib, ib2, ie, k1, k2, k3, num, newelm, limexp;
-    int k1;
-    int k2, k3;
-    double delta1, delta2, delta3;
+    double delta1, delta2, delta3, e0, e1, e1abs, e2, e3, epmach, epsinf;
-    int ib, ie;
-    double epmach, ss, epsinf;
+    double oflow, ss, res;
-    int newelm, ib2, limexp;
-    double res;
-    int num;
-    double err1, err2, err3, tol1, tol2, tol3;
+    double errA, err1, err2, err3, tol1, tol2, tol3;
 /* ***begin prologue  dqelg
 ***refer to  dqagie,dqagoe,dqagpe,dqagse
 ***revision date  830518   (yymmdd)
 ***keywords  epsilon algorithm, convergence acceleration,
           e1       element in the epsilon table is based
           e2
           e3                 e0
                        e3    e1    new
                              e2
-          newelm - number of elements to be computed in the new
+          newelm - number of elements to be computed in the new diagonal
-                   diagonal
-          error  - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
+          errA   - errA = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
-          result - the element in the new diagonal with least value
+          result - the element in the new diagonal with least value of errA
-                   of error
           machine dependent constants
           ---------------------------
           epmach is the largest relative spacing.
     epstab[*n + 2] = epstab[*n];
     newelm = (*n - 1) / 2;
     epstab[*n] = oflow;
     num = *n;
     k1 = *n;
-    i__1 = newelm;
-    for (i__ = 1; i__ <= i__1; ++i__) {
+    for (i__ = 1; i__ <= newelm; ++i__) {
 	k2 = k1 - 1;
 	k3 = k1 - 2;
 	res = epstab[k1 + 2];
 	e0 = epstab[k3];
 	e1 = epstab[k2];
 	e2 = res;
 	e1abs = fabs(e1);
 	delta2 = e2 - e1;
 	err2 = fabs(delta2);
-/* Computing MAX */
-	d__1 = fabs(e2);
-	tol2 = fmax2(d__1,e1abs) * epmach;
+	tol2 = fmax2(fabs(e2), e1abs) * epmach;
 	delta3 = e1 - e0;
 	err3 = fabs(delta3);
-/* Computing MAX */
+	tol3 = fmax2(e1abs, fabs(e0)) * epmach;
-	d__1 = e1abs, d__2 = fabs(e0);
+	if (err2 <= tol2 && err3 <= tol3) {
+	    /*           if e0, e1 and e2 are equal to within machine
-	tol3 = fmax2(d__1,d__2) * epmach;
+			 accuracy, convergence is assumed. */
-	if (err2 > tol2 || err3 > tol3) {
+	    *result = res;/*		result = e2 */
+	    *abserr = err2 + err3;/*	abserr = fabs(e1-e0)+fabs(e2-e1) */
-	    goto L10;
+	    goto L100;	/* ***jump out of do-loop */
+	}
-/*           if e0, e1 and e2 are equal to within machine */
-/*           accuracy, convergence is assumed. */
-/*           result = e2 */
-/*           abserr = fabs(e1-e0)+fabs(e2-e1) */
-	*result = res;
-	*abserr = err2 + err3;
-/* ***jump out of do-loop */
-	goto L100;
-L10:
 	e3 = epstab[k1];
 	epstab[k1] = e1;
 	delta1 = e1 - e3;
 	err1 = fabs(delta1);
-/* Computing MAX */
-	d__1 = e1abs, d__2 = fabs(e3);
-	tol1 = fmax2(d__1,d__2) * epmach;
+	tol1 = fmax2(e1abs, fabs(e3)) * epmach;
-/*           if two elements are very close to each other, omit */
+/*           if two elements are very close to each other, omit
-/*           a part of the table by adjusting the value of n */
+	     a part of the table by adjusting the value of n */
-	if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3) {
+	if (err1 > tol1 && err2 > tol2 && err3 > tol3) {
-	    goto L20;
-	}
-	ss = 1. / delta1 + 1. / delta2 - 1. / delta3;
+	    ss = 1. / delta1 + 1. / delta2 - 1. / delta3;
-	epsinf = (d__1 = ss * e1, fabs(d__1));
+	    epsinf = fabs(ss * e1);
-/*           test to detect irregular behaviour in the table, and */
+/*           test to detect irregular behaviour in the table, and
-/*           eventually omit a part of the table adjusting the value */
+	     eventually omit a part of the table adjusting the value of n. */
-/*           of n. */
-	if (epsinf > 1e-4) {
+	    if (epsinf > 1e-4) {
-	    goto L30;
+		goto L30;
+	    }
+	}
-L20:
 	*n = i__ + i__ - 1;
-/* ***jump out of do-loop */
+	goto L50;/* ***jump out of do-loop */
-	goto L50;
-/*           compute a new element and eventually adjust */
-/*           the value of result. */
+L30:/* compute a new element and eventually adjust the value of result. */
-L30:
 	res = e1 + 1. / ss;
 	epstab[k1] = res;
 	k1 += -2;
-	error = err2 + (d__1 = res - e2, fabs(d__1)) + err3;
+	errA = err2 + fabs(res - e2) + err3;
-	if (error > *abserr) {
+	if (errA <= *abserr) {
+	    *abserr = errA;
-	    goto L40;
+	    *result = res;
+	}
-	*abserr = error;
-	*result = res;
-L40:
-	;
+    }
 /*           shift the table. */
 L50:
     if (*n == limexp) {
 	*n = (limexp / 2 << 1) - 1;
+    }
-    ib = 1;
-    if (num / 2 << 1 == num) {
+    if (num / 2 << 1 == num) ib = 2; else ib = 1;
-	ib = 2;
-    }
     ie = newelm + 1;
-    i__1 = ie;
-    for (i__ = 1; i__ <= i__1; ++i__) {
+    for (i__ = 1; i__ <= ie; ++i__) {
 	ib2 = ib + 2;
 	epstab[ib] = epstab[ib2];
 	ib = ib2;
-/* L60: */
+    }
-    if (num == *n) {
+    if (num != *n) {
-	goto L80;
-    }
-    indx = num - *n + 1;
+	indx = num - *n + 1;
-    i__1 = *n;
-    for (i__ = 1; i__ <= i__1; ++i__) {
+	for (i__ = 1; i__ <= *n; ++i__) {
-	epstab[i__] = epstab[indx];
+	    epstab[i__] = epstab[indx];
-	++indx;
+	    ++indx;
-/* L70: */
+	}
+    }
-L80:
+    /*L80:*/
     if (*nres >= 4) {
-	goto L90;
+	/* L90: */
+	*abserr = fabs(*result - res3la[3]) +
+	          fabs(*result - res3la[2]) +
+	          fabs(*result - res3la[1]);
+	res3la[1] = res3la[2];
+	res3la[2] = res3la[3];
+	res3la[3] = *result;
+    } else {
+	res3la[*nres] = *result;
+	*abserr = oflow;
+    }
-    res3la[*nres] = *result;
-    *abserr = oflow;
-    goto L100;
-/*           compute error estimate */
-L90:
-    *abserr = (d__1 = *result - res3la[3], fabs(d__1)) + (d__2 = *result -
-	    res3la[2], fabs(d__2)) + (d__3 = *result - res3la[1], fabs(d__3));
-    res3la[1] = res3la[2];
-    res3la[2] = res3la[3];
-    res3la[3] = *result;
-L100:
-/* Computing MAX */
+L100:/* compute error estimate */
-    d__1 = *abserr, d__2 = epmach * 5. * fabs(*result);
+    *abserr = fmax2(*abserr, epmach * 5. * fabs(*result));
-    *abserr = fmax2(d__1,d__2);
     return;
 } /* rdqelg_ */
-static void  rdqk21(IntStruct IS, double *a, double *b, double *result,
+static void  rdqk21(integr_fn f, void *ex, double *a, double *b, double *result,
 		    double *abserr, double *resabs, double *resasc)
+{
     /* Initialized data */
     static double wg[5] = { .066671344308688137593568809893332,
 	    .134709217311473325928054001771707,
 	    .142775938577060080797094273138717,
 	    .147739104901338491374841515972068,
 	    .149445554002916905664936468389821 };
-    /* System generated locals */
-    double d__1, d__2, d__3;
     /* Local variables */
+    double fv1[10], fv2[10], vec[21];
     double absc, resg, resk, fsum, fval1, fval2;
-    int jtwm1, j;
     double hlgth, centr, reskh, uflow;
     double fc, epmach, dhlgth;
-    double fv1[10], fv2[10], vec[21];
-    int jtw;
+    int j, jtw, jtwm1;
 /* ***begin prologue  dqk21
 ***date written   800101   (yymmdd)
 ***revision date  830518   (yymmdd)
 ***category no.  h2a1a2
     centr = (*a + *b) * .5;
     hlgth = (*b - *a) * .5;
     dhlgth = fabs(hlgth);
-/*           compute the 21-point kronrod approximation to */
+/*           compute the 21-point kronrod approximation to
-/*           the integral, and estimate the absolute error. */
+	     the integral, and estimate the absolute error. */
     resg = 0.;
     vec[0] = centr;
     for (j = 1; j <= 5; ++j) {
 	jtw = j << 1;
+    }
     for (j = 1; j <= 5; ++j) {
 	jtwm1 = (j << 1) - 1;
 	absc = hlgth * xgk[jtwm1 - 1];
 	vec[(j << 1) + 9] = centr - absc;
-/* L6: */
 	vec[(j << 1) + 10] = centr + absc;
+    }
-    Rintfn(vec, 21, IS);
+    f(vec, 21, ex);
     fc = vec[0];
     resk = wgk[10] * fc;
     *resabs = fabs(resk);
     for (j = 1; j <= 5; ++j) {
 	jtw = j << 1;
 	resk += wgk[jtwm1 - 1] * fsum;
 	*resabs += wgk[jtwm1 - 1] * (fabs(fval1) + fabs(fval2));
 /* L15: */
+    }
     reskh = resk * .5;
-    *resasc = wgk[10] * (d__1 = fc - reskh, fabs(d__1));
+    *resasc = wgk[10] * fabs(fc - reskh);
     for (j = 1; j <= 10; ++j) {
-	*resasc += wgk[j - 1] * ((d__1 = fv1[j - 1] - reskh, fabs(d__1)) + (
+	*resasc += wgk[j - 1] * (fabs(fv1[j - 1] - reskh) +
-		d__2 = fv2[j - 1] - reskh, fabs(d__2)));
+				 fabs(fv2[j - 1] - reskh));
 /* L20: */
+    }
     *result = resk * hlgth;
     *resabs *= dhlgth;
     *resasc *= dhlgth;
-    *abserr = (d__1 = (resk - resg) * hlgth, fabs(d__1));
+    *abserr = fabs((resk - resg) * hlgth);
     if (*resasc != 0. && *abserr != 0.) {
-/* Computing MIN */
-	d__3 = *abserr * 200. / *resasc;
-	d__1 = 1., d__2 = pow(d__3, 1.5);
-	*abserr = *resasc * fmin2(d__1,d__2);
+	*abserr = *resasc * fmin2(1., pow(*abserr * 200. / *resasc, 1.5));
+    }
     if (*resabs > uflow / (epmach * 50.)) {
-/* Computing MAX */
-	d__1 = epmach * 50. * *resabs;
-	*abserr = fmax2(d__1,*abserr);
+	*abserr = fmax2(epmach * 50. * *resabs, *abserr);
+    }
     return;
 } /* rdqk21_ */
 static void rdqpsrt(int *limit, int *last, int *maxerr,
 		    double *ermax, double *elist, int *iord, int *nrmax)
+{
-    /* System generated locals */
-    int i__1;
     /* Local variables */
-    int ibeg, jbnd, i__, j, k, isucc, jupbn;
+    int i, j, k, ido, jbnd, isucc, jupbn;
     double errmin, errmax;
-    int ido;
-/* ***begin prologue  dqpsrt */
+/* ***begin prologue  dqpsrt
-/* ***refer to  dqage,dqagie,dqagpe,dqawse */
+ ***refer to  dqage,dqagie,dqagpe,dqawse
-/* ***routines called  (none) */
+ ***routines called  (none)
-/* ***revision date  810101   (yymmdd) */
+ ***revision date  810101   (yymmdd)
-/* ***keywords  sequential sorting */
+ ***keywords  sequential sorting
-/* ***author  piessens,robert,appl. math. & progr. div. - k.u.leuven */
+ ***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
-/*           de doncker,elise,appl. math. & progr. div. - k.u.leuven */
+           de doncker,elise,appl. math. & progr. div. - k.u.leuven
-/* ***purpose  this routine maintains the descending ordering in the */
+ ***purpose  this routine maintains the descending ordering in the
-/*            list of the local error estimated resulting from the */
+            list of the local error estimated resulting from the
-/*            interval subdivision process. at each call two error */
+            interval subdivision process. at each call two error
-/*            estimates are inserted using the sequential search */
+            estimates are inserted using the sequential search
-/*            method, top-down for the largest error estimate and */
+            method, top-down for the largest error estimate and
-/*            bottom-up for the smallest error estimate. */
+            bottom-up for the smallest error estimate.
-/* ***description */
+ ***description
-/*           ordering routine */
+           ordering routine
-/*           standard fortran subroutine */
+           standard fortran subroutine
-/*           double precision version */
+           double precision version
-/*           parameters (meaning at output) */
+           parameters (meaning at output)
-/*              limit  - int */
+              limit  - int
-/*                       maximum number of error estimates the list */
+                       maximum number of error estimates the list
-/*                       can contain */
+                       can contain
-/*              last   - int */
+              last   - int
-/*                       number of error estimates currently in the list */
+                       number of error estimates currently in the list
-/*              maxerr - int */
+              maxerr - int
-/*                       maxerr points to the nrmax-th largest error */
+                       maxerr points to the nrmax-th largest error
-/*                       estimate currently in the list */
+                       estimate currently in the list
-/*              ermax  - double precision */
+              ermax  - double precision
-/*                       nrmax-th largest error estimate */
+                       nrmax-th largest error estimate
-/*                       ermax = elist(maxerr) */
+                       ermax = elist(maxerr)
-/*              elist  - double precision */
+              elist  - double precision
-/*                       vector of dimension last containing */
+                       vector of dimension last containing
-/*                       the error estimates */
+                       the error estimates
-/*              iord   - int */
+              iord   - int
-/*                       vector of dimension last, the first k elements */
+                       vector of dimension last, the first k elements
-/*                       of which contain pointers to the error */
+                       of which contain pointers to the error
-/*                       estimates, such that */
+                       estimates, such that
-/*                       elist(iord(1)),...,  elist(iord(k)) */
+                       elist(iord(1)),...,  elist(iord(k))
-/*                       form a decreasing sequence, with */
+                       form a decreasing sequence, with
-/*                       k = last if last.le.(limit/2+2), and */
+                       k = last if last <= (limit/2+2), and
-/*                       k = limit+1-last otherwise */
+                       k = limit+1-last otherwise
-/*              nrmax  - int */
-/*                       maxerr = iord(nrmax) */
-/* ***end prologue  dqpsrt */
+              nrmax  - int
+                       maxerr = iord(nrmax)
+***end prologue  dqpsrt
+*/
-/*           check whether the list contains more than */
-/*           two error estimates. */
-/* ***first executable statement  dqpsrt */
     /* Parameter adjustments */
     --iord;
     --elist;
     /* Function Body */
-    if (*last > 2) {
-	goto L10;
-    }
-    iord[1] = 1;
-    iord[2] = 2;
-    goto L90;
+/*           check whether the list contains more than
+	     two error estimates. */
+    if (*last <= 2) {
+	iord[1] = 1;
+	iord[2] = 2;
+	goto Last;
+    }
-/*           this part of the routine is only executed if, due to a */
+/*           this part of the routine is only executed if, due to a
-/*           difficult integrand, subdivision increased the error */
+	     difficult integrand, subdivision increased the error
-/*           estimate. in the normal case the insert procedure should */
+	     estimate. in the normal case the insert procedure should
-/*           start after the nrmax-th largest error estimate. */
+	     start after the nrmax-th largest error estimate. */
-L10:
     errmax = elist[*maxerr];
-    if (*nrmax == 1) {
+    if (*nrmax > 1) {
-	goto L30;
-    }
-    ido = *nrmax - 1;
+	ido = *nrmax - 1;
-    i__1 = ido;
-    for (i__ = 1; i__ <= i__1; ++i__) {
+	for (i = 1; i <= ido; ++i) {
-	isucc = iord[*nrmax - 1];
+	    isucc = iord[*nrmax - 1];
+	    if (errmax <= elist[isucc])
-/* ***jump out of do-loop */
+		break; /* out of for-loop */
-	if (errmax <= elist[isucc]) {
+	    iord[*nrmax] = isucc;
+	    --(*nrmax);
-	    goto L30;
+	    /* L20: */
+	}
-	iord[*nrmax] = isucc;
-	--(*nrmax);
-/* L20: */
+    }
-/*           compute the number of elements in the list to be maintained */
+/*L30:       compute the number of elements in the list to be maintained
-/*           in descending order. this number depends on the number of */
+	     in descending order. this number depends on the number of
-/*           subdivisions still allowed. */
+	     subdivisions still allowed. */
-L30:
-    jupbn = *last;
-    if (*last > *limit / 2 + 2) {
+    if (*last > *limit / 2 + 2)
 	jupbn = *limit + 3 - *last;
-    }
+    else
+	jupbn = *last;
     errmin = elist[*last];
-/*           insert errmax by traversing the list top-down, */
+/*           insert errmax by traversing the list top-down,
-/*           starting comparison from the element elist(iord(nrmax+1)). */
+	     starting comparison from the element elist(iord(nrmax+1)). */
     jbnd = jupbn - 1;
-    ibeg = *nrmax + 1;
+    for (i = *nrmax + 1; i <= jbnd; ++i) {
-    if (ibeg > jbnd) {
+	isucc = iord[i];
-	goto L50;
+	if (errmax >= elist[isucc]) {/* ***jump out of do-loop */
-    }
+	    /* L60: insert errmin by traversing the list bottom-up. */
-    i__1 = jbnd;
+	    iord[i - 1] = *maxerr;
-    for (i__ = ibeg; i__ <= i__1; ++i__) {
+	    for (j = i, k = jbnd; j <= jbnd; j++, k--) {
-	isucc = iord[i__];
+		isucc = iord[k];
+		if (errmin < elist[isucc]) {
-/* ***jump out of do-loop */
+		    /* goto L80; ***jump out of do-loop */
+		    iord[k + 1] = *last;
+		    goto Last;
+		}
-	if (errmax >= elist[isucc]) {
+		iord[k + 1] = isucc;
+	    }
+	    iord[i] = *last;
-	    goto L60;
+	    goto Last;
+	}
-	iord[i__ - 1] = isucc;
+	iord[i - 1] = isucc;
-/* L40: */
+    }
-L50:
     iord[jbnd] = *maxerr;
     iord[jupbn] = *last;
-    goto L90;
-/*           insert errmin by traversing the list bottom-up. */
+Last:/* set maxerr and ermax. */
-L60:
-    iord[i__ - 1] = *maxerr;
-    k = jbnd;
-    i__1 = jbnd;
-    for (j = i__; j <= i__1; ++j) {
-	isucc = iord[k];
-/* ***jump out of do-loop */
-	if (errmin < elist[isucc]) {
-	    goto L80;
-	}
-	iord[k + 1] = isucc;
-	--k;
-/* L70: */
-    }
-    iord[i__] = *last;
-    goto L90;
-L80:
-    iord[k + 1] = *last;
-/*           set maxerr and ermax. */
-L90:
     *maxerr = iord[*nrmax];
     *ermax = elist[*maxerr];
     return;
 } /* rdqpsrt_ */

The R Project SVN R

(root)/trunk/src/appl/integrate.c – Rev 13845 → 19234