commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1182134 [1/3] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/analysis/function/ main/java/org/apache/commons/math/analysis/interpolation/ main/java/org/apache/commons/math/analysis/polynomials/ main/java/org/apache/...
Date Tue, 11 Oct 2011 22:55:09 GMT
Author: erans
Date: Tue Oct 11 22:55:08 2011
New Revision: 1182134

URL: http://svn.apache.org/viewvc?rev=1182134&view=rev
Log:
MATH-689
Moved arrays utilities from "MathUtils" to "MathArrays".

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathArrays.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/MathArraysTest.java   (with props)
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/function/StepFunction.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LinearInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SplineInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/UnivariateRealPeriodicInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialSplineFunction.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/NonMonotonicSequenceException.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/euclidean/threed/Vector3D.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/PowellOptimizer.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/clustering/EuclideanIntegerPoint.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/MillerUpdatingRegression.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/RegressionResults.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MultidimensionalCounter.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/Precision.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/exception/NonMonotonicSequenceExceptionTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/MathUtilsTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/function/StepFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/function/StepFunction.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/function/StepFunction.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/function/StepFunction.java Tue Oct 11 22:55:08 2011
@@ -22,7 +22,7 @@ import org.apache.commons.math.analysis.
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.NullArgumentException;
 import org.apache.commons.math.exception.NoDataException;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * <a href="http://en.wikipedia.org/wiki/Step_function">
@@ -68,10 +68,10 @@ public class StepFunction implements Uni
         if (y.length != x.length) {
             throw new DimensionMismatchException(y.length, x.length);
         }
-        MathUtils.checkOrder(x);
+        MathArrays.checkOrder(x);
 
-        abscissa = MathUtils.copyOf(x);
-        ordinate = MathUtils.copyOf(y);
+        abscissa = MathArrays.copyOf(x);
+        ordinate = MathArrays.copyOf(y);
     }
 
     /** {@inheritDoc} */

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java Tue Oct 11 22:55:08 2011
@@ -20,7 +20,7 @@ import org.apache.commons.math.analysis.
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.NoDataException;
 import org.apache.commons.math.exception.OutOfRangeException;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Function that implements the
@@ -114,8 +114,8 @@ public class BicubicSplineInterpolatingF
             throw new DimensionMismatchException(xLen, d2FdXdY.length);
         }
 
-        MathUtils.checkOrder(x);
-        MathUtils.checkOrder(y);
+        MathArrays.checkOrder(x);
+        MathArrays.checkOrder(y);
 
         xval = x.clone();
         yval = y.clone();

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolator.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolator.java Tue Oct 11 22:55:08 2011
@@ -20,7 +20,7 @@ import org.apache.commons.math.analysis.
 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.NoDataException;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Generates a bicubic interpolating function.
@@ -43,8 +43,8 @@ public class BicubicSplineInterpolator
             throw new DimensionMismatchException(xval.length, fval.length);
         }
 
-        MathUtils.checkOrder(xval);
-        MathUtils.checkOrder(yval);
+        MathArrays.checkOrder(xval);
+        MathArrays.checkOrder(yval);
 
         final int xLen = xval.length;
         final int yLen = yval.length;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LinearInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LinearInterpolator.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LinearInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LinearInterpolator.java Tue Oct 11 22:55:08 2011
@@ -21,7 +21,7 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.exception.NumberIsTooSmallException;
 import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Implements a linear function for interpolation of real univariate functions.
@@ -53,7 +53,7 @@ public class LinearInterpolator implemen
         // Number of intervals.  The number of data points is n + 1.
         int n = x.length - 1;
 
-        MathUtils.checkOrder(x);
+        MathArrays.checkOrder(x);
 
         // Slope of the lines between the datapoints.
         final double m[] = new double[n];

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java Tue Oct 11 22:55:08 2011
@@ -28,6 +28,7 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.exception.util.LocalizedFormats;
 import org.apache.commons.math.util.FastMath;
 import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Implements the <a href="http://en.wikipedia.org/wiki/Local_regression">
@@ -216,7 +217,7 @@ public class LoessInterpolator
         checkAllFiniteReal(yval);
         checkAllFiniteReal(weights);
 
-        MathUtils.checkOrder(xval);
+        MathArrays.checkOrder(xval);
 
         if (n == 1) {
             return new double[]{yval[0]};

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java Tue Oct 11 22:55:08 2011
@@ -18,7 +18,7 @@ package org.apache.commons.math.analysis
 
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.NoDataException;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
 import org.apache.commons.math.optimization.fitting.PolynomialFitter;
 import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
@@ -87,8 +87,8 @@ public class SmoothingPolynomialBicubicS
             }
         }
 
-        MathUtils.checkOrder(xval);
-        MathUtils.checkOrder(yval);
+        MathArrays.checkOrder(xval);
+        MathArrays.checkOrder(yval);
 
         // For each line y[j] (0 <= j < yLen), construct a polynomial, with
         // respect to variable x, fitting array fval[][j]

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SplineInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SplineInterpolator.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SplineInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SplineInterpolator.java Tue Oct 11 22:55:08 2011
@@ -21,7 +21,7 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.exception.NumberIsTooSmallException;
 import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Computes a natural (also known as "free", "unclamped") cubic spline interpolation for the data set.
@@ -77,7 +77,7 @@ public class SplineInterpolator implemen
         // Number of intervals.  The number of data points is n + 1.
         int n = x.length - 1;
 
-        MathUtils.checkOrder(x);
+        MathArrays.checkOrder(x);
 
         // Differences between knot points
         double h[] = new double[n];

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java Tue Oct 11 22:55:08 2011
@@ -20,7 +20,7 @@ import org.apache.commons.math.analysis.
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.NoDataException;
 import org.apache.commons.math.exception.OutOfRangeException;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Function that implements the
@@ -185,9 +185,9 @@ public class TricubicSplineInterpolating
             throw new DimensionMismatchException(xLen, d3FdXdYdZ.length);
         }
 
-        MathUtils.checkOrder(x);
-        MathUtils.checkOrder(y);
-        MathUtils.checkOrder(z);
+        MathArrays.checkOrder(x);
+        MathArrays.checkOrder(y);
+        MathArrays.checkOrder(z);
 
         xval = x.clone();
         yval = y.clone();

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java Tue Oct 11 22:55:08 2011
@@ -18,7 +18,7 @@ package org.apache.commons.math.analysis
 
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.NoDataException;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Generates a tricubic interpolating function.
@@ -42,9 +42,9 @@ public class TricubicSplineInterpolator
             throw new DimensionMismatchException(xval.length, fval.length);
         }
 
-        MathUtils.checkOrder(xval);
-        MathUtils.checkOrder(yval);
-        MathUtils.checkOrder(zval);
+        MathArrays.checkOrder(xval);
+        MathArrays.checkOrder(yval);
+        MathArrays.checkOrder(zval);
 
         final int xLen = xval.length;
         final int yLen = yval.length;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/UnivariateRealPeriodicInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/UnivariateRealPeriodicInterpolator.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/UnivariateRealPeriodicInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/UnivariateRealPeriodicInterpolator.java Tue Oct 11 22:55:08 2011
@@ -18,6 +18,7 @@ package org.apache.commons.math.analysis
 
 import org.apache.commons.math.analysis.UnivariateRealFunction;
 import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 import org.apache.commons.math.exception.NumberIsTooSmallException;
 
 /**
@@ -85,7 +86,7 @@ public class UnivariateRealPeriodicInter
             throw new NumberIsTooSmallException(xval.length, extend, true);
         }
 
-        MathUtils.checkOrder(xval);
+        MathArrays.checkOrder(xval);
         final double offset = xval[0];
 
         final int len = xval.length + extend * 2;
@@ -108,7 +109,7 @@ public class UnivariateRealPeriodicInter
             y[index] = yval[i];
         }
 
-        MathUtils.sortInPlace(x, y);
+        MathArrays.sortInPlace(x, y);
 
         final UnivariateRealFunction f = interpolator.interpolate(x, y);
         return new UnivariateRealFunction() {

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.java Tue Oct 11 22:55:08 2011
@@ -18,7 +18,7 @@ package org.apache.commons.math.analysis
 
 import org.apache.commons.math.analysis.UnivariateRealFunction;
 import org.apache.commons.math.util.FastMath;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.NumberIsTooSmallException;
 import org.apache.commons.math.exception.util.LocalizedFormats;
@@ -76,7 +76,7 @@ public class PolynomialFunctionLagrangeF
         coefficientsComputed = false;
 
         if (!verifyInterpolationArray(x, y, false)) {
-            MathUtils.sortInPlace(this.x, this.y);
+            MathArrays.sortInPlace(this.x, this.y);
             // Second check in case some abscissa is duplicated.
             verifyInterpolationArray(this.x, this.y, true);
         }
@@ -179,7 +179,7 @@ public class PolynomialFunctionLagrangeF
         System.arraycopy(x, 0, xNew, 0, x.length);
         System.arraycopy(y, 0, yNew, 0, y.length);
 
-        MathUtils.sortInPlace(xNew, yNew);
+        MathArrays.sortInPlace(xNew, yNew);
         // Second check in case some abscissa is duplicated.
         verifyInterpolationArray(xNew, yNew, true);
         return evaluateInternal(xNew, yNew, z);
@@ -318,6 +318,6 @@ public class PolynomialFunctionLagrangeF
             throw new NumberIsTooSmallException(LocalizedFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true);
         }
 
-        return MathUtils.checkOrder(x, MathUtils.OrderDirection.INCREASING, true, abort);
+        return MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, abort);
     }
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialSplineFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialSplineFunction.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialSplineFunction.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialSplineFunction.java Tue Oct 11 22:55:08 2011
@@ -18,7 +18,7 @@ package org.apache.commons.math.analysis
 
 import java.util.Arrays;
 
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
 import org.apache.commons.math.analysis.UnivariateRealFunction;
 import org.apache.commons.math.exception.OutOfRangeException;
@@ -109,7 +109,7 @@ public class PolynomialSplineFunction im
         if (knots.length - 1 != polynomials.length) {
             throw new DimensionMismatchException(polynomials.length, knots.length);
         }
-        MathUtils.checkOrder(knots);
+        MathArrays.checkOrder(knots);
 
         this.n = knots.length -1;
         this.knots = new double[n + 1];

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/NonMonotonicSequenceException.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/NonMonotonicSequenceException.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/NonMonotonicSequenceException.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/NonMonotonicSequenceException.java Tue Oct 11 22:55:08 2011
@@ -16,7 +16,7 @@
  */
 package org.apache.commons.math.exception;
 
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 import org.apache.commons.math.exception.util.LocalizedFormats;
 
 /**
@@ -32,7 +32,7 @@ public class NonMonotonicSequenceExcepti
     /**
      * Direction (positive for increasing, negative for decreasing).
      */
-    private final MathUtils.OrderDirection direction;
+    private final MathArrays.OrderDirection direction;
     /**
      * Whether the sequence must be strictly increasing or decreasing.
      */
@@ -58,7 +58,7 @@ public class NonMonotonicSequenceExcepti
     public NonMonotonicSequenceException(Number wrong,
                                          Number previous,
                                          int index) {
-        this(wrong, previous, index, MathUtils.OrderDirection.INCREASING, true);
+        this(wrong, previous, index, MathArrays.OrderDirection.INCREASING, true);
     }
 
     /**
@@ -75,9 +75,9 @@ public class NonMonotonicSequenceExcepti
     public NonMonotonicSequenceException(Number wrong,
                                          Number previous,
                                          int index,
-                                         MathUtils.OrderDirection direction,
+                                         MathArrays.OrderDirection direction,
                                          boolean strict) {
-        super(direction == MathUtils.OrderDirection.INCREASING ?
+        super(direction == MathArrays.OrderDirection.INCREASING ?
               (strict ?
                LocalizedFormats.NOT_STRICTLY_INCREASING_SEQUENCE :
                LocalizedFormats.NOT_INCREASING_SEQUENCE) :
@@ -95,7 +95,7 @@ public class NonMonotonicSequenceExcepti
     /**
      * @return the order direction.
      **/
-    public MathUtils.OrderDirection getDirection() {
+    public MathArrays.OrderDirection getDirection() {
         return direction;
     }
     /**

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/euclidean/threed/Vector3D.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/euclidean/threed/Vector3D.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/euclidean/threed/Vector3D.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/geometry/euclidean/threed/Vector3D.java Tue Oct 11 22:55:08 2011
@@ -26,6 +26,7 @@ import org.apache.commons.math.geometry.
 import org.apache.commons.math.geometry.Space;
 import org.apache.commons.math.util.FastMath;
 import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * This class implements vectors in a three-dimensional space.
@@ -132,9 +133,9 @@ public class Vector3D implements Seriali
      * @param u2 second base (unscaled) vector
      */
     public Vector3D(double a1, Vector3D u1, double a2, Vector3D u2) {
-        this.x = MathUtils.linearCombination(a1, u1.x, a2, u2.x);
-        this.y = MathUtils.linearCombination(a1, u1.y, a2, u2.y);
-        this.z = MathUtils.linearCombination(a1, u1.z, a2, u2.z);
+        this.x = MathArrays.linearCombination(a1, u1.x, a2, u2.x);
+        this.y = MathArrays.linearCombination(a1, u1.y, a2, u2.y);
+        this.z = MathArrays.linearCombination(a1, u1.z, a2, u2.z);
     }
 
     /** Linear constructor
@@ -149,9 +150,9 @@ public class Vector3D implements Seriali
      */
     public Vector3D(double a1, Vector3D u1, double a2, Vector3D u2,
                     double a3, Vector3D u3) {
-        this.x = MathUtils.linearCombination(a1, u1.x, a2, u2.x, a3, u3.x);
-        this.y = MathUtils.linearCombination(a1, u1.y, a2, u2.y, a3, u3.y);
-        this.z = MathUtils.linearCombination(a1, u1.z, a2, u2.z, a3, u3.z);
+        this.x = MathArrays.linearCombination(a1, u1.x, a2, u2.x, a3, u3.x);
+        this.y = MathArrays.linearCombination(a1, u1.y, a2, u2.y, a3, u3.y);
+        this.z = MathArrays.linearCombination(a1, u1.z, a2, u2.z, a3, u3.z);
     }
 
     /** Linear constructor
@@ -168,9 +169,9 @@ public class Vector3D implements Seriali
      */
     public Vector3D(double a1, Vector3D u1, double a2, Vector3D u2,
                     double a3, Vector3D u3, double a4, Vector3D u4) {
-        this.x = MathUtils.linearCombination(a1, u1.x, a2, u2.x, a3, u3.x, a4, u4.x);
-        this.y = MathUtils.linearCombination(a1, u1.y, a2, u2.y, a3, u3.y, a4, u4.y);
-        this.z = MathUtils.linearCombination(a1, u1.z, a2, u2.z, a3, u3.z, a4, u4.z);
+        this.x = MathArrays.linearCombination(a1, u1.x, a2, u2.x, a3, u3.x, a4, u4.x);
+        this.y = MathArrays.linearCombination(a1, u1.y, a2, u2.y, a3, u3.y, a4, u4.y);
+        this.z = MathArrays.linearCombination(a1, u1.z, a2, u2.z, a3, u3.z, a4, u4.z);
     }
 
     /** Get the abscissa of the vector.
@@ -422,11 +423,11 @@ public class Vector3D implements Seriali
      * algorithms to preserve accuracy and reduce cancellation effects.
      * It should be very accurate even for nearly orthogonal vectors.
      * </p>
-     * @see MathUtils#linearCombination(double, double, double, double, double, double)
+     * @see MathArrays#linearCombination(double, double, double, double, double, double)
      */
     public double dotProduct(final Vector<Euclidean3D> v) {
         final Vector3D v3 = (Vector3D) v;
-        return MathUtils.linearCombination(x, v3.x, y, v3.y, z, v3.z);
+        return MathArrays.linearCombination(x, v3.x, y, v3.y, z, v3.z);
     }
 
     /** Compute the cross-product of the instance with another vector.
@@ -435,9 +436,9 @@ public class Vector3D implements Seriali
      */
     public Vector3D crossProduct(final Vector<Euclidean3D> v) {
         final Vector3D v3 = (Vector3D) v;
-        return new Vector3D(MathUtils.linearCombination(y, v3.z, -z, v3.y),
-                            MathUtils.linearCombination(z, v3.x, -x, v3.z),
-                            MathUtils.linearCombination(x, v3.y, -y, v3.x));
+        return new Vector3D(MathArrays.linearCombination(y, v3.z, -z, v3.y),
+                            MathArrays.linearCombination(z, v3.x, -x, v3.z),
+                            MathArrays.linearCombination(x, v3.y, -y, v3.x));
     }
 
     /** {@inheritDoc} */

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java Tue Oct 11 22:55:08 2011
@@ -16,7 +16,7 @@
 package org.apache.commons.math.linear;
 
 import java.util.Arrays;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 import org.apache.commons.math.exception.ConvergenceException;
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.util.LocalizedFormats;
@@ -55,7 +55,7 @@ public class PivotingQRDecomposition {
     }
 
     public int[] getOrder() {
-        return MathUtils.copyOf(permutation);
+        return MathArrays.copyOf(permutation);
     }
 
     public PivotingQRDecomposition(RealMatrix matrix) throws ConvergenceException {

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/BOBYQAOptimizer.java Tue Oct 11 22:55:08 2011
@@ -31,7 +31,7 @@ import org.apache.commons.math.linear.Re
 import org.apache.commons.math.optimization.GoalType;
 import org.apache.commons.math.optimization.MultivariateRealOptimizer;
 import org.apache.commons.math.optimization.RealPointValuePair;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Powell's BOBYQA algorithm. This implementation is translated and
@@ -157,8 +157,8 @@ public class BOBYQAOptimizer
                            double[] upperBound,
                            double initialTrustRegionRadius,
                            double stoppingTrustRegionRadius) {
-        this.lowerBound = lowerBound == null ? null : MathUtils.copyOf(lowerBound);
-        this.upperBound = upperBound == null ? null : MathUtils.copyOf(upperBound);
+        this.lowerBound = lowerBound == null ? null : MathArrays.copyOf(lowerBound);
+        this.upperBound = upperBound == null ? null : MathArrays.copyOf(upperBound);
         this.numberOfInterpolationPoints = numberOfInterpolationPoints;
         this.initialTrustRegionRadius = initialTrustRegionRadius;
         this.stoppingTrustRegionRadius = stoppingTrustRegionRadius;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/CMAESOptimizer.java Tue Oct 11 22:55:08 2011
@@ -36,7 +36,7 @@ import org.apache.commons.math.optimizat
 import org.apache.commons.math.optimization.RealPointValuePair;
 import org.apache.commons.math.random.MersenneTwister;
 import org.apache.commons.math.random.RandomGenerator;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * <p>An implementation of the active Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
@@ -395,9 +395,9 @@ public class CMAESOptimizer extends
                 int[] arindex = sortedIndices(fitness);
                 // Calculate new xmean, this is selection and recombination
                 RealMatrix xold = xmean; // for speed up of Eq. (2) and (3)
-                RealMatrix bestArx = selectColumns(arx, MathUtils.copyOf(arindex, mu));
+                RealMatrix bestArx = selectColumns(arx, MathArrays.copyOf(arindex, mu));
                 xmean = bestArx.multiply(weights);
-                RealMatrix bestArz = selectColumns(arz, MathUtils.copyOf(arindex, mu));
+                RealMatrix bestArz = selectColumns(arz, MathArrays.copyOf(arindex, mu));
                 RealMatrix zmean = bestArz.multiply(weights);
                 boolean hsig = updateEvolutionPaths(zmean, xold);
                 if (diagonalOnly <= 0) {
@@ -711,7 +711,7 @@ public class CMAESOptimizer extends
                 // prepare vectors, compute negative updating matrix Cneg
                 int[] arReverseIndex = reverse(arindex);
                 RealMatrix arzneg
-                    = selectColumns(arz, MathUtils.copyOf(arReverseIndex, mu));
+                    = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu));
                 RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
                 int[] idxnorms = sortedIndices(arnorms.getRow(0));
                 RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/PowellOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/PowellOptimizer.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/PowellOptimizer.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/direct/PowellOptimizer.java Tue Oct 11 22:55:08 2011
@@ -18,7 +18,7 @@
 package org.apache.commons.math.optimization.direct;
 
 import org.apache.commons.math.util.FastMath;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 import org.apache.commons.math.analysis.UnivariateRealFunction;
 import org.apache.commons.math.analysis.MultivariateRealFunction;
 import org.apache.commons.math.exception.NumberIsTooSmallException;
@@ -141,7 +141,7 @@ public class PowellOptimizer
             double alphaMin = 0;
 
             for (int i = 0; i < n; i++) {
-                final double[] d = MathUtils.copyOf(direc[i]);
+                final double[] d = MathArrays.copyOf(direc[i]);
 
                 fX2 = fVal;
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/clustering/EuclideanIntegerPoint.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/clustering/EuclideanIntegerPoint.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/clustering/EuclideanIntegerPoint.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/clustering/EuclideanIntegerPoint.java Tue Oct 11 22:55:08 2011
@@ -20,7 +20,7 @@ package org.apache.commons.math.stat.clu
 import java.io.Serializable;
 import java.util.Collection;
 
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * A simple implementation of {@link Clusterable} for points with integer coordinates.
@@ -54,7 +54,7 @@ public class EuclideanIntegerPoint imple
 
     /** {@inheritDoc} */
     public double distanceFrom(final EuclideanIntegerPoint p) {
-        return MathUtils.distance(point, p.getPoint());
+        return MathArrays.distance(point, p.getPoint());
     }
 
     /** {@inheritDoc} */

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/MillerUpdatingRegression.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/MillerUpdatingRegression.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/MillerUpdatingRegression.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/MillerUpdatingRegression.java Tue Oct 11 22:55:08 2011
@@ -20,6 +20,7 @@ import java.util.Arrays;
 import org.apache.commons.math.exception.util.LocalizedFormats;
 import org.apache.commons.math.util.FastMath;
 import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * <p>This class is a concrete implementation of the {@link UpdatingMultipleLinearRegression} interface.</p>
@@ -186,7 +187,7 @@ public class MillerUpdatingRegression im
                     x.length, nvars);
         }
         if (!this.hasIntercept) {
-            include(MathUtils.copyOf(x, x.length), 1.0, y);
+            include(MathArrays.copyOf(x, x.length), 1.0, y);
         } else {
             double[] tmp = new double[x.length + 1];
             System.arraycopy(x, 0, tmp, 1, x.length);
@@ -918,7 +919,7 @@ public class MillerUpdatingRegression im
      * @return int[] with the current order of the regressors
      */
     public int[] getOrderOfRegressors(){
-        return MathUtils.copyOf(vorder);
+        return MathArrays.copyOf(vorder);
     }
 
     /**

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/RegressionResults.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/RegressionResults.java?rev=1182134&r1=1182133&r2=1182134&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/RegressionResults.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/regression/RegressionResults.java Tue Oct 11 22:55:08 2011
@@ -19,7 +19,7 @@ package org.apache.commons.math.stat.reg
 import java.io.Serializable;
 import java.util.Arrays;
 import org.apache.commons.math.util.FastMath;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.MathArrays;
 
 /**
  * Results of a Multiple Linear Regression model fit.
@@ -98,10 +98,10 @@ public class RegressionResults implement
             final boolean containsConstant,
             final boolean copyData) {
         if (copyData) {
-            this.parameters = MathUtils.copyOf(parameters);
+            this.parameters = MathArrays.copyOf(parameters);
             this.varCovData = new double[varcov.length][];
             for (int i = 0; i < varcov.length; i++) {
-                this.varCovData[i] = MathUtils.copyOf(varcov[i]);
+                this.varCovData[i] = MathArrays.copyOf(varcov[i]);
             }
         } else {
             this.parameters = parameters;
@@ -170,7 +170,7 @@ public class RegressionResults implement
         if (this.parameters == null) {
             return null;
         }
-        return MathUtils.copyOf(parameters);
+        return MathArrays.copyOf(parameters);
     }
 
     /**

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathArrays.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathArrays.java?rev=1182134&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathArrays.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathArrays.java Tue Oct 11 22:55:08 2011
@@ -0,0 +1,925 @@
+/*
+ * 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.util;
+
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Comparator;
+import java.util.Collections;
+
+import org.apache.commons.math.exception.DimensionMismatchException;
+import org.apache.commons.math.exception.MathInternalError;
+import org.apache.commons.math.exception.NonMonotonicSequenceException;
+import org.apache.commons.math.exception.NullArgumentException;
+
+/**
+ * Arrays utilities.
+ *
+ * @since 3.0
+ * @version $Id$
+ */
+public class MathArrays {
+    /** Factor used for splitting double numbers: n = 2^27 + 1 (i.e. {@value}). */
+    private static final int SPLIT_FACTOR = 0x8000001;
+
+    /**
+     * Private constructor.
+     */
+    private MathArrays() {}
+
+    /**
+     * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
+     *
+     * @param p1 the first point
+     * @param p2 the second point
+     * @return the L<sub>1</sub> distance between the two points
+     */
+    public static double distance1(double[] p1, double[] p2) {
+        double sum = 0;
+        for (int i = 0; i < p1.length; i++) {
+            sum += FastMath.abs(p1[i] - p2[i]);
+        }
+        return sum;
+    }
+
+    /**
+     * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
+     *
+     * @param p1 the first point
+     * @param p2 the second point
+     * @return the L<sub>1</sub> distance between the two points
+     */
+    public static int distance1(int[] p1, int[] p2) {
+      int sum = 0;
+      for (int i = 0; i < p1.length; i++) {
+          sum += FastMath.abs(p1[i] - p2[i]);
+      }
+      return sum;
+    }
+
+    /**
+     * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
+     *
+     * @param p1 the first point
+     * @param p2 the second point
+     * @return the L<sub>2</sub> distance between the two points
+     */
+    public static double distance(double[] p1, double[] p2) {
+        double sum = 0;
+        for (int i = 0; i < p1.length; i++) {
+            final double dp = p1[i] - p2[i];
+            sum += dp * dp;
+        }
+        return FastMath.sqrt(sum);
+    }
+
+    /**
+     * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
+     *
+     * @param p1 the first point
+     * @param p2 the second point
+     * @return the L<sub>2</sub> distance between the two points
+     */
+    public static double distance(int[] p1, int[] p2) {
+      double sum = 0;
+      for (int i = 0; i < p1.length; i++) {
+          final double dp = p1[i] - p2[i];
+          sum += dp * dp;
+      }
+      return FastMath.sqrt(sum);
+    }
+
+    /**
+     * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
+     *
+     * @param p1 the first point
+     * @param p2 the second point
+     * @return the L<sub>&infin;</sub> distance between the two points
+     */
+    public static double distanceInf(double[] p1, double[] p2) {
+        double max = 0;
+        for (int i = 0; i < p1.length; i++) {
+            max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
+        }
+        return max;
+    }
+
+    /**
+     * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
+     *
+     * @param p1 the first point
+     * @param p2 the second point
+     * @return the L<sub>&infin;</sub> distance between the two points
+     */
+    public static int distanceInf(int[] p1, int[] p2) {
+        int max = 0;
+        for (int i = 0; i < p1.length; i++) {
+            max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
+        }
+        return max;
+    }
+
+    /**
+     * Specification of ordering direction.
+     */
+    public static enum OrderDirection {
+        /** Constant for increasing direction. */
+        INCREASING,
+        /** Constant for decreasing direction. */
+        DECREASING
+    }
+
+    /**
+     * Check that an array is monotonically increasing or decreasing.
+     *
+     * @param val Values.
+     * @param dir Ordering direction.
+     * @param strict Whether the order should be strict.
+     * @return {@code true} if sorted, {@code false} otherwise.
+     */
+    public static boolean isMonotonic(Comparable[] val,
+                                      OrderDirection dir,
+                                      boolean strict) {
+        Comparable previous = val[0];
+        final int max = val.length;
+        int comp;
+        for (int i = 1; i < max; i++) {
+            switch (dir) {
+            case INCREASING:
+                comp = -val[i].compareTo(previous);
+                if (strict) {
+                    if (0 <= comp) {
+                        return false;
+                    }
+                } else {
+                    if ( comp > 0) {
+                        return false;
+                    }
+                }
+                break;
+            case DECREASING:
+                comp = val[i].compareTo(previous);
+                if (strict) {
+                    if (comp >= 0) {
+                        return false;
+                    }
+                } else {
+                    if (comp > 0) {
+                       return false;
+                    }
+                }
+                break;
+            default:
+                // Should never happen.
+                throw new MathInternalError();
+            }
+
+            previous = val[i];
+        }
+        return true;
+    }
+
+    /**
+     * Check that an array is monotonically increasing or decreasing.
+     *
+     * @param val Values.
+     * @param dir Ordering direction.
+     * @param strict Whether the order should be strict.
+     * @return {@code true} if sorted, {@code false} otherwise.
+     */
+    public static boolean isMonotonic(double[] val,
+                                      OrderDirection dir,
+                                      boolean strict) {
+        return checkOrder(val, dir, strict, false);
+    }
+
+    /**
+     * Check that the given array is sorted.
+     *
+     * @param val Values.
+     * @param dir Ordering direction.
+     * @param strict Whether the order should be strict.
+     * @param abort Whether to throw an exception if the check fails.
+     * @return {@code true} if the array is sorted.
+     * @throws NonMonotonicSequenceException if the array is not sorted
+     * and {@code abort} is {@code true}.
+     */
+    public static boolean checkOrder(double[] val, OrderDirection dir,
+                                     boolean strict, boolean abort) {
+        double previous = val[0];
+        final int max = val.length;
+
+        int index;
+        ITEM:
+        for (index = 1; index < max; index++) {
+            switch (dir) {
+            case INCREASING:
+                if (strict) {
+                    if (val[index] <= previous) {
+                        break ITEM;
+                    }
+                } else {
+                    if (val[index] < previous) {
+                        break ITEM;
+                    }
+                }
+                break;
+            case DECREASING:
+                if (strict) {
+                    if (val[index] >= previous) {
+                        break ITEM;
+                    }
+                } else {
+                    if (val[index] > previous) {
+                        break ITEM;
+                    }
+                }
+                break;
+            default:
+                // Should never happen.
+                throw new MathInternalError();
+            }
+
+            previous = val[index];
+        }
+
+        if (index == max) {
+            // Loop completed.
+            return true;
+        }
+
+        // Loop early exit means wrong ordering.
+        if (abort) {
+            throw new NonMonotonicSequenceException(val[index], previous, index, dir, strict);
+        } else {
+            return false;
+        }
+    }
+
+    /**
+     * Check that the given array is sorted.
+     *
+     * @param val Values.
+     * @param dir Ordering direction.
+     * @param strict Whether the order should be strict.
+     * @throws NonMonotonicSequenceException if the array is not sorted.
+     * @since 2.2
+     */
+    public static void checkOrder(double[] val, OrderDirection dir,
+                                  boolean strict) {
+        checkOrder(val, dir, strict, true);
+    }
+
+    /**
+     * Check that the given array is sorted in strictly increasing order.
+     *
+     * @param val Values.
+     * @throws NonMonotonicSequenceException if the array is not sorted.
+     * @since 2.2
+     */
+    public static void checkOrder(double[] val) {
+        checkOrder(val, OrderDirection.INCREASING, true);
+    }
+
+    /**
+     * Returns the Cartesian norm (2-norm), handling both overflow and underflow.
+     * Translation of the minpack enorm subroutine.
+     *
+     * The redistribution policy for MINPACK is available
+     * <a href="http://www.netlib.org/minpack/disclaimer">here</a>, for
+     * convenience, it is reproduced below.</p>
+     *
+     * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
+     * <tr><td>
+     *    Minpack Copyright Notice (1999) University of Chicago.
+     *    All rights reserved
+     * </td></tr>
+     * <tr><td>
+     * Redistribution and use in source and binary forms, with or without
+     * modification, are permitted provided that the following conditions
+     * are met:
+     * <ol>
+     *  <li>Redistributions of source code must retain the above copyright
+     *      notice, this list of conditions and the following disclaimer.</li>
+     * <li>Redistributions in binary form must reproduce the above
+     *     copyright notice, this list of conditions and the following
+     *     disclaimer in the documentation and/or other materials provided
+     *     with the distribution.</li>
+     * <li>The end-user documentation included with the redistribution, if any,
+     *     must include the following acknowledgment:
+     *     {@code This product includes software developed by the University of
+     *           Chicago, as Operator of Argonne National Laboratory.}
+     *     Alternately, this acknowledgment may appear in the software itself,
+     *     if and wherever such third-party acknowledgments normally appear.</li>
+     * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
+     *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
+     *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
+     *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
+     *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
+     *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
+     *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
+     *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
+     *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
+     *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
+     *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
+     *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
+     *     BE CORRECTED.</strong></li>
+     * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
+     *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
+     *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
+     *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
+     *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
+     *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
+     *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
+     *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
+     *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
+     *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
+     * <ol></td></tr>
+     * </table>
+     *
+     * @param v Vector of doubles.
+     * @return the 2-norm of the vector.
+     * @since 2.2
+     */
+    public static double safeNorm(double[] v) {
+        double rdwarf = 3.834e-20;
+        double rgiant = 1.304e+19;
+        double s1 = 0;
+        double s2 = 0;
+        double s3 = 0;
+        double x1max = 0;
+        double x3max = 0;
+        double floatn = (double) v.length;
+        double agiant = rgiant / floatn;
+        for (int i = 0; i < v.length; i++) {
+            double xabs = Math.abs(v[i]);
+            if (xabs < rdwarf || xabs > agiant) {
+                if (xabs > rdwarf) {
+                    if (xabs > x1max) {
+                        double r = x1max / xabs;
+                        s1= 1 + s1 * r * r;
+                        x1max = xabs;
+                    } else {
+                        double r = xabs / x1max;
+                        s1 += r * r;
+                    }
+                } else {
+                    if (xabs > x3max) {
+                        double r = x3max / xabs;
+                        s3= 1 + s3 * r * r;
+                        x3max = xabs;
+                    } else {
+                        if (xabs != 0) {
+                            double r = xabs / x3max;
+                            s3 += r * r;
+                        }
+                    }
+                }
+            } else {
+                s2 += xabs * xabs;
+            }
+        }
+        double norm;
+        if (s1 != 0) {
+            norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
+        } else {
+            if (s2 == 0) {
+                norm = x3max * Math.sqrt(s3);
+            } else {
+                if (s2 >= x3max) {
+                    norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
+                } else {
+                    norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
+                }
+            }
+        }
+        return norm;
+    }
+
+    /**
+     * Sort an array in ascending order in place and perform the same reordering
+     * of entries on other arrays. For example, if
+     * {@code x = [3, 1, 2], y = [1, 2, 3]} and {@code z = [0, 5, 7]}, then
+     * {@code sortInPlace(x, y, z)} will update {@code x} to {@code [1, 2, 3]},
+     * {@code y} to {@code [2, 3, 1]} and {@code z} to {@code [5, 7, 0]}.
+     *
+     * @param x Array to be sorted and used as a pattern for permutation
+     * of the other arrays.
+     * @param yList Set of arrays whose permutations of entries will follow
+     * those performed on {@code x}.
+     * @throws DimensionMismatchException if any {@code y} is not the same
+     * size as {@code x}.
+     * @throws NullArgumentException if {@code x} or any {@code y} is null.
+     * @since 3.0
+     */
+    public static void sortInPlace(double[] x,
+                                   double[] ... yList) {
+        sortInPlace(x, OrderDirection.INCREASING, yList);
+    }
+
+    /**
+     * Sort an array in place and perform the same reordering of entries on
+     * other arrays.  This method works the same as
+     * {@link #sortInPlace(double[], double[] ...)}, but allows the order of
+     * the sort to be provided in the {@code dir} parameter.
+     *
+     * @param x Array to be sorted and used as a pattern for permutation
+     * of the other arrays.
+     * @param dir Order direction.
+     * @param yList Set of arrays whose permutations of entries will follow
+     * those performed on {@code x}.
+     * @throws DimensionMismatchException if any {@code y} is not the same
+     * size as {@code x}.
+     * @throws NullArgumentException if {@code x} or any {@code y} is null
+     * @since 3.0
+     */
+    public static void sortInPlace(double[] x,
+                                   final OrderDirection dir,
+                                   double[] ... yList) {
+        if (x == null) {
+            throw new NullArgumentException();
+        }
+
+        final int len = x.length;
+        final List<Pair<Double, double[]>> list
+            = new ArrayList<Pair<Double, double[]>>(len);
+
+        final int yListLen = yList.length;
+        for (int i = 0; i < len; i++) {
+            final double[] yValues = new double[yListLen];
+            for (int j = 0; j < yListLen; j++) {
+                double[] y = yList[j];
+                if (y == null) {
+                    throw new NullArgumentException();
+                }
+                if (y.length != len) {
+                    throw new DimensionMismatchException(y.length, len);
+                }
+                yValues[j] = y[i];
+            }
+            list.add(new Pair<Double, double[]>(x[i], yValues));
+        }
+
+        final Comparator<Pair<Double, double[]>> comp
+            = new Comparator<Pair<Double, double[]>>() {
+            public int compare(Pair<Double, double[]> o1,
+                               Pair<Double, double[]> o2) {
+                int val;
+                switch (dir) {
+                case INCREASING:
+                    val = o1.getKey().compareTo(o2.getKey());
+                break;
+                case DECREASING:
+                    val = o2.getKey().compareTo(o1.getKey());
+                break;
+                default:
+                    // Should never happen.
+                    throw new MathInternalError();
+                }
+                return val;
+            }
+        };
+
+        Collections.sort(list, comp);
+
+        for (int i = 0; i < len; i++) {
+            final Pair<Double, double[]> e = list.get(i);
+            x[i] = e.getKey();
+            final double[] yValues = e.getValue();
+            for (int j = 0; j < yListLen; j++) {
+                yList[j][i] = yValues[j];
+            }
+        }
+    }
+
+    /**
+     * Creates a copy of the {@code source} array.
+     *
+     * @param source Array to be copied.
+     * @return the copied array.
+     */
+     public static int[] copyOf(int[] source) {
+         return copyOf(source, source.length);
+     }
+
+    /**
+     * Creates a copy of the {@code source} array.
+     *
+     * @param source Array to be copied.
+     * @return the copied array.
+     */
+     public static double[] copyOf(double[] source) {
+         return copyOf(source, source.length);
+     }
+
+    /**
+     * Creates a copy of the {@code source} array.
+     *
+     * @param source Array to be copied.
+     * @param len Number of entries to copy. If smaller then the source
+     * length, the copy will be truncated, if larger it will padded with
+     * zeroes.
+     * @return the copied array.
+     */
+    public static int[] copyOf(int[] source, int len) {
+         final int[] output = new int[len];
+         System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length));
+         return output;
+     }
+
+    /**
+     * Creates a copy of the {@code source} array.
+     *
+     * @param source Array to be copied.
+     * @param len Number of entries to copy. If smaller then the source
+     * length, the copy will be truncated, if larger it will padded with
+     * zeroes.
+     * @return the copied array.
+     */
+    public static double[] copyOf(double[] source, int len) {
+         final double[] output = new double[len];
+         System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length));
+         return output;
+     }
+
+    /**
+     * Compute a linear combination accurately.
+     * This method computes the sum of the products
+     * <code>a<sub>i</sub> b<sub>i</sub></code> to high accuracy.
+     * It does so by using specific multiplication and addition algorithms to
+     * preserve accuracy and reduce cancellation effects.
+     * <br/>
+     * It is based on the 2005 paper
+     * <a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
+     * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
+     * and Shin'ichi Oishi published in SIAM J. Sci. Comput.
+     *
+     * @param a Factors.
+     * @param b Factors.
+     * @return <code>&Sigma;<sub>i</sub> a<sub>i</sub> b<sub>i</sub></code>.
+     */
+    public static double linearCombination(final double[] a, final double[] b) {
+        final int len = a.length;
+        if (len != b.length) {
+            throw new DimensionMismatchException(len, b.length);
+        }
+
+        final double[] prodHigh = new double[len];
+        double prodLowSum = 0;
+
+        for (int i = 0; i < len; i++) {
+            final double ai = a[i];
+            final double ca = SPLIT_FACTOR * ai;
+            final double aHigh = ca - (ca - ai);
+            final double aLow = ai - aHigh;
+
+            final double bi = b[i];
+            final double cb = SPLIT_FACTOR * bi;
+            final double bHigh = cb - (cb - bi);
+            final double bLow = bi - bHigh;
+            prodHigh[i] = ai * bi;
+            final double prodLow = aLow * bLow - (((prodHigh[i] -
+                                                    aHigh * bHigh) -
+                                                   aLow * bHigh) -
+                                                  aHigh * bLow);
+            prodLowSum += prodLow;
+        }
+
+
+        final double prodHighCur = prodHigh[0];
+        double prodHighNext = prodHigh[1];
+        double sHighPrev = prodHighCur + prodHighNext;
+        double sPrime = sHighPrev - prodHighNext;
+        double sLowSum = (prodHighNext - (sHighPrev - sPrime)) + (prodHighCur - sPrime);
+
+        final int lenMinusOne = len - 1;
+        for (int i = 1; i < lenMinusOne; i++) {
+            prodHighNext = prodHigh[i + 1];
+            final double sHighCur = sHighPrev + prodHighNext;
+            sPrime = sHighCur - prodHighNext;
+            sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime);
+            sHighPrev = sHighCur;
+        }
+
+        double result = sHighPrev + (prodLowSum + sLowSum);
+
+        if (Double.isNaN(result)) {
+            // either we have split infinite numbers or some coefficients were NaNs,
+            // just rely on the naive implementation and let IEEE754 handle this
+            result = 0;
+            for (int i = 0; i < len; ++i) {
+                result += a[i] * b[i];
+            }
+        }
+
+        return result;
+    }
+
+    /**
+     * Compute a linear combination accurately.
+     * <p>
+     * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
+     * a<sub>2</sub>&times;b<sub>2</sub> to high accuracy. It does
+     * so by using specific multiplication and addition algorithms to
+     * preserve accuracy and reduce cancellation effects. It is based
+     * on the 2005 paper <a
+     * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
+     * Accurate Sum and Dot Product</a> by Takeshi Ogita,
+     * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
+     * </p>
+     * @param a1 first factor of the first term
+     * @param b1 second factor of the first term
+     * @param a2 first factor of the second term
+     * @param b2 second factor of the second term
+     * @return a<sub>1</sub>&times;b<sub>1</sub> +
+     * a<sub>2</sub>&times;b<sub>2</sub>
+     * @see #linearCombination(double, double, double, double, double, double)
+     * @see #linearCombination(double, double, double, double, double, double, double, double)
+     */
+    public static double linearCombination(final double a1, final double b1,
+                                           final double a2, final double b2) {
+
+        // the code below is split in many additions/subtractions that may
+        // appear redundant. However, they should NOT be simplified, as they
+        // use IEEE754 floating point arithmetic rounding properties.
+        // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
+        // The variable naming conventions are that xyzHigh contains the most significant
+        // bits of xyz and xyzLow contains its least significant bits. So theoretically
+        // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
+        // be represented in only one double precision number so we preserve two numbers
+        // to hold it as long as we can, combining the high and low order bits together
+        // only at the end, after cancellation may have occurred on high order bits
+
+        // split a1 and b1 as two 26 bits numbers
+        final double ca1        = SPLIT_FACTOR * a1;
+        final double a1High     = ca1 - (ca1 - a1);
+        final double a1Low      = a1 - a1High;
+        final double cb1        = SPLIT_FACTOR * b1;
+        final double b1High     = cb1 - (cb1 - b1);
+        final double b1Low      = b1 - b1High;
+
+        // accurate multiplication a1 * b1
+        final double prod1High  = a1 * b1;
+        final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
+
+        // split a2 and b2 as two 26 bits numbers
+        final double ca2        = SPLIT_FACTOR * a2;
+        final double a2High     = ca2 - (ca2 - a2);
+        final double a2Low      = a2 - a2High;
+        final double cb2        = SPLIT_FACTOR * b2;
+        final double b2High     = cb2 - (cb2 - b2);
+        final double b2Low      = b2 - b2High;
+
+        // accurate multiplication a2 * b2
+        final double prod2High  = a2 * b2;
+        final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
+
+        // accurate addition a1 * b1 + a2 * b2
+        final double s12High    = prod1High + prod2High;
+        final double s12Prime   = s12High - prod2High;
+        final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
+
+        // final rounding, s12 may have suffered many cancellations, we try
+        // to recover some bits from the extra words we have saved up to now
+        double result = s12High + (prod1Low + prod2Low + s12Low);
+
+        if (Double.isNaN(result)) {
+            // either we have split infinite numbers or some coefficients were NaNs,
+            // just rely on the naive implementation and let IEEE754 handle this
+            result = a1 * b1 + a2 * b2;
+        }
+
+        return result;
+    }
+
+    /**
+     * Compute a linear combination accurately.
+     * <p>
+     * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
+     * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
+     * to high accuracy. It does so by using specific multiplication and
+     * addition algorithms to preserve accuracy and reduce cancellation effects.
+     * It is based on the 2005 paper <a
+     * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
+     * Accurate Sum and Dot Product</a> by Takeshi Ogita,
+     * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
+     * </p>
+     * @param a1 first factor of the first term
+     * @param b1 second factor of the first term
+     * @param a2 first factor of the second term
+     * @param b2 second factor of the second term
+     * @param a3 first factor of the third term
+     * @param b3 second factor of the third term
+     * @return a<sub>1</sub>&times;b<sub>1</sub> +
+     * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
+     * @see #linearCombination(double, double, double, double)
+     * @see #linearCombination(double, double, double, double, double, double, double, double)
+     */
+    public static double linearCombination(final double a1, final double b1,
+                                           final double a2, final double b2,
+                                           final double a3, final double b3) {
+
+        // the code below is split in many additions/subtractions that may
+        // appear redundant. However, they should NOT be simplified, as they
+        // do use IEEE754 floating point arithmetic rounding properties.
+        // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
+        // The variables naming conventions are that xyzHigh contains the most significant
+        // bits of xyz and xyzLow contains its least significant bits. So theoretically
+        // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
+        // be represented in only one double precision number so we preserve two numbers
+        // to hold it as long as we can, combining the high and low order bits together
+        // only at the end, after cancellation may have occurred on high order bits
+
+        // split a1 and b1 as two 26 bits numbers
+        final double ca1        = SPLIT_FACTOR * a1;
+        final double a1High     = ca1 - (ca1 - a1);
+        final double a1Low      = a1 - a1High;
+        final double cb1        = SPLIT_FACTOR * b1;
+        final double b1High     = cb1 - (cb1 - b1);
+        final double b1Low      = b1 - b1High;
+
+        // accurate multiplication a1 * b1
+        final double prod1High  = a1 * b1;
+        final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
+
+        // split a2 and b2 as two 26 bits numbers
+        final double ca2        = SPLIT_FACTOR * a2;
+        final double a2High     = ca2 - (ca2 - a2);
+        final double a2Low      = a2 - a2High;
+        final double cb2        = SPLIT_FACTOR * b2;
+        final double b2High     = cb2 - (cb2 - b2);
+        final double b2Low      = b2 - b2High;
+
+        // accurate multiplication a2 * b2
+        final double prod2High  = a2 * b2;
+        final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
+
+        // split a3 and b3 as two 26 bits numbers
+        final double ca3        = SPLIT_FACTOR * a3;
+        final double a3High     = ca3 - (ca3 - a3);
+        final double a3Low      = a3 - a3High;
+        final double cb3        = SPLIT_FACTOR * b3;
+        final double b3High     = cb3 - (cb3 - b3);
+        final double b3Low      = b3 - b3High;
+
+        // accurate multiplication a3 * b3
+        final double prod3High  = a3 * b3;
+        final double prod3Low   = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
+
+        // accurate addition a1 * b1 + a2 * b2
+        final double s12High    = prod1High + prod2High;
+        final double s12Prime   = s12High - prod2High;
+        final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
+
+        // accurate addition a1 * b1 + a2 * b2 + a3 * b3
+        final double s123High   = s12High + prod3High;
+        final double s123Prime  = s123High - prod3High;
+        final double s123Low    = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);
+
+        // final rounding, s123 may have suffered many cancellations, we try
+        // to recover some bits from the extra words we have saved up to now
+        double result = s123High + (prod1Low + prod2Low + prod3Low + s12Low + s123Low);
+
+        if (Double.isNaN(result)) {
+            // either we have split infinite numbers or some coefficients were NaNs,
+            // just rely on the naive implementation and let IEEE754 handle this
+            result = a1 * b1 + a2 * b2 + a3 * b3;
+        }
+
+        return result;
+    }
+
+    /**
+     * Compute a linear combination accurately.
+     * <p>
+     * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
+     * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub> +
+     * a<sub>4</sub>&times;b<sub>4</sub>
+     * to high accuracy. It does so by using specific multiplication and
+     * addition algorithms to preserve accuracy and reduce cancellation effects.
+     * It is based on the 2005 paper <a
+     * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
+     * Accurate Sum and Dot Product</a> by Takeshi Ogita,
+     * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
+     * </p>
+     * @param a1 first factor of the first term
+     * @param b1 second factor of the first term
+     * @param a2 first factor of the second term
+     * @param b2 second factor of the second term
+     * @param a3 first factor of the third term
+     * @param b3 second factor of the third term
+     * @param a4 first factor of the third term
+     * @param b4 second factor of the third term
+     * @return a<sub>1</sub>&times;b<sub>1</sub> +
+     * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub> +
+     * a<sub>4</sub>&times;b<sub>4</sub>
+     * @see #linearCombination(double, double, double, double)
+     * @see #linearCombination(double, double, double, double, double, double)
+     */
+    public static double linearCombination(final double a1, final double b1,
+                                           final double a2, final double b2,
+                                           final double a3, final double b3,
+                                           final double a4, final double b4) {
+
+        // the code below is split in many additions/subtractions that may
+        // appear redundant. However, they should NOT be simplified, as they
+        // do use IEEE754 floating point arithmetic rounding properties.
+        // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
+        // The variables naming conventions are that xyzHigh contains the most significant
+        // bits of xyz and xyzLow contains its least significant bits. So theoretically
+        // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
+        // be represented in only one double precision number so we preserve two numbers
+        // to hold it as long as we can, combining the high and low order bits together
+        // only at the end, after cancellation may have occurred on high order bits
+
+        // split a1 and b1 as two 26 bits numbers
+        final double ca1        = SPLIT_FACTOR * a1;
+        final double a1High     = ca1 - (ca1 - a1);
+        final double a1Low      = a1 - a1High;
+        final double cb1        = SPLIT_FACTOR * b1;
+        final double b1High     = cb1 - (cb1 - b1);
+        final double b1Low      = b1 - b1High;
+
+        // accurate multiplication a1 * b1
+        final double prod1High  = a1 * b1;
+        final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
+
+        // split a2 and b2 as two 26 bits numbers
+        final double ca2        = SPLIT_FACTOR * a2;
+        final double a2High     = ca2 - (ca2 - a2);
+        final double a2Low      = a2 - a2High;
+        final double cb2        = SPLIT_FACTOR * b2;
+        final double b2High     = cb2 - (cb2 - b2);
+        final double b2Low      = b2 - b2High;
+
+        // accurate multiplication a2 * b2
+        final double prod2High  = a2 * b2;
+        final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
+
+        // split a3 and b3 as two 26 bits numbers
+        final double ca3        = SPLIT_FACTOR * a3;
+        final double a3High     = ca3 - (ca3 - a3);
+        final double a3Low      = a3 - a3High;
+        final double cb3        = SPLIT_FACTOR * b3;
+        final double b3High     = cb3 - (cb3 - b3);
+        final double b3Low      = b3 - b3High;
+
+        // accurate multiplication a3 * b3
+        final double prod3High  = a3 * b3;
+        final double prod3Low   = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
+
+        // split a4 and b4 as two 26 bits numbers
+        final double ca4        = SPLIT_FACTOR * a4;
+        final double a4High     = ca4 - (ca4 - a4);
+        final double a4Low      = a4 - a4High;
+        final double cb4        = SPLIT_FACTOR * b4;
+        final double b4High     = cb4 - (cb4 - b4);
+        final double b4Low      = b4 - b4High;
+
+        // accurate multiplication a4 * b4
+        final double prod4High  = a4 * b4;
+        final double prod4Low   = a4Low * b4Low - (((prod4High - a4High * b4High) - a4Low * b4High) - a4High * b4Low);
+
+        // accurate addition a1 * b1 + a2 * b2
+        final double s12High    = prod1High + prod2High;
+        final double s12Prime   = s12High - prod2High;
+        final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
+
+        // accurate addition a1 * b1 + a2 * b2 + a3 * b3
+        final double s123High   = s12High + prod3High;
+        final double s123Prime  = s123High - prod3High;
+        final double s123Low    = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);
+
+        // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4
+        final double s1234High  = s123High + prod4High;
+        final double s1234Prime = s1234High - prod4High;
+        final double s1234Low   = (prod4High - (s1234High - s1234Prime)) + (s123High - s1234Prime);
+
+        // final rounding, s1234 may have suffered many cancellations, we try
+        // to recover some bits from the extra words we have saved up to now
+        double result = s1234High + (prod1Low + prod2Low + prod3Low + prod4Low + s12Low + s123Low + s1234Low);
+
+        if (Double.isNaN(result)) {
+            // either we have split infinite numbers or some coefficients were NaNs,
+            // just rely on the naive implementation and let IEEE754 handle this
+            result = a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4;
+        }
+
+        return result;
+    }
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathArrays.java
------------------------------------------------------------------------------
    svn:eol-style = native



Mime
View raw message