Mailing-List: contact commits-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: dev@commons.apache.org
Content-Type: text/plain; charset="utf-8"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Subject: svn commit: r1152276 - in /commons/proper/math/trunk/src:
 main/java/org/apache/commons/math/analysis/solvers/ site/xdoc/
 site/xdoc/userguide/ test/java/org/apache/commons/math/analysis/solvers/
Date: Fri, 29 Jul 2011 15:44:23 -0000
To: commits@commons.apache.org
From: luc@apache.org
Message-Id: <20110729154423.5AAF023889BB@eris.apache.org>

Author: luc
Date: Fri Jul 29 15:44:21 2011
New Revision: 1152276

URL: http://svn.apache.org/viewvc?rev=1152276&view=rev
Log:
Added a Brent-like solver that has higher (user specified) order and does
bracket selection on the result: BracketingNthOrderBrentSolver.

JIRA: MATH-635

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java   (with props)
Modified:
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java?rev=1152276&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java Fri Jul 29 15:44:21 2011
@@ -0,0 +1,397 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis.solvers;
+
+
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.commons.math.exception.MathInternalError;
+import org.apache.commons.math.exception.NoBracketingException;
+import org.apache.commons.math.exception.NumberIsTooSmallException;
+import org.apache.commons.math.util.FastMath;
+import org.apache.commons.math.util.MathUtils;
+
+/**
+ * This class implements a modification of the <a
+ * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>.
+ * <p>
+ * The changes with respect to the original Brent algorithm are:
+ * <ul>
+ *   <li>the returned value is chosen in the current interval according
+ *   to user specified {@link AllowedSolutions},</li>
+ *   <li>the maximal order for the invert polynomial root search is
+ *   user-specified instead of being invert quadratic only</li>
+ * </ul>
+ * </p>
+ * The given interval must bracket the root.
+ *
+ * @version $Id$
+ */
+public class BracketingNthOrderBrentSolver
+    extends AbstractUnivariateRealSolver
+    implements BracketedUnivariateRealSolver<UnivariateRealFunction> {
+
+    /** Default absolute accuracy. */
+    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
+
+    /** Default maximal order. */
+    private static final int DEFAULT_MAXIMAL_ORDER = 5;
+
+    /** Maximal aging triggering an attempt to balance the bracketing interval. */
+    private static final int MAXIMAL_AGING = 2;
+
+    /** Reduction factor for attempts to balance the bracketing interval. */
+    private static final double REDUCTION_FACTOR = 1.0 / 16.0;
+
+    /** Maximal order. */
+    private final int maximalOrder;
+
+    /** The kinds of solutions that the algorithm may accept. */
+    private AllowedSolutions allowed;
+
+    /**
+     * Construct a solver with default accuracy and maximal order (1e-6 and 5 respectively)
+     */
+    public BracketingNthOrderBrentSolver() {
+        this(DEFAULT_ABSOLUTE_ACCURACY, DEFAULT_MAXIMAL_ORDER);
+    }
+
+    /**
+     * Construct a solver.
+     *
+     * @param absoluteAccuracy Absolute accuracy.
+     * @param maximalOrder maximal order.
+     * @exception NumberIsTooSmallException if maximal order is lower than 2
+     */
+    public BracketingNthOrderBrentSolver(final double absoluteAccuracy,
+                                         final int maximalOrder)
+        throws NumberIsTooSmallException {
+        super(absoluteAccuracy);
+        if (maximalOrder < 2) {
+            throw new NumberIsTooSmallException(maximalOrder, 2, true);
+        }
+        this.maximalOrder = maximalOrder;
+        this.allowed = AllowedSolutions.ANY_SIDE;
+    }
+
+    /**
+     * Construct a solver.
+     *
+     * @param relativeAccuracy Relative accuracy.
+     * @param absoluteAccuracy Absolute accuracy.
+     * @param maximalOrder maximal order.
+     * @exception NumberIsTooSmallException if maximal order is lower than 2
+     */
+    public BracketingNthOrderBrentSolver(final double relativeAccuracy,
+                                         final double absoluteAccuracy,
+                                         final int maximalOrder)
+        throws NumberIsTooSmallException {
+        super(relativeAccuracy, absoluteAccuracy);
+        if (maximalOrder < 2) {
+            throw new NumberIsTooSmallException(maximalOrder, 2, true);
+        }
+        this.maximalOrder = maximalOrder;
+        this.allowed = AllowedSolutions.ANY_SIDE;
+    }
+
+    /**
+     * Construct a solver.
+     *
+     * @param relativeAccuracy Relative accuracy.
+     * @param absoluteAccuracy Absolute accuracy.
+     * @param functionValueAccuracy Function value accuracy.
+     * @param maximalOrder maximal order.
+     * @exception NumberIsTooSmallException if maximal order is lower than 2
+     */
+    public BracketingNthOrderBrentSolver(final double relativeAccuracy,
+                                         final double absoluteAccuracy,
+                                         final double functionValueAccuracy,
+                                         final int maximalOrder)
+        throws NumberIsTooSmallException {
+        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
+        if (maximalOrder < 2) {
+            throw new NumberIsTooSmallException(maximalOrder, 2, true);
+        }
+        this.maximalOrder = maximalOrder;
+        this.allowed = AllowedSolutions.ANY_SIDE;
+    }
+
+    /** Get the maximal order.
+     * @return maximal order
+     */
+    public int getMaximalOrder() {
+        return maximalOrder;
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    protected double doSolve() {
+
+        // prepare arrays with the first points
+        final double[] x = new double[maximalOrder + 1];
+        final double[] y = new double[maximalOrder + 1];
+        x[0] = getMin();
+        x[1] = getStartValue();
+        x[2] = getMax();
+        verifySequence(x[0], x[1], x[2]);
+
+        // evaluate initial guess
+        y[1] = computeObjectiveValue(x[1]);
+        if (MathUtils.equals(y[1], 0.0, 1)) {
+            // return the initial guess if it is a perfect root.
+            return x[1];
+        }
+
+        // evaluate first  endpoint
+        y[0] = computeObjectiveValue(x[0]);
+        if (MathUtils.equals(y[0], 0.0, 1)) {
+            // return the first endpoint if it is a perfect root.
+            return x[0];
+        }
+
+        int nbPoints;
+        int signChangeIndex;
+        if (y[0] * y[1] < 0) {
+
+            // reduce interval if it brackets the root
+            nbPoints        = 2;
+            signChangeIndex = 1;
+
+        } else {
+
+            // evaluate second endpoint
+            y[2] = computeObjectiveValue(x[2]);
+            if (MathUtils.equals(y[2], 0.0, 1)) {
+                // return the second endpoint if it is a perfect root.
+                return x[2];
+            }
+
+            if (y[1] * y[2] < 0) {
+                // use all computed point as a start sampling array for solving
+                nbPoints        = 3;
+                signChangeIndex = 2;
+            } else {
+                throw new NoBracketingException(x[0], x[2], y[0], y[2]);
+            }
+
+        }
+
+        // prepare a work array for inverse polynomial interpolation
+        final double[] tmpX = new double[x.length];
+
+        // current tightest bracketing of the root
+        double xA    = x[signChangeIndex - 1];
+        double yA    = y[signChangeIndex - 1];
+        double absYA = FastMath.abs(yA);
+        int agingA   = 0;
+        double xB    = x[signChangeIndex];
+        double yB    = y[signChangeIndex];
+        double absYB = FastMath.abs(yB);
+        int agingB   = 0;
+
+        // search loop
+        while (true) {
+
+            // check convergence of bracketing interval
+            final double xTol = getAbsoluteAccuracy() +
+                                getRelativeAccuracy() * FastMath.max(FastMath.abs(xA), FastMath.abs(xB));
+            if (((xB - xA) <= xTol) || (FastMath.max(absYA, absYB) < getFunctionValueAccuracy())) {
+                switch (allowed) {
+                case ANY_SIDE :
+                    return absYA < absYB ? xA : xB;
+                case LEFT_SIDE :
+                    return xA;
+                case RIGHT_SIDE :
+                    return xB;
+                case BELOW_SIDE :
+                    return (yA <= 0) ? xA : xB;
+                case ABOVE_SIDE :
+                    return (yA <  0) ? xB : xA;
+                default :
+                    // this should never happen
+                    throw new MathInternalError(null);
+                }
+            }
+
+            // target for the next evaluation point
+            double targetY;
+            if (agingA >= MAXIMAL_AGING) {
+                // we keep updating the high bracket, try to compensate this
+                targetY = -REDUCTION_FACTOR * yB;
+            } else if (agingB >= MAXIMAL_AGING) {
+                // we keep updating the low bracket, try to compensate this
+                targetY = -REDUCTION_FACTOR * yA;
+            } else {
+                // bracketing is balanced, try to find the root itself
+                targetY = 0;
+            }
+
+            // make a few attempts to guess a root,
+            double nextX;
+            int start = 0;
+            int end   = nbPoints;
+            do {
+
+                // guess a value for current target, using inverse polynomial interpolation
+                System.arraycopy(x, start, tmpX, start, end - start);
+                nextX = guessX(targetY, tmpX, y, start, end);
+
+                if (!((nextX > xA) && (nextX < xB))) {
+                    // the guessed root is not strictly inside of the tightest bracketing interval
+
+                    // the guessed root is either not strictly inside the interval or it
+                    // is a NaN (which occurs when some sampling points share the same y)
+                    // we try again with a lower interpolation order
+                    if (signChangeIndex - start >= end - signChangeIndex) {
+                        // we have more points before the sign change, drop the lowest point
+                        ++start;
+                    } else {
+                        // we have more points after sign change, drop the highest point
+                        --end;
+                    }
+
+                    // we need to do one more attempt
+                    nextX = Double.NaN;
+
+                }
+
+            } while (Double.isNaN(nextX) && (end - start > 1));
+
+            if (Double.isNaN(nextX)) {
+                // fall back to bisection
+                nextX = xA + 0.5 * (xB - xA);
+                start = signChangeIndex - 1;
+                end   = signChangeIndex;
+            }
+
+            // evaluate the function at the guessed root
+            final double nextY = computeObjectiveValue(nextX);
+            if (MathUtils.equals(nextY, 0.0, 1)) {
+                // we have found an exact root, since it is not an approximation
+                // we don't need to bother about the allowed solutions setting
+                return nextX;
+            }
+
+            if ((nbPoints > 2) && (end - start != nbPoints)) {
+
+                // we have been forced to ignore some points to keep bracketing,
+                // they are probably too far from the root, drop them from now on
+                nbPoints = end - start;
+                System.arraycopy(x, start, x, 0, nbPoints);
+                System.arraycopy(y, start, y, 0, nbPoints);
+                signChangeIndex -= start;
+
+            } else  if (nbPoints == x.length) {
+
+                // we have to drop one point in order to insert the new one
+                nbPoints--;
+
+                // keep the tightest bracketing interval as centered as possible
+                if (signChangeIndex >= (x.length + 1) / 2) {
+                    // we drop the lowest point, we have to shift the arrays and the index
+                    System.arraycopy(x, 1, x, 0, nbPoints);
+                    System.arraycopy(y, 1, y, 0, nbPoints);
+                    --signChangeIndex;
+                }
+
+            }
+
+            // insert the last computed point
+            //(by construction, we know it lies inside the tightest bracketing interval)
+            System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex);
+            x[signChangeIndex] = nextX;
+            System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex);
+            y[signChangeIndex] = nextY;
+            ++nbPoints;
+
+            // update the bracketing interval
+            if (nextY * yA <= 0) {
+                // the sign change occurs before the inserted point
+                xB = nextX;
+                yB = nextY;
+                absYB = FastMath.abs(yB);
+                ++agingA;
+                agingB = 0;
+            } else {
+                // the sign change occurs after the inserted point
+                xA = nextX;
+                yA = nextY;
+                absYA = FastMath.abs(yA);
+                agingA = 0;
+                ++agingB;
+
+                // update the sign change index
+                signChangeIndex++;
+
+            }
+
+        }
+
+    }
+
+    /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation.
+     * <p>
+     * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q
+     * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>),
+     * Q(y<sub>i</sub>) = x<sub>i</sub>.
+     * </p>
+     * @param targetY target value for y
+     * @param x reference points abscissas for interpolation,
+     * note that this array <em>is</em> modified during computation
+     * @param y reference points ordinates for interpolation
+     * @param start start index of the points to consider (inclusive)
+     * @param end end index of the points to consider (exclusive)
+     * @return guessed root (will be a NaN if two points share the same y)
+     */
+    private double guessX(final double targetY, final double[] x, final double[] y,
+                          final int start, final int end) {
+
+        // compute Q Newton coefficients by divided differences
+        for (int i = start; i < end - 1; ++i) {
+            final int delta = i + 1 - start;
+            for (int j = end - 1; j > i; --j) {
+                x[j] = (x[j] - x[j-1]) / (y[j] - y[j - delta]);
+            }
+        }
+
+        // evaluate Q(targetY)
+        double x0 = 0;
+        for (int j = end - 1; j >= start; --j) {
+            x0 = x[j] + x0 * (targetY - y[j]);
+        }
+
+        return x0;
+
+    }
+
+    /** {@inheritDoc} */
+    public double solve(int maxEval, UnivariateRealFunction f, double min,
+                        double max, AllowedSolutions allowedSolutions) {
+        this.allowed = allowedSolutions;
+        return super.solve(maxEval, f, min, max);
+    }
+
+    /** {@inheritDoc} */
+    public double solve(int maxEval, UnivariateRealFunction f, double min,
+                        double max, double startValue,
+                        AllowedSolutions allowedSolutions) {
+        this.allowed = allowedSolutions;
+        return super.solve(maxEval, f, min, max, startValue);
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1152276&r1=1152275&r2=1152276&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Jul 29 15:44:21 2011
@@ -52,6 +52,10 @@ The <action> type attribute can be add,u
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="luc" type="add" issue="MATH-635" >
+        Added a Brent-like solver that has higher (user specified) order and does
+        bracket selection on the result: BracketingNthOrderBrentSolver.
+      </action>
       <action dev="luc" type="add" >
         Added a few shortcut methods and predicates to Dfp (abs, isZero,
         negativeOrNull, strictlyNegative, positiveOrNull, strictlyPositive).

Modified: commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml?rev=1152276&r1=1152275&r2=1152276&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml Fri Jul 29 15:44:21 2011
@@ -74,6 +74,12 @@
               <td>no</td>
             </tr>
             <tr>
+              <td><a href="../apidocs/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolver.html">bracketing n<sup>th</sup> order Brent</a></td>
+              <td>variable order, guaranteed</td>
+              <td>yes</td>
+              <td>yes</td>
+            </tr>
+            <tr>
               <td><a href="../apidocs/org/apache/commons/math/analysis/solvers/IllinoisSolver.html">Illinois Method</a></td>
               <td><a href="../apidocs/org/apache/commons/math/analysis/UnivariateRealFunction.html">univariate real-valued functions</a></td>
               <td>super-linear, guaranteed</td>
@@ -206,7 +212,8 @@
         <source>UnivariateRealFunction function = // some user defined function object
 final double relativeAccuracy = 1.0e-12;
 final double absoluteAccuracy = 1.0e-8;
-UnivariateRealSolver solver   = new PegasusSolver(relativeAccuracy, absoluteAccuracy);
+final int    maxOrder         = 5;
+UnivariateRealSolver solver   = new BracketingNthOrderBrentSolver(relativeAccuracy, absoluteAccuracy, maxOrder);
 double c = solver.solve(100, function, 1.0, 5.0, AllowedSolutions.LEFT_SIDE);</source>
         <p>
           Force bracketing, by refining a base solution found by a non-bracketing solver:
@@ -221,9 +228,9 @@ double c = UnivariateRealSolverUtils.for
                                                baseRoot, 1.0, 5.0, AllowedSolutions.LEFT_SIDE);
 </source>
         <p>
-          The <code>BrentSolve</code> uses the Brent-Dekker algorithm which is
-          fast and robust.  This algorithm is recommended for most users.  If there
-          are multiple roots in the interval, or there is a large domain of indeterminacy, the
+          The <code>BrentSolver</code> uses the Brent-Dekker algorithm which is
+          fast and robust.  If there are multiple roots in the interval,
+          or there is a large domain of indeterminacy, the
           algorithm will converge to a random root in the interval without
           indication that there are problems.  Interestingly, the examined text
           book implementations all disagree in details of the convergence
@@ -233,6 +240,15 @@ double c = UnivariateRealSolverUtils.for
           algorithm.
         </p>
         <p>
+          The <code>BracketingNthOrderBrentSolver</code> uses an extension of the
+          Brent-Dekker algorithm which uses inverse n<sup>th</sup> order polynomial
+          interpolation instead of inverse quadratic interpolation, and which allows
+          selection of the side of the convergence interval for result bracketing.
+          This is now the recommended algorithm for most users since it has the
+          largest order, doesn't require derivatives, has guaranteed convergence
+          and allows result bracket selection.
+        </p>
+        <p>
           The <code>SecantSolver</code> uses a straightforward secant
           algorithm which does not bracket the search and therefore does not
           guarantee convergence.  It may be faster than Brent on some well-behaved

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java?rev=1152276&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java Fri Jul 29 15:44:21 2011
@@ -0,0 +1,181 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.analysis.solvers;
+
+import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
+import org.apache.commons.math.analysis.QuinticFunction;
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.commons.math.exception.NumberIsTooSmallException;
+import org.apache.commons.math.exception.TooManyEvaluationsException;
+import org.apache.commons.math.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test case for {@link BracketingNthOrderBrentSolver bracketing n<sup>th</sup> order Brent} solver.
+ *
+ * @version $Id$
+ */
+public final class BracketingNthOrderBrentSolverTest extends BaseSecantSolverAbstractTest {
+    /** {@inheritDoc} */
+    protected UnivariateRealSolver getSolver() {
+        return new BracketingNthOrderBrentSolver();
+    }
+
+    /** {@inheritDoc} */
+    protected int[] getQuinticEvalCounts() {
+        return new int[] {1, 3, 8, 1, 9, 4, 8, 1, 12, 1, 14};
+    }
+
+    @Test(expected=NumberIsTooSmallException.class)
+    public void testInsufficientOrder1() {
+        new BracketingNthOrderBrentSolver(1.0e-10, 1);
+    }
+
+    @Test(expected=NumberIsTooSmallException.class)
+    public void testInsufficientOrder2() {
+        new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1);
+    }
+
+    @Test(expected=NumberIsTooSmallException.class)
+    public void testInsufficientOrder3() {
+        new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1.0e-10, 1);
+    }
+
+    @Test
+    public void testConstructorsOK() {
+        Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 2).getMaximalOrder());
+        Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 2).getMaximalOrder());
+        Assert.assertEquals(2, new BracketingNthOrderBrentSolver(1.0e-10, 1.0e-10, 1.0e-10, 2).getMaximalOrder());
+    }
+
+    @Test
+    public void testConvergenceOnFunctionAccuracy() {
+        BracketingNthOrderBrentSolver solver =
+                new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-10, 0.001, 3);
+        QuinticFunction f = new QuinticFunction();
+        double result = solver.solve(20, f, 0.2, 0.9, 0.4, AllowedSolutions.BELOW_SIDE);
+        Assert.assertEquals(0, f.value(result), solver.getFunctionValueAccuracy());
+        Assert.assertTrue(f.value(result) <= 0);
+        Assert.assertTrue(result - 0.5 > solver.getAbsoluteAccuracy());
+        result = solver.solve(20, f, -0.9, -0.2,  -0.4, AllowedSolutions.ABOVE_SIDE);
+        Assert.assertEquals(0, f.value(result), solver.getFunctionValueAccuracy());
+        Assert.assertTrue(f.value(result) >= 0);
+        Assert.assertTrue(result + 0.5 < -solver.getAbsoluteAccuracy());
+    }
+
+    @Test
+    public void testFasterThanNewton() {
+        // the following test functions come from Beny Neta's paper:
+        // "Several New Methods for solving Equations"
+        // intern J. Computer Math Vol 23 pp 265-282
+        // available here: http://www.math.nps.navy.mil/~bneta/SeveralNewMethods.PDF
+        // the reference roots have been computed by the Dfp solver to more than
+        // 80 digits and checked with emacs (only the first 20 digits are reproduced here)
+        compare(new TestFunction(0.0, -2, 2) {
+            public double value(double x)      { return FastMath.sin(x) - 0.5 * x; }
+            public double derivative(double x) { return FastMath.cos(x) - 0.5; }
+        });
+        compare(new TestFunction(6.3087771299726890947, -5, 10) {
+            public double value(double x)      { return FastMath.pow(x, 5) + x - 10000; }
+            public double derivative(double x) { return 5 * FastMath.pow(x, 4) + 1; }
+        });
+        compare(new TestFunction(9.6335955628326951924, 0.001, 10) {
+            public double value(double x)      { return FastMath.sqrt(x) - 1 / x - 3; }
+            public double derivative(double x) { return 0.5 / FastMath.sqrt(x) + 1 / (x * x); }
+        });
+        compare(new TestFunction(2.8424389537844470678, -5, 5) {
+            public double value(double x)      { return FastMath.exp(x) + x - 20; }
+            public double derivative(double x) { return FastMath.exp(x) + 1; }
+        });
+        compare(new TestFunction(8.3094326942315717953, 0.001, 10) {
+            public double value(double x)      { return FastMath.log(x) + FastMath.sqrt(x) - 5; }
+            public double derivative(double x) { return 1 / x + 0.5 / FastMath.sqrt(x); }
+        });
+        compare(new TestFunction(1.4655712318767680266, -0.5, 1.5) {
+            public double value(double x)      { return (x - 1) * x * x - 1; }
+            public double derivative(double x) { return (3 * x - 2) * x; }
+        });
+
+    }
+
+    private void compare(TestFunction f) {
+        compare(f, f.getRoot(), f.getMin(), f.getMax());
+    }
+
+    private void compare(DifferentiableUnivariateRealFunction f,
+                         double root, double min, double max) {
+        NewtonSolver newton = new NewtonSolver(1.0e-12);
+        BracketingNthOrderBrentSolver bracketing =
+                new BracketingNthOrderBrentSolver(1.0e-12, 1.0e-12, 1.0e-18, 5);
+        double resultN;
+        try {
+            resultN = newton.solve(100, f, min, max);
+        } catch (TooManyEvaluationsException tmee) {
+            resultN = Double.NaN;
+        }
+        double resultB;
+        try {
+            resultB = bracketing.solve(100, f, min, max);
+        } catch (TooManyEvaluationsException tmee) {
+            resultB = Double.NaN;
+        }
+        Assert.assertEquals(root, resultN, newton.getAbsoluteAccuracy());
+        Assert.assertEquals(root, resultB, bracketing.getAbsoluteAccuracy());
+        Assert.assertTrue(bracketing.getEvaluations() < newton.getEvaluations());
+    }
+
+    private static abstract class TestFunction implements DifferentiableUnivariateRealFunction {
+
+        private final double root;
+        private final double min;
+        private final double max;
+
+        protected TestFunction(final double root, final double min, final double max) {
+            this.root = root;
+            this.min  = min;
+            this.max  = max;
+        }
+
+        public double getRoot() {
+            return root;
+        }
+
+        public double getMin() {
+            return min;
+        }
+
+        public double getMax() {
+            return max;
+        }
+
+        public abstract double value(double x);
+
+        public abstract double derivative(double x);
+
+        public UnivariateRealFunction derivative() {
+            return new UnivariateRealFunction() {
+                public double value(double x) {
+                     return derivative(x);
+                }
+            };
+        }
+
+    }
+
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BracketingNthOrderBrentSolverTest.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision