Changeset 180

Timestamp:

03/11/08 04:48:03 (9 months ago)

Author:

Don Clugston

Message:

Added most of the test cases from Alefeld's paper. This clearly shows that the remaining instances of poor performance occur when a successful cubic interpolation forces secant interpolation to be used from that point on.

Files:

• trunk/mathextra/Bracket.d (modified) (16 diffs)

trunk/mathextra/Bracket.d

@@ -1 +1 @@
-/** Algorithms for finding roots and extrema of functions using bracketing.
- *
-
+alefeld's 
- * reverting to parabolic or secant interpolation where necessary.
- * This implementation improves worst-case performance by a factor of more than
- * 100, by performing bisection in _binary space_ (instead of a simple bisection)
- * when slow convergence is encountered.
+ *
+ *  The alefeld algorithm performs badly when:
+ *    - the range in binary space is so large that a naive bisection causes
+ *       catastrophic cancellation; or
+ *    - after a cubic fit, two points have the same y value, so that it then uses secant every time.
- */
-module Bracket;
@@ -80 +85 @@
-unittest{
-    int numCalls;
+    int numProblems=0;
-    
-    void testbinaryChop(real delegate(real) f, real x1, real x2) {
@@ -101 +107 @@
-    void testToms(real delegate(real) f, real x1, real x2) {
-        numCalls=0;
+        ++numProblems;
-        auto result = toms_solve!(real)(f, x1, x2, f(x1), f(x2),(real a, real b){ 
-         return b==nextUp(a) || a==nextUp(b); }, 300 );
-        printf("TOMS Num calls = %d\n",numCalls);
-
-
-
-
-
+
+
+
+
+
+
+
-        }
-*/      
-    }
-    void testBoth(real delegate(real) f, real x1, real x2) {
@@ -158 +166 @@
-    // IE this is 3X faster in the worst case
-    n=3;
+    numProblems=0;
-    foreach(k; nvals) {
-        n=k;
-        testToms(&power, -1, 10);
-    }
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-}
-
@@ -168 +293 @@
-/*
-
-(real x)(
-real q =  sin(x) - x/2;
-for (int i=1; i<20; ++i) q+=(2*i-5.0)*(2*i-5.0)/((x-i*i)*(x-i*i)*(x-i*i));
-return q;
-}
-range PI/2, PI
-n*n + 1e-9 .. (n+1)*(n+1) - 1e-9. n=1(1)19.
@@ -320 +440 @@
-/+  
-SolveResult!(T) toms_solve(T)(T delegate(T) f, T ax, T bx, T fax, T fbx, bool delegate(T,T) tol, uint max_iter)
-         SUBROUTINE RROOT(NPROB,NEPS,EPS,A,B,ROOT)
-         INTEGER NPROB,ITNUM,ISIGN,NEPS
-         DOUBLE PRECISION A,B,FA,FB,C,U,FU,MU,A0,B0,TOL,D,FD
-         DOUBLE PRECISION PROF,E,FE,EPS,ROOT
-         EXTERNAL ISIGN
-         PARAMETER (MU=0.5D0)
-
-// Initialization. set the number of iteration as 0. call subroutine
@@ -331 +445 @@
-// dumb values for the variables "e" and "fe".
-
-
-
-
-
-
+
+
+
+
-         
+
+    for (;;) {
-// iteration starts. the enclosing interval before executing the
-// iteration is recorded as [a0, b0].
-
-
-
+
-
-// updates the number of iteration.
-
-
-
-
-
-
-
-
-
-
-
+
+
+
-         
-// for the first iteration, secant step is taken.
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-         
-// starting with the second iteration, in the first two steps, either
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-// call subroutine "brackt" to get a shrinked enclosing interval as
@@ -393 +530 @@
-// if the criterion is satisfied or the exact solution is obtained.
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
-}
-+/
@@ -487 +564 @@
-{
-   if ((((a-b)==a)) || ((b-a)==b)) {
+       printf("x");
-       // Catastrophic cancellation
-//       assert(a!=0 && b!=0);
@@ -500 +578 @@
-        return c;
- }
+       printf("s");
-   T tol = T.epsilon * 5;
-   T c = a - (fa / (fb - fa)) * (b - a);
@@ -550 +629 @@
-      // Oops, failure, try a secant step:
-      c = secant_interpolate(a, b, fa, fb);
-
+ else printf("p");
-   return c;
-}
@@ -612 +691 @@
-             // Out of bounds step, fall back to quadratic interpolation:          
-             c = quadratic_interpolate!(T)(a, b, d, fa, fb, fd, 3);
-
+else printf("c");
-      } else {
+          printf("=");
-         c = quadratic_interpolate!(T)(a, b, d, fa, fb, fd, 2);
-      }
@@ -633 +713 @@
-             // Out of bounds step, fall back to quadratic interpolation:          
-             c = quadratic_interpolate!(T)(a, b, d, fa, fb, fd, 3);
-
+else printf("c");
-      } else {
+          printf("=");
-         c = quadratic_interpolate!(T)(a, b, d, fa, fb, fd, 2);
-      }
@@ -661 +742 @@
-//         c = a + (b - a) / 2;
-       if ((a-b)==a || (b-a)==b) {
+      printf("b");
-        // DAC: Using ieeeMean here improves worst-case performance by a factor of ~300.
-        // (from >32800 to ~110 for 80-bit reals).
@@ -673 +755 @@
-            }
-       } else {
+      printf("m");
-            c = a + (b - a) / 2;
-       }       
@@ -691 +774 @@
-      // DAC: Also do a binary chop if we're not within a factor of 2 yet -- ie
-      // if we don't yet know what the exponent is.
-
-
+
+
+
-         continue;
+      printf("B");
+//      if( (a==0 || b==0 || (fabs(a)>0.5*fabs(b) && fabs(b)>0.5*fabs(a))) &&  (b - a) < 0.5 * (b0 - a0))
+//         continue;
-
-      // There is a very nasty case where the range contains zero.