--- a/scripts/general/quadl.m	Sun Nov 08 19:58:33 2015 -0800
+++ b/scripts/general/quadl.m	Sun Nov 08 20:07:37 2015 -0800
@@ -21,6 +21,7 @@
 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol})
 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
+## @deftypefnx {Function File} {[@var{q}, @var{nfun}] =} quadl (@dots{})
 ##
 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
 ## an adaptive Lobatto rule.
@@ -32,142 +33,99 @@
 ## @var{a} and @var{b} are the lower and upper limits of integration.  Both
 ## limits must be finite.
 ##
-## The optional argument @var{tol} defines the relative tolerance with which
-## to perform the integration.  The default value is @code{eps}.
+## The optional argument @var{tol} defines the absolute tolerance with which
+## to perform the integration.  The default value is @math{1e-6}.
 ##
 ## The algorithm used by @code{quadl} involves recursively subdividing the
 ## integration interval.  If @var{trace} is defined then for each subinterval
-## display: (1) the left end of the subinterval, (2) the length of the
-## subinterval, (3) the approximation of the integral over the subinterval.
+## display: (1) the total number of function evaluations, (2) the left end of
+## the subinterval, (3) the length of the subinterval, (4) the approximation of
+## the integral over the subinterval.
 ##
 ## Additional arguments @var{p1}, etc., are passed directly to the function
 ## @var{f}.  To use default values for @var{tol} and @var{trace}, one may pass
 ## empty matrices ([]).
 ##
+## The result of the integration is returned in @var{q}.
+##
+## The optional output @var{nfun} indicates the total number of function
+## evaluations performed.
+##
 ## Reference: @nospell{W. Gander and W. Gautschi}, @cite{Adaptive Quadrature -
 ## Revisited}, BIT Vol. 40, No. 1, March 2000, pp. 84--101.
 ## @url{http://www.inf.ethz.ch/personal/gander/}
 ## @seealso{quad, quadv, quadgk, quadcc, trapz, dblquad, triplequad}
 ## @end deftypefn
 
-##   Author: Walter Gautschi
-##   Date: 08/03/98
-##   Reference: Gander, Computermathematik, Birkhaeuser, 1992.
+## Original Author: Walter Gautschi
+## Date: 08/03/98
+## Reference: Gander, Computermathematik, Birkhaeuser, 1992.
 
 ## 2003-08-05 Shai Ayal
 ##   * permission from author to release as GPL
-## 2004-02-10 Paul Kienzle
-##   * renamed to quadl for compatibility
-##   * replace global variable terminate2 with local function need_warning
-##   * add paper ref to docs
 
-function q = quadl (f, a, b, tol = [], trace = false, varargin)
+function [q, nfun] = quadl (f, a, b, tol = [], trace = false, varargin)
 
   if (nargin < 3)
     print_usage ();
   endif
 
   if (isa (a, "single") || isa (b, "single"))
-    myeps = eps ("single");
+    eps = eps ("single");
   else
-    myeps = eps;
+    eps = eps ("double");
   endif
   if (isempty (tol))
-    tol = myeps;
+    tol = 1e-6;
+  elseif (! isscalar (tol) || tol < 0)
+    error ("quadl: TOL must be a scalar >=0");
+  elseif (tol < eps)
+    tol = eps;
   endif
   if (isempty (trace))
     trace = false;
   endif
-  if (tol < myeps)
-    tol = myeps;
-  endif
 
-  ## Track whether recursion has occurred
-  global __quadl_recurse_done__;
-  __quadl_recurse_done__ = false;
-  ## Track whether warning about machine precision has been issued
-  global __quadl_need_warning__;
-  __quadl_need_warning__ = true;
-
-  m = (a+b)/2;
-  h = (b-a)/2;
-  alpha = sqrt (2/3);
-  beta = 1/sqrt (5);
-
-  x1 = .942882415695480;
-  x2 = .641853342345781;
-  x3 = .236383199662150;
-
-  x = [a, m-x1*h, m-alpha*h, m-x2*h, m-beta*h, m-x3*h, m, m+x3*h, ...
-       m+beta*h, m+x2*h, m+alpha*h, m+x1*h, b];
-
-  y = feval (f, x, varargin{:});
+  y = feval (f, [a, b], varargin{:});
+  nfun = 1;
 
   fa = y(1);
-  fb = y(13);
-
-  i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9)));
+  fb = y(2);
 
-  i1 = (h/1470)*(   77*(y(1)+y(13))
-                 + 432*(y(3)+y(11))
-                 + 625*(y(5)+y(9))
-                 + 672*y(7));
+  h = b - a;
 
-  is = h*( .0158271919734802*(y(1)+y(13))
-          +.0942738402188500*(y(2)+y(12))
-          + .155071987336585*(y(3)+y(11))
-          + .188821573960182*(y(4)+y(10))
-          + .199773405226859*(y(5)+y(9))
-          + .224926465333340*(y(6)+y(8))
-          + .242611071901408*y(7));
+  [q, nfun, hmin] = adaptlobstp (f, a, b, fa, fb, Inf, nfun, abs (h),
+                                 tol, trace, varargin{:});
 
-  s = sign (is);
-  if (s == 0)
-    s = 1;
+  if (nfun > 10_000)
+    warning ("quadl: maximum iteration count reached -- possible singular integral");
+  elseif (any (! isfinite (q(:))))
+    warning ("quadl: infinite or NaN function evaluations were returned");
+  elseif (hmin < (b - a) * eps)
+    warning ("quadl: minimum step size reached -- possible singular integral");
   endif
-  erri1 = abs (i1-is);
-  erri2 = abs (i2-is);
-  if (erri2 != 0)
-    R = erri1/erri2;
-  else
-    R = 1;
-  endif
-  if (R > 0 && R < 1)
-    tol /= R;
-  endif
-  is = s * abs (is) * tol/myeps;
-  if (is == 0)
-    is = b-a;
-  endif
-
-  q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:});
 
 endfunction
 
-## ADAPTLOBSTP  Recursive function used by QUADL.
-##
-##   Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to
-##   approximate the integral of F(X) from A to B to
-##   an appropriate relative error.  The argument 'F' is
-##   a string containing the name of f.  The remaining
-##   arguments are generated by ADAPTLOB or by recursion.
-##
-##   Walter Gautschi, 08/03/98
+function [q, nfun, hmin] = adaptlobstp (f, a, b, fa, fb, q0, nfun, hmin,
+                                        tol, trace, varargin)
+  persistent alpha = sqrt (2/3);
+  persistent beta = 1 / sqrt (5);
 
-function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
-  global __quadl_recurse_done__;
-  global __quadl_need_warning__;
+  if (nfun > 10_000)
+    q = q0;
+    return;
+  endif
 
-  h = (b-a)/2;
-  m = (a+b)/2;
-  alpha = sqrt (2/3);
-  beta = 1 / sqrt (5);
-  mll = m-alpha*h;
-  ml  = m-beta*h;
-  mr  = m+beta*h;
-  mrr = m+alpha*h;
+  h = (b - a) / 2;
+  m = (a + b) / 2;
+  mll = m - alpha*h;
+  ml  = m - beta*h;
+  mr  = m + beta*h;
+  mrr = m + alpha*h;
   x = [mll, ml, m, mr, mrr];
   y = feval (f, x, varargin{:});
+  nfun += 1;
   fmll = y(1);
   fml  = y(2);
   fm   = y(3);
@@ -175,42 +133,60 @@
   fmrr = y(5);
   i2 = (h/6)*(fa + fb + 5*(fml+fmr));
   i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm);
-  if ((is+(i1-i2) == is || mll <= a || b <= mrr) && __quadl_recurse_done__)
-    if ((m <= a || b <= m) && __quadl_need_warning__)
-      warning ("quadl: interval contains no more machine number");
-      warning ("quadl: required tolerance may not be met");
-      __quadl_need_warning__ = false;
-    endif
+
+  if (abs (b - a) < hmin)
+    hmin = abs (b - a);
+  endif
+
+  if (trace)
+    disp ([nfun, a, b-a, i1]);
+  endif
+
+  ## Force at least one adaptive step (nfun > 2 test).
+  if ((abs (i1-i2) < tol || mll <= a || b <= mrr) && nfun > 2)
     q = i1;
-    if (trace)
-      disp ([a, b-a, q]);
-    endif
   else
-    __quadl_recurse_done__ = true;
-    q = (  adaptlobstp (f, a  , mll, fa  , fmll, is, trace, varargin{:})
-         + adaptlobstp (f, mll, ml , fmll, fml , is, trace, varargin{:})
-         + adaptlobstp (f, ml , m  , fml , fm  , is, trace, varargin{:})
-         + adaptlobstp (f, m  , mr , fm  , fmr , is, trace, varargin{:})
-         + adaptlobstp (f, mr , mrr, fmr , fmrr, is, trace, varargin{:})
-         + adaptlobstp (f, mrr, b  , fmrr, fb  , is, trace, varargin{:}));
+    q = zeros (6, 1);
+    [q(1), nfun, hmin] = adaptlobstp (f, a  , mll, fa  , fmll, q0/6, nfun, hmin,
+                                      tol, trace, varargin{:});
+    [q(2), nfun, hmin] = adaptlobstp (f, mll, ml , fmll, fml , q0/6, nfun, hmin,
+                                      tol, trace, varargin{:});
+    [q(3), nfun, hmin] = adaptlobstp (f, ml , m  , fml , fm  , q0/6, nfun, hmin,
+                                      tol, trace, varargin{:});
+    [q(4), nfun, hmin] = adaptlobstp (f, m  , mr , fm  , fmr , q0/6, nfun, hmin,
+                                      tol, trace, varargin{:});
+    [q(5), nfun, hmin] = adaptlobstp (f, mr , mrr, fmr , fmrr, q0/6, nfun, hmin,
+                                      tol, trace, varargin{:});
+    [q(6), nfun, hmin] = adaptlobstp (f, mrr, b  , fmrr, fb  , q0/6, nfun, hmin,
+                                      tol, trace, varargin{:});
+    q = sum (q);
   endif
+
 endfunction
 
 
 ## basic functionality
-%!assert (quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16)
+%!assert (quadl (@(x) sin (x), 0, pi), 2, 5e-15)
 
 ## the values here are very high so it may be unavoidable that this fails
-%!assert (quadl (@(x) sin (3*x).*cosh (x).*sinh (x),10,15),
-%!         2.588424538641647e+10, -1.1e-14)
+%!assert (quadl (@(x) sin (3*x).*cosh (x).*sinh (x),10,15, 1e-3),
+%!        2.588424538641647e+10, 1e-3)
 
 ## extra parameters
 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3),
-%!        cos(2)/3 - cos(5)/3, -3e-16)
+%!        cos(2)/3 - cos(5)/3, 1e-15)
 
 ## test different tolerances.
-%!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, []),
-%!        (60 + sin(4) - sin(64))/12, -0.3)
-%!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []),
-%!        (60 + sin(4) - sin(64))/12, -0.1)
+%!test
+%! [q, nfun1] = quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.5, []);
+%! assert (q, (60 + sin(4) - sin(64))/12, 0.5);
+%! [q, nfun2] = quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []);
+%! assert (q, (60 + sin(4) - sin(64))/12, 0.1);
+%! assert (nfun2 > nfun1);
 
+## Test input validation
+%!error quadl ()
+%!error quadl (@sin)
+%!error quadl (@sin,1)
+%!error <TOL must be a scalar> quadl (@sin, 0, 1, ones (2,2))
+