Fwd: quadl test request for Matlab

[John W. Eaton]

On 24-Aug-2011, Dmitri A. Sergatskov wrote:

| With tolerance above 0.25 the recursion never happens because in the
| (is+(i1-i2) == is) condition "is" = 1/eps and condition is true.
| One can force the recursion to happen at least once by doing some dirty trick
| like that:
| 
| 
| --- /usr/share/octave/3.4.0/m/general/quadl.m   2011-04-08
| 11:54:40.000000000 -0500
| +++ quadl.m     2011-08-24 20:43:08.626997473 -0500
| @@ -57,6 +57,8 @@
| 
|  function Q = quadl (f, a, b, tol, trace, varargin)
|   need_warning (1);
| +  global iter;
| +  iter = 0;
|   if (nargin < 4)
|     tol = [];
|   endif
| @@ -142,6 +144,7 @@
|  ##   Walter Gautschi, 08/03/98
| 
|  function Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
| +  global iter;
|   h = (b-a)/2;
|   m = (a+b)/2;
|   alpha = sqrt(2/3);
| @@ -159,7 +162,7 @@
|   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)
| +  if ((iter != 0)&&(is+(i1-i2) == is || mll <= a || b <= mrr))
|     if ((m <= a || b <= m) && need_warning ())
|       warning ("quadl: interval contains no more machine number");
|       warning ("quadl: required tolerance may not be met");
| @@ -170,6 +173,7 @@
|       disp ([a, b-a, Q]);
|     endif
|   else
| +    iter++;
|     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{:})
| 
| 
| I do not claim it to be a proper fix...

I also do not claim that this is a proper fix, but here is a way to do
it without the global variable which could conflict with some user's
global:

diff --git a/scripts/general/quadl.m b/scripts/general/quadl.m
--- a/scripts/general/quadl.m
+++ b/scripts/general/quadl.m
@@ -64,6 +64,7 @@
 
 function q = quadl (f, a, b, tol, trace, varargin)
   need_warning (1);
+  tried_at_least_once (0);
   if (nargin < 4)
     tol = [];
   endif
@@ -166,7 +167,7 @@
   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)
+  if ((is+(i1-i2) == is || mll <= a || b <= mrr) && tried_at_least_once ())
     if ((m <= a || b <= m) && need_warning ())
       warning ("quadl: interval contains no more machine number");
       warning ("quadl: required tolerance may not be met");
@@ -195,6 +196,14 @@
   endif
 endfunction
 
+function r = tried_at_least_once (v)
+  persistent count = [];
+  if (nargin == 0)
+    r = count++;
+  else
+    count = v;
+  endif
+endfunction
 
 ## basic functionality
 %!assert( quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16)


jwe