commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1039083 [1/4] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/analysis/ main/java/org/apache/commons/math/analysis/polynomials/ main/java/org/apache/commons/math/analysis/solvers/ main/java/org/apache/commons/math/di...
Date Thu, 25 Nov 2010 16:22:01 GMT
Author: erans
Date: Thu Nov 25 16:22:00 2010
New Revision: 1039083

URL: http://svn.apache.org/viewvc?rev=1039083&view=rev
Log:
MATH-439
Refactored the "solvers" package. Implementations refer to number of
evaluation of the objective function (instead of the number of iterations).
New interfaces and base classes.
"NewtonSolver" fits in the design without resorting to a cast.
Created class "MullerSolver2" to contain the code of the method named "solve2"
in class "MullerSolver".
Removed "UnivariateRealSolverFactory" and "UnivariateRealSolverFactoryImpl".
Default solver in "UnivariateRealSolverUtils" is explicitely instantiated.
"AbstractContinuousDistribution": Type of exception thrown changed in
"UnivariateRealSolverUtils".
Factored out duplicate code (in "GaussNewtonOptimizerTest" and
"LevenbergMarquardtOptimizerTest"): class "Circle" is now called
"CircleVectorial". Also factored out the "Circle" class from
"NonLinearConjugateGradientOptimizerTest": class is named "CircleScalar".
Created "SecantSolverTest", moving there all the tests for the class
"SecantSolver" that were located in class "BrentSolverTest".
Created new interface and base class for polynomial functions solvers
("LaguerreSolver") so that the function type is now checked at compile time.
Removed deprecated exceptions (MATH-441).
Javadoc clean-up.
Lowered tolerance values in some unit tests.
Tests upgraded to Junit 4 (MATH-423).


Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractDifferentiableUnivariateRealSolver.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractPolynomialSolver.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractUnivariateRealSolver.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseAbstractUnivariateRealSolver.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseUnivariateRealSolver.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/DifferentiableUnivariateRealSolver.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver2.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/PolynomialSolver.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/NoBracketingException.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/MullerSolver2Test.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/SecantSolverTest.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/general/CircleScalar.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/general/CircleVectorial.java   (with props)
Removed:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverFactory.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverFactoryImpl.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverFactoryImplTest.java
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/FunctionUtils.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionNewtonForm.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BisectionSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/NewtonSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/RiddersSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/SecantSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverImpl.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverUtils.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/LeastSquaresConverter.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/fitting/GaussianFitter.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/fitting/ParametricGaussianFunction.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/NonLinearConjugateGradientOptimizer.java
    commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BisectionSolverTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/BrentSolverTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/LaguerreSolverTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/MullerSolverTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/NewtonSolverTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/RiddersSolverTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverUtilsTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizerTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizerTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/general/NonLinearConjugateGradientOptimizerTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/FunctionUtils.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/FunctionUtils.java?rev=1039083&r1=1039082&r2=1039083&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/FunctionUtils.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/FunctionUtils.java Thu Nov 25 16:22:00 2010
@@ -27,6 +27,11 @@ import org.apache.commons.math.analysis.
  */
 public class FunctionUtils {
     /**
+     * Class only contains static methods.
+     */
+    private FunctionUtils() {}
+
+    /**
      * Compose functions.
      *
      * @param f List of functions.
@@ -142,6 +147,7 @@ public class FunctionUtils {
      *
      * @param f Binary function.
      * @param fixed Value to which the first argument of {@code f} is set.
+     * @return a unary function.
      */
     public static UnivariateRealFunction fix1stArgument(final BivariateRealFunction f,
                                                         final double fixed) {
@@ -157,6 +163,7 @@ public class FunctionUtils {
      *
      * @param f Binary function.
      * @param fixed Value to which the second argument of {@code f} is set.
+     * @return a unary function.
      */
     public static UnivariateRealFunction fix2ndArgument(final BivariateRealFunction f,
                                                         final double fixed) {

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionNewtonForm.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionNewtonForm.java?rev=1039083&r1=1039082&r2=1039083&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionNewtonForm.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionNewtonForm.java Thu Nov 25 16:22:00 2010
@@ -16,7 +16,6 @@
  */
 package org.apache.commons.math.analysis.polynomials;
 
-import org.apache.commons.math.exception.NullArgumentException;
 import org.apache.commons.math.exception.NoDataException;
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.analysis.UnivariateRealFunction;
@@ -69,7 +68,8 @@ public class PolynomialFunctionNewtonFor
      *
      * @param a Coefficients in Newton form formula.
      * @param c Centers.
-     * @throws NullArgumentException if any argument is {@code null}.
+     * @throws org.apache.commons.math.exception.NullArgumentException if
+     * any argument is {@code null}.
      * @throws NoDataException if any array has zero length.
      * @throws DimensionMismatchException if the size difference between
      * {@code a} and {@code c} is not equal to 1.
@@ -154,7 +154,8 @@ public class PolynomialFunctionNewtonFor
      * @param c Centers.
      * @param z Point at which the function value is to be computed.
      * @return the function value.
-     * @throws NullArgumentException if any argument is {@code null}.
+     * @throws org.apache.commons.math.exception.NullArgumentException if
+     * any argument is {@code null}.
      * @throws NoDataException if any array has zero length.
      * @throws DimensionMismatchException if the size difference between
      * {@code a} and {@code c} is not equal to 1.
@@ -202,7 +203,8 @@ public class PolynomialFunctionNewtonFor
      *
      * @param a the coefficients in Newton form formula
      * @param c the centers
-     * @throws NullArgumentException if any argument is {@code null}.
+     * @throws org.apache.commons.math.exception.NullArgumentException if
+     * any argument is {@code null}.
      * @throws NoDataException if any array has zero length.
      * @throws DimensionMismatchException if the size difference between
      * {@code a} and {@code c} is not equal to 1.

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractDifferentiableUnivariateRealSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractDifferentiableUnivariateRealSolver.java?rev=1039083&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractDifferentiableUnivariateRealSolver.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractDifferentiableUnivariateRealSolver.java Thu Nov 25 16:22:00 2010
@@ -0,0 +1,81 @@
+/*
+ * 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.UnivariateRealFunction;
+
+/**
+ * Provide a default implementation for several functions useful to generic
+ * solvers.
+ *
+ * @version $Revision$ $Date$
+ * @since 3.0
+ */
+public abstract class AbstractDifferentiableUnivariateRealSolver
+    extends BaseAbstractUnivariateRealSolver<DifferentiableUnivariateRealFunction>
+    implements DifferentiableUnivariateRealSolver {
+    /** Derivative of the function to solve. */
+    private UnivariateRealFunction functionDerivative;
+
+    /**
+     * Construct a solver with given absolute accuracy.
+     *
+     * @param absoluteAccuracy Maximum absolute error.
+     */
+    protected AbstractDifferentiableUnivariateRealSolver(final double absoluteAccuracy) {
+        super(absoluteAccuracy);
+    }
+
+    /**
+     * Construct a solver with given accuracies.
+     *
+     * @param relativeAccuracy Maximum relative error.
+     * @param absoluteAccuracy Maximum absolute error.
+     * @param functionValueAccuracy Maximum function value error.
+     */
+    protected AbstractDifferentiableUnivariateRealSolver(final double relativeAccuracy,
+                                                         final double absoluteAccuracy,
+                                                         final double functionValueAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
+    }
+
+    /**
+     * Compute the objective function value.
+     *
+     * @param point Point at which the objective function must be evaluated.
+     * @return the objective function value at specified point.
+     * @throws org.apache.commons.math.exception.TooManyEvaluationsException
+     * if the maximal number of evaluations is exceeded.
+     */
+    protected double computeDerivativeObjectiveValue(double point) {
+        incrementEvaluationCount();
+        return functionDerivative.value(point);
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    protected void setup(DifferentiableUnivariateRealFunction f,
+                         double min, double max,
+                         double startValue) {
+        super.setup(f, min, max, startValue);
+        functionDerivative = f.derivative();
+    }
+}

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

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractPolynomialSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractPolynomialSolver.java?rev=1039083&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractPolynomialSolver.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractPolynomialSolver.java Thu Nov 25 16:22:00 2010
@@ -0,0 +1,82 @@
+/*
+ * 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.polynomials.PolynomialFunction;
+
+/**
+ * Base class for solvers.
+ *
+ * @version $Revision$ $Date$
+ * @since 3.0
+ */
+public abstract class AbstractPolynomialSolver
+    extends BaseAbstractUnivariateRealSolver<PolynomialFunction>
+    implements PolynomialSolver {
+    /** Function. */
+    private PolynomialFunction polynomialFunction;
+
+    /**
+     * Construct a solver with given absolute accuracy.
+     *
+     * @param absoluteAccuracy Maximum absolute error.
+     */
+    protected AbstractPolynomialSolver(final double absoluteAccuracy) {
+        super(absoluteAccuracy);
+    }
+    /**
+     * Construct a solver with given accuracies.
+     *
+     * @param relativeAccuracy Maximum relative error.
+     * @param absoluteAccuracy Maximum absolute error.
+     */
+    protected AbstractPolynomialSolver(final double relativeAccuracy,
+                                       final double absoluteAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy);
+    }
+    /**
+     * Construct a solver with given accuracies.
+     *
+     * @param relativeAccuracy Maximum relative error.
+     * @param absoluteAccuracy Maximum absolute error.
+     * @param functionValueAccuracy Maximum function value error.
+     */
+    protected AbstractPolynomialSolver(final double relativeAccuracy,
+                                       final double absoluteAccuracy,
+                                       final double functionValueAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    protected void setup(PolynomialFunction f,
+                         double min, double max,
+                         double startValue) {
+        super.setup(f, min, max, startValue);
+        polynomialFunction = f;
+    }
+
+    /**
+     * @return the coefficients of the polynomial function.
+     */
+    protected double[] getCoefficients() {
+        return polynomialFunction.getCoefficients();
+    }
+}

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

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractUnivariateRealSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractUnivariateRealSolver.java?rev=1039083&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractUnivariateRealSolver.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/AbstractUnivariateRealSolver.java Thu Nov 25 16:22:00 2010
@@ -0,0 +1,61 @@
+/*
+ * 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;
+
+/**
+ * Base class for solvers.
+ *
+ * @version $Revision$ $Date$
+ * @since 3.0
+ */
+public abstract class AbstractUnivariateRealSolver
+    extends BaseAbstractUnivariateRealSolver<UnivariateRealFunction>
+    implements UnivariateRealSolver {
+    /**
+     * Construct a solver with given absolute accuracy.
+     *
+     * @param absoluteAccuracy Maximum absolute error.
+     */
+    protected AbstractUnivariateRealSolver(final double absoluteAccuracy) {
+        super(absoluteAccuracy);
+    }
+    /**
+     * Construct a solver with given accuracies.
+     *
+     * @param relativeAccuracy Maximum relative error.
+     * @param absoluteAccuracy Maximum absolute error.
+     */
+    protected AbstractUnivariateRealSolver(final double relativeAccuracy,
+                                           final double absoluteAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy);
+    }
+    /**
+     * Construct a solver with given accuracies.
+     *
+     * @param relativeAccuracy Maximum relative error.
+     * @param absoluteAccuracy Maximum absolute error.
+     * @param functionValueAccuracy Maximum function value error.
+     */
+    protected AbstractUnivariateRealSolver(final double relativeAccuracy,
+                                           final double absoluteAccuracy,
+                                           final double functionValueAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
+    }
+}

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

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseAbstractUnivariateRealSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseAbstractUnivariateRealSolver.java?rev=1039083&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseAbstractUnivariateRealSolver.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseAbstractUnivariateRealSolver.java Thu Nov 25 16:22:00 2010
@@ -0,0 +1,299 @@
+/*
+ * 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.util.Incrementor;
+import org.apache.commons.math.exception.MaxCountExceededException;
+import org.apache.commons.math.exception.TooManyEvaluationsException;
+import org.apache.commons.math.exception.NullArgumentException;
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+
+/**
+ * Provide a default implementation for several functions useful to generic
+ * solvers.
+ *
+ * @param <FUNC> Type of function to solve.
+ *
+ * @version $Revision: 1030464 $ $Date: 2010-11-03 14:46:04 +0100 (Wed, 03 Nov 2010) $
+ * @since 2.0
+ */
+public abstract class BaseAbstractUnivariateRealSolver<FUNC extends UnivariateRealFunction>
+    implements BaseUnivariateRealSolver<FUNC> {
+    /** Default absolute accuracy */
+    public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
+    /** Default relative accuracy. */
+    public static final double DEFAULT_RELATIVE_ACCURACY = 1e-14;
+    /** Default function value accuracy. */
+    public static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15;
+    /** Function value accuracy. */
+    private final double functionValueAccuracy;
+    /** Absolute accuracy. */
+    private final double absoluteAccuracy;
+    /** Relative accuracy. */
+    private final double relativeAccuracy;
+    /** Evaluations counter. */
+    private final Incrementor evaluations = new Incrementor();
+    /** Lower end of search interval. */
+    private double searchMin;
+    /** Higher end of search interval. */
+    private double searchMax;
+    /** Initial guess. */
+    private double searchStart;
+    /** Function to solve. */
+    private FUNC function;
+
+    /**
+     * Construct a solver with given absolute accuracy.
+     *
+     * @param absoluteAccuracy Maximum absolute error.
+     */
+    protected BaseAbstractUnivariateRealSolver(final double absoluteAccuracy) {
+        this(DEFAULT_RELATIVE_ACCURACY,
+             absoluteAccuracy,
+             DEFAULT_FUNCTION_VALUE_ACCURACY);
+    }
+
+    /**
+     * Construct a solver with given accuracies.
+     *
+     * @param relativeAccuracy Maximum relative error.
+     * @param absoluteAccuracy Maximum absolute error.
+     */
+    protected BaseAbstractUnivariateRealSolver(final double relativeAccuracy,
+                                               final double absoluteAccuracy) {
+        this(relativeAccuracy,
+             absoluteAccuracy,
+             DEFAULT_FUNCTION_VALUE_ACCURACY);
+    }
+
+    /**
+     * Construct a solver with given accuracies.
+     *
+     * @param relativeAccuracy Maximum relative error.
+     * @param absoluteAccuracy Maximum absolute error.
+     * @param functionValueAccuracy Maximum function value error.
+     */
+    protected BaseAbstractUnivariateRealSolver(final double relativeAccuracy,
+                                               final double absoluteAccuracy,
+                                               final double functionValueAccuracy) {
+        this.absoluteAccuracy = absoluteAccuracy;
+        this.relativeAccuracy = relativeAccuracy;
+        this.functionValueAccuracy = functionValueAccuracy;
+    }
+
+    /** {@inheritDoc} */
+    public void setMaxEvaluations(int maxEvaluations) {
+        evaluations.setMaximalCount(maxEvaluations);
+    }
+    /** {@inheritDoc} */
+    public int getMaxEvaluations() {
+        return evaluations.getMaximalCount();
+    }
+    /** {@inheritDoc} */
+    public int getEvaluations() {
+        return evaluations.getCount();
+    }
+    /**
+     * @return the lower end of the search interval.
+     */
+    public double getMin() {
+        return searchMin;
+    }
+    /**
+     * @return the higher end of the search interval.
+     */
+    public double getMax() {
+        return searchMax;
+    }
+    /**
+     * @return the initial guess.
+     */
+    public double getStartValue() {
+        return searchStart;
+    }
+    /**
+     * {@inheritDoc}
+     */
+    public double getAbsoluteAccuracy() {
+        return absoluteAccuracy;
+    }
+    /**
+     * {@inheritDoc}
+     */
+    public double getRelativeAccuracy() {
+        return relativeAccuracy;
+    }
+    /**
+     * {@inheritDoc}
+     */
+    public double getFunctionValueAccuracy() {
+        return functionValueAccuracy;
+    }
+
+    /**
+     * Compute the objective function value.
+     *
+     * @param point Point at which the objective function must be evaluated.
+     * @return the objective function value at specified point.
+     * @throws TooManyEvaluationsException if the maximal number of evaluations
+     * is exceeded.
+     */
+    protected double computeObjectiveValue(double point) {
+        incrementEvaluationCount();
+        return function.value(point);
+    }
+
+    /**
+     * Prepare for computation.
+     * Subclasses must call this method if they override any of the
+     * {@code solve} methods.
+     *
+     * @param f Function to solve.
+     * @param min Lower bound for the interval.
+     * @param max Upper bound for the interval.
+     * @param startValue Start value to use.
+     */
+    protected void setup(FUNC f,
+                         double min, double max,
+                         double startValue) {
+        // Checks.
+        if (f == null) {
+            throw new NullArgumentException();
+        }
+
+        // Reset.
+        searchMin = min;
+        searchMax = max;
+        searchStart = startValue;
+        function = f;
+        evaluations.resetCount();
+    }
+
+    /** {@inheritDoc} */
+    public double solve(FUNC f, double min, double max, double startValue) {
+        // Initialization.
+        setup(f, min, max, startValue);
+
+        // Perform computation.
+        return doSolve();
+    }
+
+    /** {@inheritDoc} */
+    public double solve(FUNC f, double min, double max) {
+        return solve(f, min, max, min + 0.5 * (max - min));
+    }
+
+    /** {@inheritDoc} */
+    public double solve(FUNC f, double startValue) {
+        return solve(f, Double.NaN, Double.NaN, startValue);
+    }
+
+    /**
+     * Method for implementing actual optimization algorithms in derived
+     * classes.
+     *
+     * @return the root.
+     * @throws TooManyEvaluationsException if the maximal number of evaluations
+     * is exceeded.
+     */
+    protected abstract double doSolve();
+
+    /**
+     * Check whether the function takes opposite signs at the endpoints.
+     *
+     * @param lower Lower endpoint.
+     * @param upper Upper endpoint.
+     * @return {@code true} if the function values have opposite signs at the
+     * given points.
+     */
+    protected boolean isBracketing(final double lower,
+                                   final double upper) {
+        return UnivariateRealSolverUtils.isBracketing(function, lower, upper);
+    }
+
+    /**
+     * Check whether the arguments form a (strictly) increasing sequence.
+     *
+     * @param start First number.
+     * @param mid Second number.
+     * @param end Third number.
+     * @return {@code true} if the arguments form an increasing sequence.
+     */
+    protected boolean isSequence(final double start,
+                                 final double mid,
+                                 final double end) {
+        return UnivariateRealSolverUtils.isSequence(start, mid, end);
+    }
+
+    /**
+     * Check that the endpoints specify an interval.
+     *
+     * @param lower Lower endpoint.
+     * @param upper Upper endpoint.
+     * @throws org.apache.commons.math.exception.NumberIsTooLargeException
+     * if {@code lower >= upper}.
+     */
+    protected void verifyInterval(final double lower,
+                                  final double upper) {
+        UnivariateRealSolverUtils.verifyInterval(lower, upper);
+    }
+
+    /**
+     * Check that {@code lower < initial < upper}.
+     *
+     * @param lower Lower endpoint.
+     * @param initial Initial value.
+     * @param upper Upper endpoint.
+     * @throws org.apache.commons.math.exception.NumberIsTooLargeException
+     * if {@code lower >= initial} or {@code initial >= upper}.
+     */
+    protected void verifySequence(final double lower,
+                                  final double initial,
+                                  final double upper) {
+        UnivariateRealSolverUtils.verifySequence(lower, initial, upper);
+    }
+
+    /**
+     * Check that the endpoints specify an interval and the function takes
+     * opposite signs at the endpoints.
+     *
+     * @param lower Lower endpoint.
+     * @param upper Upper endpoint.
+     * @throws org.apache.commons.math.exception.NoBracketingException if
+     * the function has the same sign at the endpoints.
+     */
+    protected void verifyBracketing(final double lower,
+                                    final double upper) {
+        UnivariateRealSolverUtils.verifyBracketing(function, lower, upper);
+    }
+
+    /**
+     * Increment the evaluation count by one.
+     * Method {@link #computeObjectiveValue(double)} calls this method internally.
+     * It is provided for subclasses that do not exclusively use
+     * {@code computeObjectiveValue} to solve the function.
+     * See e.g. {@link AbstractDifferentiableUnivariateRealSolver}.
+     */
+    protected void incrementEvaluationCount() {
+        try {
+            evaluations.incrementCount();
+        } catch (MaxCountExceededException e) {
+            throw new TooManyEvaluationsException(e.getMax());
+        }
+    }
+}

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

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseUnivariateRealSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseUnivariateRealSolver.java?rev=1039083&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseUnivariateRealSolver.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BaseUnivariateRealSolver.java Thu Nov 25 16:22:00 2010
@@ -0,0 +1,114 @@
+/*
+ * 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;
+
+
+/**
+ * Interface for (univariate real) rootfinding algorithms.
+ * Implementations will search for only one zero in the given interval.
+ *
+ * @param <FUNC> Type of function to solve.
+ *
+ * @version $Revision$ $Date$
+ * @since 3.0
+ */
+public interface BaseUnivariateRealSolver<FUNC extends UnivariateRealFunction> {
+    /**
+     * Set the maximal number of function evaluations.
+     *
+     * @param maxEvaluations Maximal number of function evaluations.
+     */
+    void setMaxEvaluations(int maxEvaluations);
+
+    /**
+     * Get the maximal number of function evaluations.
+     *
+     * @return the maximal number of function evaluations.
+     */
+    int getMaxEvaluations();
+
+    /**
+     * Get the number of evaluations of the objective function.
+     * The number of evaluations corresponds to the last call to the
+     * {@code optimize} method. It is 0 if the method has not been
+     * called yet.
+     *
+     * @return the number of evaluations of the objective function.
+     */
+    int getEvaluations();
+
+    /**
+     * @return the absolute accuracy.
+     */
+    double getAbsoluteAccuracy();
+    /**
+     * @return the relative accuracy.
+     */
+    double getRelativeAccuracy();
+    /**
+     * @return the function value accuracy.
+     */
+    double getFunctionValueAccuracy();
+
+    /**
+     * Solve for a zero root in the given interval.
+     * A solver may require that the interval brackets a single zero root.
+     * Solvers that do require bracketing should be able to handle the case
+     * where one of the endpoints is itself a root.
+     *
+     * @param f Function to solve.
+     * @param min Lower bound for the interval.
+     * @param max Upper bound for the interval.
+     * @return a value where the function is zero.
+     * @throws IllegalArgumentException if {@code min > max} or the endpoints
+     * do not satisfy the requirements specified by the solver.
+     * @since 2.0
+     */
+    double solve(FUNC f, double min, double max);
+
+    /**
+     * Solve for a zero in the given interval, start at {@code startValue}.
+     * A solver may require that the interval brackets a single zero root.
+     * Solvers that do require bracketing should be able to handle the case
+     * where one of the endpoints is itself a root.
+     *
+     * @param f Function to solve.
+     * @param min Lower bound for the interval.
+     * @param max Upper bound for the interval.
+     * @param startValue Start value to use.
+     * @return a value where the function is zero.
+     * @throws IllegalArgumentException if {@code min > max} or the arguments
+     * do not satisfy the requirements specified by the solver.
+     * @since 2.0
+     */
+    double solve(FUNC f, double min, double max, double startValue);
+
+    /**
+     * Solve for a zero in the vicinity of {@code startValue}.
+     * A solver may require that the interval brackets a single zero root.
+     *
+     * @param f Function to solve.
+     * @param startValue Start value to use.
+     * @return a value where the function is zero.
+     * @throws IllegalArgumentException if {@code min > max} or the arguments
+     * do not satisfy the requirements specified by the solver.
+     * @since 2.0
+     */
+    double solve(FUNC f, double startValue);
+}

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

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BisectionSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BisectionSolver.java?rev=1039083&r1=1039082&r2=1039083&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BisectionSolver.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BisectionSolver.java Thu Nov 25 16:22:00 2010
@@ -16,9 +16,6 @@
  */
 package org.apache.commons.math.analysis.solvers;
 
-import org.apache.commons.math.MaxIterationsExceededException;
-import org.apache.commons.math.analysis.UnivariateRealFunction;
-import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.util.FastMath;
 
 /**
@@ -29,39 +26,54 @@ import org.apache.commons.math.util.Fast
  *
  * @version $Revision$ $Date$
  */
-public class BisectionSolver extends UnivariateRealSolverImpl {
+public class BisectionSolver extends AbstractUnivariateRealSolver {
+    /** Default absolute accuracy. */
+    public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
 
     /**
+     * Construct a solver with default accuracy.
+     */
+    public BisectionSolver() {
+        this(DEFAULT_ABSOLUTE_ACCURACY);
+    }
+    /**
      * Construct a solver.
      *
+     * @param absoluteAccuracy Absolute accuracy.
      */
-    public BisectionSolver() {
-        super(100, 1E-6);
+    public BisectionSolver(double absoluteAccuracy) {
+        super(absoluteAccuracy);
     }
-
-    /** {@inheritDoc} */
-    public double solve(final UnivariateRealFunction f, double min, double max, double initial)
-        throws MaxIterationsExceededException, MathUserException {
-        return solve(f, min, max);
+    /**
+     * Construct a solver.
+     *
+     * @param relativeAccuracy Relative accuracy.
+     * @param absoluteAccuracy Absolute accuracy.
+     */
+    public BisectionSolver(double relativeAccuracy,
+                           double absoluteAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy);
     }
 
-    /** {@inheritDoc} */
-    public double solve(final UnivariateRealFunction f, double min, double max)
-        throws MaxIterationsExceededException, MathUserException {
-
-        clearResult();
-        verifyInterval(min,max);
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    protected double doSolve() {
+        double min = getMin();
+        double max = getMax();
+        verifyInterval(min, max);
+        final double absoluteAccuracy = getAbsoluteAccuracy();
         double m;
         double fm;
         double fmin;
 
-        int i = 0;
-        while (i < maximalIterationCount) {
+        while (true) {
             m = UnivariateRealSolverUtils.midpoint(min, max);
-           fmin = f.value(min);
-           fm = f.value(m);
+            fmin = computeObjectiveValue(min);
+            fm = computeObjectiveValue(m);
 
-            if (fm * fmin > 0.0) {
+            if (fm * fmin > 0) {
                 // max and m bracket the root.
                 min = m;
             } else {
@@ -71,12 +83,8 @@ public class BisectionSolver extends Uni
 
             if (FastMath.abs(max - min) <= absoluteAccuracy) {
                 m = UnivariateRealSolverUtils.midpoint(min, max);
-                setResult(m, i);
                 return m;
             }
-            ++i;
         }
-
-        throw new MaxIterationsExceededException(maximalIterationCount);
     }
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java?rev=1039083&r1=1039082&r2=1039083&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java Thu Nov 25 16:22:00 2010
@@ -17,297 +17,214 @@
 package org.apache.commons.math.analysis.solvers;
 
 
-import org.apache.commons.math.MathRuntimeException;
-import org.apache.commons.math.MaxIterationsExceededException;
-import org.apache.commons.math.analysis.UnivariateRealFunction;
-import org.apache.commons.math.exception.MathUserException;
-import org.apache.commons.math.exception.util.LocalizedFormats;
+import org.apache.commons.math.exception.NoBracketingException;
 import org.apache.commons.math.util.FastMath;
+import org.apache.commons.math.util.MathUtils;
 
 /**
- * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
- * Brent algorithm</a> for  finding zeros of real univariate functions.
- * <p>
- * The function should be continuous but not necessarily smooth.</p>
+ * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
+ * Brent algorithm</a> for finding zeros of real univariate functions.
+ * The function should be continuous but not necessarily smooth.
+ * The {@code solve} method returns a zero {@code x} of the function {@code f}
+ * in the given interval {@code [a, b]} to within a tolerance
+ * {@code 6 eps abs(x) + t} where {@code eps} is the relative accuracy and
+ * {@code t} is the absolute accuracy.
+ * The given interval must bracket the root.
  *
  * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
  */
-public class BrentSolver extends UnivariateRealSolverImpl {
-
-    /**
-     * Default absolute accuracy
-     * @since 2.1
-     */
-    public static final double DEFAULT_ABSOLUTE_ACCURACY = 1E-6;
-
-    /** Default maximum number of iterations
-     * @since 2.1
-     */
-    public static final int DEFAULT_MAXIMUM_ITERATIONS = 100;
-
+public class BrentSolver extends AbstractUnivariateRealSolver {
     /** Serializable version identifier */
     private static final long serialVersionUID = 7694577816772532779L;
+    /** Default absolute accuracy. */
+    public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
 
     /**
-     * Construct a solver with default properties.
+     * Construct a solver with default accuracies.
      */
     public BrentSolver() {
-        super(DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
+        this(DEFAULT_ABSOLUTE_ACCURACY);
     }
-
     /**
-     * Construct a solver with the given absolute accuracy.
+     * Construct a solver.
      *
-     * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
-     * @since 2.1
+     * @param absoluteAccuracy Absolute accuracy.
      */
     public BrentSolver(double absoluteAccuracy) {
-        super(DEFAULT_MAXIMUM_ITERATIONS, absoluteAccuracy);
+        super(absoluteAccuracy);
     }
-
     /**
-     * Contstruct a solver with the given maximum iterations and absolute accuracy.
+     * Construct a solver.
      *
-     * @param maximumIterations maximum number of iterations
-     * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
-     * @since 2.1
+     * @param relativeAccuracy Relative accuracy.
+     * @param absoluteAccuracy Absolute accuracy.
      */
-    public BrentSolver(int maximumIterations, double absoluteAccuracy) {
-        super(maximumIterations, absoluteAccuracy);
+    public BrentSolver(double relativeAccuracy,
+                       double absoluteAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy);
     }
-
     /**
-     * Find a zero in the given interval with an initial guess.
-     * <p>Throws <code>IllegalArgumentException</code> if the values of the
-     * function at the three points have the same sign (note that it is
-     * allowed to have endpoints with the same sign if the initial point has
-     * opposite sign function-wise).</p>
+     * Construct a solver.
      *
-     * @param f function to solve.
-     * @param min the lower bound for the interval.
-     * @param max the upper bound for the interval.
-     * @param initial the start value to use (must be set to min if no
-     * initial point is known).
-     * @return the value where the function is zero
-     * @throws MaxIterationsExceededException the maximum iteration count is exceeded
-     * @throws MathUserException if an error occurs evaluating  the function
-     * @throws IllegalArgumentException if initial is not between min and max
-     * (even if it <em>is</em> a root)
+     * @param relativeAccuracy Relative accuracy.
+     * @param absoluteAccuracy Absolute accuracy.
+     * @param functionValueAccuracy Function value accuracy.
      */
-    public double solve(final UnivariateRealFunction f,
-                        final double min, final double max, final double initial)
-        throws MaxIterationsExceededException, MathUserException {
-
-        clearResult();
-        if ((initial < min) || (initial > max)) {
-            throw MathRuntimeException.createIllegalArgumentException(
-                  LocalizedFormats.INVALID_INTERVAL_INITIAL_VALUE_PARAMETERS,
-                  min, initial, max);
-        }
+    public BrentSolver(double relativeAccuracy,
+                       double absoluteAccuracy,
+                       double functionValueAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    protected double doSolve() {
+        double min = getMin();
+        double max = getMax();
+        final double initial = getStartValue();
+        final double functionValueAccuracy = getFunctionValueAccuracy();
 
-        // return the initial guess if it is good enough
-        double yInitial = f.value(initial);
+        verifySequence(min, initial, max);
+
+        // Return the initial guess if it is good enough.
+        double yInitial = computeObjectiveValue(initial);
         if (FastMath.abs(yInitial) <= functionValueAccuracy) {
-            setResult(initial, 0);
-            return result;
+            return initial;
         }
 
-        // return the first endpoint if it is good enough
-        double yMin = f.value(min);
+        // Return the first endpoint if it is good enough.
+        double yMin = computeObjectiveValue(min);
         if (FastMath.abs(yMin) <= functionValueAccuracy) {
-            setResult(min, 0);
-            return result;
+            return min;
         }
 
-        // reduce interval if min and initial bracket the root
+        // Reduce interval if min and initial bracket the root.
         if (yInitial * yMin < 0) {
-            return solve(f, min, yMin, initial, yInitial, min, yMin);
+            return brent(min, initial, yMin, yInitial);
         }
 
-        // return the second endpoint if it is good enough
-        double yMax = f.value(max);
+        // Return the second endpoint if it is good enough.
+        double yMax = computeObjectiveValue(max);
         if (FastMath.abs(yMax) <= functionValueAccuracy) {
-            setResult(max, 0);
-            return result;
+            return max;
         }
 
-        // reduce interval if initial and max bracket the root
+        // Reduce interval if initial and max bracket the root.
         if (yInitial * yMax < 0) {
-            return solve(f, initial, yInitial, max, yMax, initial, yInitial);
+            return brent(initial, max, yInitial, yMax);
         }
 
-        throw MathRuntimeException.createIllegalArgumentException(
-              LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, yMin, yMax);
-
+        throw new NoBracketingException(min, max, yMin, yMax);
     }
 
     /**
-     * Find a zero in the given interval.
-     * <p>
-     * Requires that the values of the function at the endpoints have opposite
-     * signs. An <code>IllegalArgumentException</code> is thrown if this is not
-     * the case.</p>
+     * Search for a zero inside the provided interval.
+     * This implemenation is based on the algorithm described at page 58 of
+     * the book
+     * <quote>
+     *  <b>Algorithms for Minimization Without Derivatives</b>
+     *  <it>Richard P. Brent</it>
+     *  Dover 0-486-41998-3
+     * </quote>
      *
-     * @param f the function to solve
-     * @param min the lower bound for the interval.
-     * @param max the upper bound for the interval.
-     * @return the value where the function is zero
-     * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
-     * @throws MathUserException if an error occurs evaluating the function
-     * @throws IllegalArgumentException if min is not less than max or the
-     * signs of the values of the function at the endpoints are not opposites
-     */
-    public double solve(final UnivariateRealFunction f,
-                        final double min, final double max)
-        throws MaxIterationsExceededException, MathUserException {
-
-        clearResult();
-        verifyInterval(min, max);
-
-        double ret = Double.NaN;
-
-        double yMin = f.value(min);
-        double yMax = f.value(max);
-
-        // Verify bracketing
-        double sign = yMin * yMax;
-        if (sign > 0) {
-            // check if either value is close to a zero
-            if (FastMath.abs(yMin) <= functionValueAccuracy) {
-                setResult(min, 0);
-                ret = min;
-            } else if (FastMath.abs(yMax) <= functionValueAccuracy) {
-                setResult(max, 0);
-                ret = max;
-            } else {
-                // neither value is close to zero and min and max do not bracket root.
-                throw MathRuntimeException.createIllegalArgumentException(
-                        LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, yMin, yMax);
-            }
-        } else if (sign < 0){
-            // solve using only the first endpoint as initial guess
-            ret = solve(f, min, yMin, max, yMax, min, yMin);
-        } else {
-            // either min or max is a root
-            if (yMin == 0.0) {
-                ret = min;
-            } else {
-                ret = max;
+     * @param lo Lower bound of the search interval.
+     * @param hi Higher bound of the search interval.
+     * @param fLo Function value at the lower bound of the search interval.
+     * @param fHi Function value at the higher bound of the search interval.
+     * @return the value where the function is zero.
+     */
+    private double brent(double lo, double hi,
+                         double fLo, double fHi) {
+        double a = lo;
+        double fa = fLo;
+        double b = hi;
+        double fb = fHi;
+        double c = a;
+        double fc = fa;
+        double d = b - a;
+        double e = d;
+
+        final double t = getAbsoluteAccuracy();
+        final double eps = getRelativeAccuracy();
+
+        while (true) {
+            if (FastMath.abs(fc) < FastMath.abs(fb)) {
+                a = b;
+                b = c;
+                c = a;
+                fa = fb;
+                fb = fc;
+                fc = fa;
+            }
+
+            final double tol = 2 * eps * FastMath.abs(b) + t;
+            final double m = 0.5 * (c - b);
+
+            if (FastMath.abs(m) <= tol ||
+                MathUtils.equals(fb, 0))  {
+                return b;
             }
-        }
-
-        return ret;
-    }
-
-    /**
-     * Find a zero starting search according to the three provided points.
-     * @param f the function to solve
-     * @param x0 old approximation for the root
-     * @param y0 function value at the approximation for the root
-     * @param x1 last calculated approximation for the root
-     * @param y1 function value at the last calculated approximation
-     * for the root
-     * @param x2 bracket point (must be set to x0 if no bracket point is
-     * known, this will force starting with linear interpolation)
-     * @param y2 function value at the bracket point.
-     * @return the value where the function is zero
-     * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
-     * @throws MathUserException if an error occurs evaluating the function
-     */
-    private double solve(final UnivariateRealFunction f,
-                         double x0, double y0,
-                         double x1, double y1,
-                         double x2, double y2)
-    throws MaxIterationsExceededException, MathUserException {
-
-        double delta = x1 - x0;
-        double oldDelta = delta;
-
-        int i = 0;
-        while (i < maximalIterationCount) {
-            if (FastMath.abs(y2) < FastMath.abs(y1)) {
-                // use the bracket point if is better than last approximation
-                x0 = x1;
-                x1 = x2;
-                x2 = x0;
-                y0 = y1;
-                y1 = y2;
-                y2 = y0;
-            }
-            if (FastMath.abs(y1) <= functionValueAccuracy) {
-                // Avoid division by very small values. Assume
-                // the iteration has converged (the problem may
-                // still be ill conditioned)
-                setResult(x1, i);
-                return result;
-            }
-            double dx = x2 - x1;
-            double tolerance =
-                FastMath.max(relativeAccuracy * FastMath.abs(x1), absoluteAccuracy);
-            if (FastMath.abs(dx) <= tolerance) {
-                setResult(x1, i);
-                return result;
-            }
-            if ((FastMath.abs(oldDelta) < tolerance) ||
-                    (FastMath.abs(y0) <= FastMath.abs(y1))) {
+            if (FastMath.abs(e) < tol ||
+                FastMath.abs(fa) <= FastMath.abs(fb)) {
                 // Force bisection.
-                delta = 0.5 * dx;
-                oldDelta = delta;
+                d = m;
+                e = d;
             } else {
-                double r3 = y1 / y0;
+                double s = fb / fa;
                 double p;
-                double p1;
-                // the equality test (x0 == x2) is intentional,
-                // it is part of the original Brent's method,
-                // it should NOT be replaced by proximity test
-                if (x0 == x2) {
+                double q;
+                // The equality test (a == c) is intentional,
+                // it is part of the original Brent's method and
+                // it should NOT be replaced by proximity test.
+                if (a == c) {
                     // Linear interpolation.
-                    p = dx * r3;
-                    p1 = 1.0 - r3;
+                    p = 2 * m * s;
+                    q = 1 - s;
                 } else {
                     // Inverse quadratic interpolation.
-                    double r1 = y0 / y2;
-                    double r2 = y1 / y2;
-                    p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
-                    p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
+                    q = fa / fc;
+                    final double r = fb / fc;
+                    p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
+                    q = (q - 1) * (r - 1) * (s - 1);
                 }
-                if (p > 0.0) {
-                    p1 = -p1;
+                if (p > 0) {
+                    q = -q;
                 } else {
                     p = -p;
                 }
-                if (2.0 * p >= 1.5 * dx * p1 - FastMath.abs(tolerance * p1) ||
-                        p >= FastMath.abs(0.5 * oldDelta * p1)) {
+                s = e;
+                e = d;
+                if (p >= 1.5 * m * q - FastMath.abs(tol * q) ||
+                    p >= FastMath.abs(0.5 * s * q)) {
                     // Inverse quadratic interpolation gives a value
                     // in the wrong direction, or progress is slow.
                     // Fall back to bisection.
-                    delta = 0.5 * dx;
-                    oldDelta = delta;
+                    d = m;
+                    e = d;
                 } else {
-                    oldDelta = delta;
-                    delta = p / p1;
+                    d = p / q;
                 }
             }
-            // Save old X1, Y1
-            x0 = x1;
-            y0 = y1;
-            // Compute new X1, Y1
-            if (FastMath.abs(delta) > tolerance) {
-                x1 = x1 + delta;
-            } else if (dx > 0.0) {
-                x1 = x1 + 0.5 * tolerance;
-            } else if (dx <= 0.0) {
-                x1 = x1 - 0.5 * tolerance;
+            a = b;
+            fa = fb;
+
+            if (FastMath.abs(d) > tol) {
+                b += d;
+            } else if (m > 0) {
+                b += tol;
+            } else {
+                b -= tol;
             }
-            y1 = f.value(x1);
-            if ((y1 > 0) == (y2 > 0)) {
-                x2 = x0;
-                y2 = y0;
-                delta = x1 - x0;
-                oldDelta = delta;
+            fb = computeObjectiveValue(b);
+            if ((fb > 0 && fc > 0) ||
+                (fb <= 0 && fc <= 0)) {
+                c = a;
+                fc = fa;
+                d = b - a;
+                e = d;
             }
-            i++;
         }
-        throw new MaxIterationsExceededException(maximalIterationCount);
     }
 }

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/DifferentiableUnivariateRealSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/DifferentiableUnivariateRealSolver.java?rev=1039083&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/DifferentiableUnivariateRealSolver.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/DifferentiableUnivariateRealSolver.java Thu Nov 25 16:22:00 2010
@@ -0,0 +1,29 @@
+/*
+ * 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;
+
+
+/**
+ * Interface for (univariate real) rootfinding algorithms.
+ * Implementations will search for only one zero in the given interval.
+ *
+ * @version $Revision: 1034896 $ $Date: 2010-11-13 23:27:34 +0100 (Sat, 13 Nov 2010) $
+ */
+public interface DifferentiableUnivariateRealSolver
+    extends BaseUnivariateRealSolver<DifferentiableUnivariateRealFunction> {}

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

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java?rev=1039083&r1=1039082&r2=1039083&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java Thu Nov 25 16:22:00 2010
@@ -16,304 +16,335 @@
  */
 package org.apache.commons.math.analysis.solvers;
 
-import org.apache.commons.math.ConvergenceException;
-import org.apache.commons.math.MathRuntimeException;
-import org.apache.commons.math.MaxIterationsExceededException;
-import org.apache.commons.math.analysis.UnivariateRealFunction;
-import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
 import org.apache.commons.math.complex.Complex;
-import org.apache.commons.math.exception.MathUserException;
+import org.apache.commons.math.exception.NoBracketingException;
+import org.apache.commons.math.exception.NullArgumentException;
+import org.apache.commons.math.exception.NoDataException;
 import org.apache.commons.math.exception.util.LocalizedFormats;
 import org.apache.commons.math.util.FastMath;
 
 /**
  * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
  * Laguerre's Method</a> for root finding of real coefficient polynomials.
- * For reference, see <b>A First Course in Numerical Analysis</b>,
- * ISBN 048641454X, chapter 8.
- * <p>
+ * For reference, see
+ * <quote>
+ *  <b>A First Course in Numerical Analysis</b>
+ *  ISBN 048641454X, chapter 8.
+ * </quote>
  * Laguerre's method is global in the sense that it can start with any initial
- * approximation and be able to solve all roots from that point.</p>
+ * approximation and be able to solve all roots from that point.
+ * The algorithm requires a bracketing condition.
  *
  * @version $Revision$ $Date$
  * @since 1.2
  */
-public class LaguerreSolver extends UnivariateRealSolverImpl {
+public class LaguerreSolver extends AbstractPolynomialSolver {
+    /** Default absolute accuracy. */
+    public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
+    /** Complex solver. */
+    protected ComplexSolver complexSolver = new ComplexSolver();
 
     /**
-     * Construct a solver.
+     * Construct a solver with default accuracies.
      */
     public LaguerreSolver() {
-        super(100, 1E-6);
+        this(DEFAULT_ABSOLUTE_ACCURACY);
     }
-
     /**
-     * Find a real root in the given interval with initial value.
-     * <p>
-     * Requires bracketing condition.</p>
+     * Construct a solver.
      *
-     * @param f function to solve (must be polynomial)
-     * @param min the lower bound for the interval
-     * @param max the upper bound for the interval
-     * @param initial the start value to use
-     * @return the point at which the function value is zero
-     * @throws ConvergenceException if the maximum iteration count is exceeded
-     * or the solver detects convergence problems otherwise
-     * @throws MathUserException if an error occurs evaluating the function
-     * @throws IllegalArgumentException if any parameters are invalid
+     * @param absoluteAccuracy Absolute accuracy.
      */
-    public double solve(final UnivariateRealFunction f,
-                        final double min, final double max, final double initial)
-        throws ConvergenceException, MathUserException {
-
-        // check for zeros before verifying bracketing
-        if (f.value(min) == 0.0) {
-            return min;
-        }
-        if (f.value(max) == 0.0) {
-            return max;
-        }
-        if (f.value(initial) == 0.0) {
-            return initial;
-        }
-
-        verifyBracketing(min, max, f);
-        verifySequence(min, initial, max);
-        if (isBracketing(min, initial, f)) {
-            return solve(f, min, initial);
-        } else {
-            return solve(f, initial, max);
-        }
-
+    public LaguerreSolver(double absoluteAccuracy) {
+        super(absoluteAccuracy);
     }
-
     /**
-     * Find a real root in the given interval.
-     * <p>
-     * Despite the bracketing condition, the root returned by solve(Complex[],
-     * Complex) may not be a real zero inside [min, max]. For example,
-     * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
-     * another initial value, or, as we did here, call solveAll() to obtain
-     * all roots and pick up the one that we're looking for.</p>
+     * Construct a solver.
+     *
+     * @param relativeAccuracy Relative accuracy.
+     * @param absoluteAccuracy Absolute accuracy.
+     */
+    public LaguerreSolver(double relativeAccuracy,
+                          double absoluteAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy);
+    }
+    /**
+     * Construct a solver.
      *
-     * @param f the function to solve
-     * @param min the lower bound for the interval
-     * @param max the upper bound for the interval
-     * @return the point at which the function value is zero
-     * @throws ConvergenceException if the maximum iteration count is exceeded
-     * or the solver detects convergence problems otherwise
-     * @throws MathUserException if an error occurs evaluating the function
-     * @throws IllegalArgumentException if any parameters are invalid
+     * @param relativeAccuracy Relative accuracy.
+     * @param absoluteAccuracy Absolute accuracy.
+     * @param functionValueAccuracy Function value accuracy.
      */
-    public double solve(final UnivariateRealFunction f,
-                        final double min, final double max)
-        throws ConvergenceException, MathUserException {
+    public LaguerreSolver(double relativeAccuracy,
+                          double absoluteAccuracy,
+                          double functionValueAccuracy) {
+        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
+    }
 
-        // check function type
-        if (!(f instanceof PolynomialFunction)) {
-            throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL);
+    /**
+     * {@inheritDoc}
+     */
+    @Override
+    public double doSolve() {
+        double min = getMin();
+        double max = getMax();
+        double initial = getStartValue();
+        final double functionValueAccuracy = getFunctionValueAccuracy();
+
+        verifySequence(min, initial, max);
+
+        // Return the initial guess if it is good enough.
+        double yInitial = computeObjectiveValue(initial);
+        if (FastMath.abs(yInitial) <= functionValueAccuracy) {
+            return initial;
         }
 
-        // check for zeros before verifying bracketing
-        if (f.value(min) == 0.0) { return min; }
-        if (f.value(max) == 0.0) { return max; }
-        verifyBracketing(min, max, f);
+        // Return the first endpoint if it is good enough.
+        double yMin = computeObjectiveValue(min);
+        if (FastMath.abs(yMin) <= functionValueAccuracy) {
+            return min;
+        }
 
-        double coefficients[] = ((PolynomialFunction) f).getCoefficients();
-        Complex c[] = new Complex[coefficients.length];
-        for (int i = 0; i < coefficients.length; i++) {
-            c[i] = new Complex(coefficients[i], 0.0);
+        // Reduce interval if min and initial bracket the root.
+        if (yInitial * yMin < 0) {
+            return laguerre(min, initial, yMin, yInitial);
         }
-        Complex initial = new Complex(0.5 * (min + max), 0.0);
-        Complex z = solve(c, initial);
-        if (isRootOK(min, max, z)) {
-            setResult(z.getReal(), iterationCount);
-            return result;
+
+        // Return the second endpoint if it is good enough.
+        double yMax = computeObjectiveValue(max);
+        if (FastMath.abs(yMax) <= functionValueAccuracy) {
+            return max;
         }
 
-        // solve all roots and select the one we're seeking
-        Complex[] root = solveAll(c, initial);
-        for (int i = 0; i < root.length; i++) {
-            if (isRootOK(min, max, root[i])) {
-                setResult(root[i].getReal(), iterationCount);
-                return result;
-            }
+        // Reduce interval if initial and max bracket the root.
+        if (yInitial * yMax < 0) {
+            return laguerre(initial, max, yInitial, yMax);
         }
 
-        // should never happen
-        throw new ConvergenceException();
+        throw new NoBracketingException(min, max, yMin, yMax);
     }
 
     /**
-     * Returns true iff the given complex root is actually a real zero
-     * in the given interval, within the solver tolerance level.
+     * Find a real root in the given interval.
      *
-     * @param min the lower bound for the interval
-     * @param max the upper bound for the interval
-     * @param z the complex root
-     * @return true iff z is the sought-after real zero
-     */
-    protected boolean isRootOK(double min, double max, Complex z) {
-        double tolerance = FastMath.max(relativeAccuracy * z.abs(), absoluteAccuracy);
-        return (isSequence(min, z.getReal(), max)) &&
-               (FastMath.abs(z.getImaginary()) <= tolerance ||
-                z.abs() <= functionValueAccuracy);
-    }
-
-    /**
-     * Find all complex roots for the polynomial with the given coefficients,
-     * starting from the given initial value.
+     * Despite the bracketing condition, the root returned by
+     * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may
+     * not be a real zero inside {@code [min, max]}.
+     * For example, <code>p(x) = x<sup>3</sup> + 1,</code>
+     * with {@code min = -2}, {@code max = 2}, {@code initial = 0}.
+     * When it occurs, this code calls
+     * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)}
+     * in order to obtain all roots and picks up one real root.
      *
-     * @param coefficients the polynomial coefficients array
-     * @param initial the start value to use
-     * @return the point at which the function value is zero
-     * @throws ConvergenceException if the maximum iteration count is exceeded
-     * or the solver detects convergence problems otherwise
-     * @throws MathUserException if an error occurs evaluating the function
-     * @throws IllegalArgumentException if any parameters are invalid
+     * @param lo Lower bound of the search interval.
+     * @param hi Higher bound of the search interval.
+     * @param fLo Function value at the lower bound of the search interval.
+     * @param fHi Function value at the higher bound of the search interval.
+     * @return the point at which the function value is zero.
      */
-    public Complex[] solveAll(double coefficients[], double initial) throws
-        ConvergenceException, MathUserException {
-
+    public double laguerre(double lo, double hi,
+                           double fLo, double fHi) {
+        double result = Double.NaN;
+        double coefficients[] = getCoefficients();
         Complex c[] = new Complex[coefficients.length];
-        Complex z = new Complex(initial, 0.0);
-        for (int i = 0; i < c.length; i++) {
-            c[i] = new Complex(coefficients[i], 0.0);
+        for (int i = 0; i < coefficients.length; i++) {
+            c[i] = new Complex(coefficients[i], 0);
+        }
+        Complex initial = new Complex(0.5 * (lo + hi), 0);
+        Complex z = complexSolver.solve(c, initial);
+        if (complexSolver.isRoot(lo, hi, z)) {
+            return z.getReal();
+        } else {
+            double r = Double.NaN;
+            // Solve all roots and select the one we are seeking.
+            Complex[] root = complexSolver.solveAll(c, initial);
+            for (int i = 0; i < root.length; i++) {
+                if (complexSolver.isRoot(lo, hi, root[i])) {
+                    r = root[i].getReal();
+                    break;
+                }
+            }
+            return r;
         }
-        return solveAll(c, z);
     }
 
     /**
-     * Find all complex roots for the polynomial with the given coefficients,
-     * starting from the given initial value.
-     *
-     * @param coefficients the polynomial coefficients array
-     * @param initial the start value to use
-     * @return the point at which the function value is zero
-     * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
-     * or the solver detects convergence problems otherwise
-     * @throws MathUserException if an error occurs evaluating the function
-     * @throws IllegalArgumentException if any parameters are invalid
+     * Class for searching all (complex) roots.
      */
-    public Complex[] solveAll(Complex coefficients[], Complex initial) throws
-        MaxIterationsExceededException, MathUserException {
-
-        int n = coefficients.length - 1;
-        int iterationCount = 0;
-        if (n < 1) {
-            throw MathRuntimeException.createIllegalArgumentException(
-                  LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n);
-        }
-        Complex c[] = new Complex[n+1];    // coefficients for deflated polynomial
-        for (int i = 0; i <= n; i++) {
-            c[i] = coefficients[i];
-        }
-
-        // solve individual root successively
-        Complex root[] = new Complex[n];
-        for (int i = 0; i < n; i++) {
-            Complex subarray[] = new Complex[n-i+1];
-            System.arraycopy(c, 0, subarray, 0, subarray.length);
-            root[i] = solve(subarray, initial);
-            // polynomial deflation using synthetic division
-            Complex newc = c[n-i];
-            Complex oldc = null;
-            for (int j = n-i-1; j >= 0; j--) {
-                oldc = c[j];
-                c[j] = newc;
-                newc = oldc.add(newc.multiply(root[i]));
+    private class ComplexSolver {
+        /**
+         * Check whether the given complex root is actually a real zero
+         * in the given interval, within the solver tolerance level.
+         *
+         * @param min Lower bound for the interval.
+         * @param max Upper bound for the interval.
+         * @param z Complex root.
+         * @return {@code true} if z is a real zero.
+         */
+        public boolean isRoot(double min, double max, Complex z) {
+            double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
+            return (isSequence(min, z.getReal(), max)) &&
+                (FastMath.abs(z.getImaginary()) <= tolerance ||
+                 z.abs() <= getFunctionValueAccuracy());
+        }
+
+        /**
+         * Find all complex roots for the polynomial with the given
+         * coefficients, starting from the given initial value.
+         *
+         * @param coefficients Polynomial coefficients.
+         * @param initial Start value.
+         * @return the point at which the function value is zero.
+         * @throws org.apache.commons.math.exception.TooManyEvaluationsException
+         * if the maximum number of evaluations is exceeded.
+         * @throws NullArgumentException if the {@code coefficients} is
+         * {@code null}.
+         * @throws NoDataException if the {@code coefficients} array is empty.
+         */
+        public Complex[] solveAll(double coefficients[], double initial) {
+            if (coefficients == null) {
+                throw new NullArgumentException();
+            }
+            Complex c[] = new Complex[coefficients.length];
+            Complex z = new Complex(initial, 0);
+            for (int i = 0; i < c.length; i++) {
+                c[i] = new Complex(coefficients[i], 0);
             }
-            iterationCount += this.iterationCount;
+            return solveAll(c, z);
         }
 
-        resultComputed = true;
-        this.iterationCount = iterationCount;
-        return root;
-    }
-
-    /**
-     * Find a complex root for the polynomial with the given coefficients,
-     * starting from the given initial value.
-     *
-     * @param coefficients the polynomial coefficients array
-     * @param initial the start value to use
-     * @return the point at which the function value is zero
-     * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
-     * or the solver detects convergence problems otherwise
-     * @throws MathUserException if an error occurs evaluating the function
-     * @throws IllegalArgumentException if any parameters are invalid
-     */
-    public Complex solve(Complex coefficients[], Complex initial) throws
-        MaxIterationsExceededException, MathUserException {
+        /**
+         * Find all complex roots for the polynomial with the given
+         * coefficients, starting from the given initial value.
+         *
+         * @param coefficients Polynomial coefficients.
+         * @param initial Start value.
+         * @return the point at which the function value is zero.
+         * @throws org.apache.commons.math.exception.TooManyEvaluationsException
+         * if the maximum number of evaluations is exceeded.
+         * @throws NullArgumentException if the {@code coefficients} is
+         * {@code null}.
+         * @throws NoDataException if the {@code coefficients} array is empty.
+         */
+        public Complex[] solveAll(Complex coefficients[], Complex initial) {
+            if (coefficients == null) {
+                throw new NullArgumentException();
+            }
+            int n = coefficients.length - 1;
+            if (n == 0) {
+                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
+            }
+            // Coefficients for deflated polynomial.
+            Complex c[] = new Complex[n + 1];
+            for (int i = 0; i <= n; i++) {
+                c[i] = coefficients[i];
+            }
 
-        int n = coefficients.length - 1;
-        if (n < 1) {
-            throw MathRuntimeException.createIllegalArgumentException(
-                  LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n);
-        }
-        Complex N  = new Complex(n,     0.0);
-        Complex N1 = new Complex(n - 1, 0.0);
-
-        int i = 1;
-        Complex pv = null;
-        Complex dv = null;
-        Complex d2v = null;
-        Complex G = null;
-        Complex G2 = null;
-        Complex H = null;
-        Complex delta = null;
-        Complex denominator = null;
-        Complex z = initial;
-        Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
-        while (i <= maximalIterationCount) {
-            // Compute pv (polynomial value), dv (derivative value), and
-            // d2v (second derivative value) simultaneously.
-            pv = coefficients[n];
-            dv = Complex.ZERO;
-            d2v = Complex.ZERO;
-            for (int j = n-1; j >= 0; j--) {
-                d2v = dv.add(z.multiply(d2v));
-                dv = pv.add(z.multiply(dv));
-                pv = coefficients[j].add(z.multiply(pv));
+            // Solve individual roots successively.
+            Complex root[] = new Complex[n];
+            for (int i = 0; i < n; i++) {
+                Complex subarray[] = new Complex[n - i + 1];
+                System.arraycopy(c, 0, subarray, 0, subarray.length);
+                root[i] = solve(subarray, initial);
+                // Polynomial deflation using synthetic division.
+                Complex newc = c[n - i];
+                Complex oldc = null;
+                for (int j = n - i - 1; j >= 0; j--) {
+                    oldc = c[j];
+                    c[j] = newc;
+                    newc = oldc.add(newc.multiply(root[i]));
+                }
             }
-            d2v = d2v.multiply(new Complex(2.0, 0.0));
 
-            // check for convergence
-            double tolerance = FastMath.max(relativeAccuracy * z.abs(),
-                                        absoluteAccuracy);
-            if ((z.subtract(oldz)).abs() <= tolerance) {
-                resultComputed = true;
-                iterationCount = i;
-                return z;
+            return root;
+        }
+
+        /**
+         * Find a complex root for the polynomial with the given coefficients,
+         * starting from the given initial value.
+         *
+         * @param coefficients Polynomial coefficients.
+         * @param initial Start value.
+         * @return the point at which the function value is zero.
+         * @throws org.apache.commons.math.exception.TooManyEvaluationsException
+         * if the maximum number of evaluations is exceeded.
+         * @throws NullArgumentException if the {@code coefficients} is
+         * {@code null}.
+         * @throws NoDataException if the {@code coefficients} array is empty.
+         */
+        public Complex solve(Complex coefficients[], Complex initial) {
+            if (coefficients == null) {
+                throw new NullArgumentException();
             }
-            if (pv.abs() <= functionValueAccuracy) {
-                resultComputed = true;
-                iterationCount = i;
-                return z;
+
+            int n = coefficients.length - 1;
+            if (n == 0) {
+                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
             }
 
-            // now pv != 0, calculate the new approximation
-            G = dv.divide(pv);
-            G2 = G.multiply(G);
-            H = G2.subtract(d2v.divide(pv));
-            delta = N1.multiply((N.multiply(H)).subtract(G2));
-            // choose a denominator larger in magnitude
-            Complex deltaSqrt = delta.sqrt();
-            Complex dplus = G.add(deltaSqrt);
-            Complex dminus = G.subtract(deltaSqrt);
-            denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
-            // Perturb z if denominator is zero, for instance,
-            // p(x) = x^3 + 1, z = 0.
-            if (denominator.equals(new Complex(0.0, 0.0))) {
-                z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
-                oldz = new Complex(Double.POSITIVE_INFINITY,
-                                   Double.POSITIVE_INFINITY);
-            } else {
-                oldz = z;
-                z = z.subtract(N.divide(denominator));
+            final double absoluteAccuracy = getAbsoluteAccuracy();
+            final double relativeAccuracy = getRelativeAccuracy();
+            final double functionValueAccuracy = getFunctionValueAccuracy();
+
+            Complex N  = new Complex(n,     0.0);
+            Complex N1 = new Complex(n - 1, 0.0);
+
+            Complex pv = null;
+            Complex dv = null;
+            Complex d2v = null;
+            Complex G = null;
+            Complex G2 = null;
+            Complex H = null;
+            Complex delta = null;
+            Complex denominator = null;
+            Complex z = initial;
+            Complex oldz = new Complex(Double.POSITIVE_INFINITY,
+                                       Double.POSITIVE_INFINITY);
+            while (true) {
+                // Compute pv (polynomial value), dv (derivative value), and
+                // d2v (second derivative value) simultaneously.
+                pv = coefficients[n];
+                dv = Complex.ZERO;
+                d2v = Complex.ZERO;
+                for (int j = n-1; j >= 0; j--) {
+                    d2v = dv.add(z.multiply(d2v));
+                    dv = pv.add(z.multiply(dv));
+                    pv = coefficients[j].add(z.multiply(pv));
+                }
+                d2v = d2v.multiply(new Complex(2.0, 0.0));
+
+                // check for convergence
+                double tolerance = FastMath.max(relativeAccuracy * z.abs(),
+                                                absoluteAccuracy);
+                if ((z.subtract(oldz)).abs() <= tolerance) {
+                    return z;
+                }
+                if (pv.abs() <= functionValueAccuracy) {
+                    return z;
+                }
+
+                // now pv != 0, calculate the new approximation
+                G = dv.divide(pv);
+                G2 = G.multiply(G);
+                H = G2.subtract(d2v.divide(pv));
+                delta = N1.multiply((N.multiply(H)).subtract(G2));
+                // choose a denominator larger in magnitude
+                Complex deltaSqrt = delta.sqrt();
+                Complex dplus = G.add(deltaSqrt);
+                Complex dminus = G.subtract(deltaSqrt);
+                denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
+                // Perturb z if denominator is zero, for instance,
+                // p(x) = x^3 + 1, z = 0.
+                if (denominator.equals(new Complex(0.0, 0.0))) {
+                    z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
+                    oldz = new Complex(Double.POSITIVE_INFINITY,
+                                       Double.POSITIVE_INFINITY);
+                } else {
+                    oldz = z;
+                    z = z.subtract(N.divide(denominator));
+                }
+                incrementEvaluationCount();
             }
-            i++;
         }
-        throw new MaxIterationsExceededException(maximalIterationCount);
     }
 }



Mime
View raw message