Return-Path: Delivered-To: apmail-commons-commits-archive@minotaur.apache.org Received: (qmail 21655 invoked from network); 16 Oct 2009 14:52:33 -0000 Received: from hermes.apache.org (HELO mail.apache.org) (140.211.11.3) by minotaur.apache.org with SMTP; 16 Oct 2009 14:52:33 -0000 Received: (qmail 59857 invoked by uid 500); 16 Oct 2009 14:52:32 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 59769 invoked by uid 500); 16 Oct 2009 14:52:32 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 59760 invoked by uid 99); 16 Oct 2009 14:52:32 -0000 Received: from nike.apache.org (HELO nike.apache.org) (192.87.106.230) by apache.org (qpsmtpd/0.29) with ESMTP; Fri, 16 Oct 2009 14:52:32 +0000 X-ASF-Spam-Status: No, hits=-2000.0 required=10.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Fri, 16 Oct 2009 14:52:19 +0000 Received: by eris.apache.org (Postfix, from userid 65534) id ADB5C23888E5; Fri, 16 Oct 2009 14:51:56 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r825919 [1/2] - in /commons/proper/math/trunk: ./ src/main/java/org/apache/commons/math/analysis/integration/ src/main/java/org/apache/commons/math/analysis/interpolation/ src/main/java/org/apache/commons/math/analysis/polynomials/ src/main... Date: Fri, 16 Oct 2009 14:51:56 -0000 To: commits@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.8 Message-Id: <20091016145156.ADB5C23888E5@eris.apache.org> X-Virus-Checked: Checked by ClamAV on apache.org Author: luc Date: Fri Oct 16 14:51:55 2009 New Revision: 825919 URL: http://svn.apache.org/viewvc?rev=825919&view=rev Log: tighten checkstyle rules: declaring multiple variables in one statement is now forbidden Modified: commons/proper/math/trunk/checkstyle.xml commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/SimpsonIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/TrapezoidIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/DividedDifferenceInterpolator.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/MicrosphereInterpolatingFunction.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/PolynomialFunctionNewtonForm.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/RiddersSolver.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockFieldMatrix.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockRealMatrix.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.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/LevenbergMarquardtOptimizer.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastCosineTransformer.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastSineTransformer.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/OpenIntToDoubleHashMap.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/OpenIntToFieldHashMap.java Modified: commons/proper/math/trunk/checkstyle.xml URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/checkstyle.xml?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/checkstyle.xml (original) +++ commons/proper/math/trunk/checkstyle.xml Fri Oct 16 14:51:55 2009 @@ -122,15 +122,24 @@ - + - - + + - + + + + + + Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/SimpsonIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/SimpsonIntegrator.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/SimpsonIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/SimpsonIntegrator.java Fri Oct 16 14:51:55 2009 @@ -66,25 +66,22 @@ final double min, final double max) throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { - int i = 1; - double s, olds, t, oldt; - clearResult(); verifyInterval(min, max); verifyIterationCount(); TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); if (minimalIterationCount == 1) { - s = (4 * qtrap.stage(f, min, max, 1) - qtrap.stage(f, min, max, 0)) / 3.0; + final double s = (4 * qtrap.stage(f, min, max, 1) - qtrap.stage(f, min, max, 0)) / 3.0; setResult(s, 1); return result; } // Simpson's rule requires at least two trapezoid stages. - olds = 0; - oldt = qtrap.stage(f, min, max, 0); - while (i <= maximalIterationCount) { - t = qtrap.stage(f, min, max, i); - s = (4 * t - oldt) / 3.0; + double olds = 0; + double oldt = qtrap.stage(f, min, max, 0); + for (int i = 1; i <= maximalIterationCount; ++i) { + final double t = qtrap.stage(f, min, max, i); + final double s = (4 * t - oldt) / 3.0; if (i >= minimalIterationCount) { final double delta = Math.abs(s - olds); final double rLimit = @@ -96,7 +93,6 @@ } olds = s; oldt = t; - i++; } throw new MaxIterationsExceededException(maximalIterationCount); } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/TrapezoidIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/TrapezoidIntegrator.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/TrapezoidIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/integration/TrapezoidIntegrator.java Fri Oct 16 14:51:55 2009 @@ -77,17 +77,15 @@ final double min, final double max, final int n) throws FunctionEvaluationException { - long i, np; - double x, spacing, sum = 0; - if (n == 0) { s = 0.5 * (max - min) * (f.value(min) + f.value(max)); return s; } else { - np = 1L << (n-1); // number of new points in this stage - spacing = (max - min) / np; // spacing between adjacent new points - x = min + 0.5 * spacing; // the first new point - for (i = 0; i < np; i++) { + final long np = 1L << (n-1); // number of new points in this stage + double sum = 0; + final double spacing = (max - min) / np; // spacing between adjacent new points + double x = min + 0.5 * spacing; // the first new point + for (long i = 0; i < np; i++) { sum += f.value(x); x += spacing; } @@ -109,16 +107,13 @@ final double min, final double max) throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { - int i = 1; - double t, oldt; - clearResult(); verifyInterval(min, max); verifyIterationCount(); - oldt = stage(f, min, max, 0); - while (i <= maximalIterationCount) { - t = stage(f, min, max, i); + double oldt = stage(f, min, max, 0); + for (int i = 1; i <= maximalIterationCount; ++i) { + final double t = stage(f, min, max, i); if (i >= minimalIterationCount) { final double delta = Math.abs(t - oldt); final double rLimit = @@ -129,7 +124,6 @@ } } oldt = t; - i++; } throw new MaxIterationsExceededException(maximalIterationCount); } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/DividedDifferenceInterpolator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/DividedDifferenceInterpolator.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/DividedDifferenceInterpolator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/DividedDifferenceInterpolator.java Fri Oct 16 14:51:55 2009 @@ -57,8 +57,6 @@ * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + * a[n](x-c[0])(x-c[1])...(x-c[n-1]) */ - double a[], c[]; - PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y); /** @@ -69,12 +67,10 @@ *

* Note x[], y[], a[] have the same length but c[]'s size is one less.

*/ - c = new double[x.length-1]; - for (int i = 0; i < c.length; i++) { - c[i] = x[i]; - } - a = computeDividedDifference(x, y); + final double[] c = new double[x.length-1]; + System.arraycopy(x, 0, c, 0, c.length); + final double[] a = computeDividedDifference(x, y); return new PolynomialFunctionNewtonForm(a, c); } @@ -94,25 +90,19 @@ * @return a fresh copy of the divided difference array * @throws DuplicateSampleAbscissaException if any abscissas coincide */ - protected static double[] computeDividedDifference(double x[], double y[]) + protected static double[] computeDividedDifference(final double x[], final double y[]) throws DuplicateSampleAbscissaException { - int i, j, n; - double divdiff[], a[], denominator; - PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y); - n = x.length; - divdiff = new double[n]; - for (i = 0; i < n; i++) { - divdiff[i] = y[i]; // initialization - } + final double[] divdiff = y.clone(); // initialization - a = new double [n]; + final int n = x.length; + final double[] a = new double [n]; a[0] = divdiff[0]; - for (i = 1; i < n; i++) { - for (j = 0; j < n-i; j++) { - denominator = x[j+i] - x[j]; + for (int i = 1; i < n; i++) { + for (int j = 0; j < n-i; j++) { + final double denominator = x[j+i] - x[j]; if (denominator == 0.0) { // This happens only when two abscissas are identical. throw new DuplicateSampleAbscissaException(x[j], j, j+i); 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=825919&r1=825918&r2=825919&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 Fri Oct 16 14:51:55 2009 @@ -277,7 +277,10 @@ // and http://en.wikipedia.org/wiki/Weighted_least_squares // (section "Weighted least squares") double sumWeights = 0; - double sumX = 0, sumXSquared = 0, sumY = 0, sumXY = 0; + double sumX = 0; + double sumXSquared = 0; + double sumY = 0; + double sumXY = 0; double denom = Math.abs(1.0 / (xval[edge] - x)); for (int k = ileft; k <= iright; ++k) { final double xk = xval[k]; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java Fri Oct 16 14:51:55 2009 @@ -162,12 +162,13 @@ // Copy data samples. samples = new HashMap(yval.length); - for (int i = 0, max = xval.length; i < max; i++) { - if (xval[i].length != dimension) { - throw new DimensionMismatchException(xval.length, yval.length); + for (int i = 0; i < xval.length; ++i) { + final double[] xvalI = xval[i]; + if ( xvalI.length != dimension) { + throw new DimensionMismatchException(xvalI.length, dimension); } - samples.put(new ArrayRealVector(xval[i]), yval[i]); + samples.put(new ArrayRealVector(xvalI), yval[i]); } microsphere = new ArrayList(microsphereElements); 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=825919&r1=825918&r2=825919&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 Fri Oct 16 14:51:55 2009 @@ -43,9 +43,14 @@ private double coefficients[]; /** - * Interpolating points (abscissas) and the function values at these points. + * Interpolating points (abscissas). */ - private double x[], y[]; + private double x[]; + + /** + * Function values at interpolating points. + */ + private double y[]; /** * Whether the polynomial coefficients are available. @@ -158,21 +163,19 @@ public static double evaluate(double x[], double y[], double z) throws DuplicateSampleAbscissaException, IllegalArgumentException { - int i, j, n, nearest = 0; - double value, c[], d[], tc, td, divider, w, dist, min_dist; - verifyInterpolationArray(x, y); - n = x.length; - c = new double[n]; - d = new double[n]; - min_dist = Double.POSITIVE_INFINITY; - for (i = 0; i < n; i++) { + int nearest = 0; + final int n = x.length; + final double[] c = new double[n]; + final double[] d = new double[n]; + double min_dist = Double.POSITIVE_INFINITY; + for (int i = 0; i < n; i++) { // initialize the difference arrays c[i] = y[i]; d[i] = y[i]; // find out the abscissa closest to z - dist = Math.abs(z - x[i]); + final double dist = Math.abs(z - x[i]); if (dist < min_dist) { nearest = i; min_dist = dist; @@ -180,19 +183,19 @@ } // initial approximation to the function value at z - value = y[nearest]; + double value = y[nearest]; - for (i = 1; i < n; i++) { - for (j = 0; j < n-i; j++) { - tc = x[j] - z; - td = x[i+j] - z; - divider = x[j] - x[i+j]; + for (int i = 1; i < n; i++) { + for (int j = 0; j < n-i; j++) { + final double tc = x[j] - z; + final double td = x[i+j] - z; + final double divider = x[j] - x[i+j]; if (divider == 0.0) { // This happens only when two abscissas are identical. throw new DuplicateSampleAbscissaException(x[i], i, i+j); } // update the difference arrays - w = (c[j+1] - d[j]) / divider; + final double w = (c[j+1] - d[j]) / divider; c[j] = tc * w; d[j] = td * w; } @@ -218,31 +221,29 @@ * @throws ArithmeticException if any abscissas coincide */ protected void computeCoefficients() throws ArithmeticException { - int i, j, n; - double c[], tc[], d, t; - n = degree() + 1; + final int n = degree() + 1; coefficients = new double[n]; - for (i = 0; i < n; i++) { + for (int i = 0; i < n; i++) { coefficients[i] = 0.0; } // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1]) - c = new double[n+1]; + final double[] c = new double[n+1]; c[0] = 1.0; - for (i = 0; i < n; i++) { - for (j = i; j > 0; j--) { + for (int i = 0; i < n; i++) { + for (int j = i; j > 0; j--) { c[j] = c[j-1] - c[j] * x[i]; } c[0] *= -x[i]; c[i+1] = 1; } - tc = new double[n]; - for (i = 0; i < n; i++) { + final double[] tc = new double[n]; + for (int i = 0; i < n; i++) { // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1]) - d = 1; - for (j = 0; j < n; j++) { + double d = 1; + for (int j = 0; j < n; j++) { if (i != j) { d *= x[i] - x[j]; } @@ -256,13 +257,13 @@ } } } - t = y[i] / d; + final double t = y[i] / d; // Lagrange polynomial is the sum of n terms, each of which is a // polynomial of degree n-1. tc[] are the coefficients of the i-th // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]). tc[n-1] = c[n]; // actually c[n] = 1 coefficients[n-1] += t * tc[n-1]; - for (j = n-2; j >= 0; j--) { + for (int j = n-2; j >= 0; j--) { tc[j] = c[j+1] + tc[j+1] * x[i]; coefficients[j] += t * tc[j]; } 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=825919&r1=825918&r2=825919&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 Fri Oct 16 14:51:55 2009 @@ -43,11 +43,15 @@ private double coefficients[]; /** - * Members of c[] are called centers of the Newton polynomial. + * Centers of the Newton polynomial. + */ + private double c[]; + + /** * When all c[i] = 0, a[] becomes normal polynomial coefficients, * i.e. a[i] = coefficients[i]. */ - private double a[], c[]; + private double a[]; /** * Whether the polynomial coefficients are available. @@ -170,16 +174,16 @@ * It also uses nested multiplication but takes O(N^2) time. */ protected void computeCoefficients() { - int i, j, n = degree(); + final int n = degree(); coefficients = new double[n+1]; - for (i = 0; i <= n; i++) { + for (int i = 0; i <= n; i++) { coefficients[i] = 0.0; } coefficients[0] = a[n]; - for (i = n-1; i >= 0; i--) { - for (j = n-i; j > 0; j--) { + for (int i = n-1; i >= 0; i--) { + for (int j = n-i; j > 0; j--) { coefficients[j] = coefficients[j-1] - c[i] * coefficients[j]; } coefficients[0] = a[i] - c[i] * coefficients[0]; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver.java Fri Oct 16 14:51:55 2009 @@ -139,39 +139,43 @@ // x1 is the last approximation and an interpolation point in (x0, x2) // x is the new root approximation and new x1 for next round // d01, d12, d012 are divided differences - double x0, x1, x2, x, oldx, y0, y1, y2, y; - double d01, d12, d012, c1, delta, xplus, xminus, tolerance; - x0 = min; y0 = f.value(x0); - x2 = max; y2 = f.value(x2); - x1 = 0.5 * (x0 + x2); y1 = f.value(x1); + double x0 = min; + double y0 = f.value(x0); + double x2 = max; + double y2 = f.value(x2); + double x1 = 0.5 * (x0 + x2); + double y1 = f.value(x1); // check for zeros before verifying bracketing - if (y0 == 0.0) { return min; } - if (y2 == 0.0) { return max; } + if (y0 == 0.0) { + return min; + } + if (y2 == 0.0) { + return max; + } verifyBracketing(min, max, f); - int i = 1; - oldx = Double.POSITIVE_INFINITY; - while (i <= maximalIterationCount) { + double oldx = Double.POSITIVE_INFINITY; + for (int i = 1; i <= maximalIterationCount; ++i) { // Muller's method employs quadratic interpolation through // x0, x1, x2 and x is the zero of the interpolating parabola. // Due to bracketing condition, this parabola must have two // real roots and we choose one in [x0, x2] to be x. - d01 = (y1 - y0) / (x1 - x0); - d12 = (y2 - y1) / (x2 - x1); - d012 = (d12 - d01) / (x2 - x0); - c1 = d01 + (x1 - x0) * d012; - delta = c1 * c1 - 4 * y1 * d012; - xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta)); - xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta)); + final double d01 = (y1 - y0) / (x1 - x0); + final double d12 = (y2 - y1) / (x2 - x1); + final double d012 = (d12 - d01) / (x2 - x0); + final double c1 = d01 + (x1 - x0) * d012; + final double delta = c1 * c1 - 4 * y1 * d012; + final double xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta)); + final double xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta)); // xplus and xminus are two roots of parabola and at least // one of them should lie in (x0, x2) - x = isSequence(x0, xplus, x2) ? xplus : xminus; - y = f.value(x); + final double x = isSequence(x0, xplus, x2) ? xplus : xminus; + final double y = f.value(x); // check for convergence - tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy); + final double tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy); if (Math.abs(x - oldx) <= tolerance) { setResult(x, i); return result; @@ -208,7 +212,6 @@ y1 = f.value(x1); oldx = Double.POSITIVE_INFINITY; } - i++; } throw new MaxIterationsExceededException(maximalIterationCount); } @@ -279,38 +282,40 @@ // x2 is the last root approximation // x is the new approximation and new x2 for next round // x0 < x1 < x2 does not hold here - double x0, x1, x2, x, oldx, y0, y1, y2, y; - double q, A, B, C, delta, denominator, tolerance; - x0 = min; y0 = f.value(x0); - x1 = max; y1 = f.value(x1); - x2 = 0.5 * (x0 + x1); y2 = f.value(x2); + double x0 = min; + double y0 = f.value(x0); + double x1 = max; + double y1 = f.value(x1); + double x2 = 0.5 * (x0 + x1); + double y2 = f.value(x2); // check for zeros before verifying bracketing if (y0 == 0.0) { return min; } if (y1 == 0.0) { return max; } verifyBracketing(min, max, f); - int i = 1; - oldx = Double.POSITIVE_INFINITY; - while (i <= maximalIterationCount) { + double oldx = Double.POSITIVE_INFINITY; + for (int i = 1; i <= maximalIterationCount; ++i) { // quadratic interpolation through x0, x1, x2 - q = (x2 - x1) / (x1 - x0); - A = q * (y2 - (1 + q) * y1 + q * y0); - B = (2*q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0; - C = (1 + q) * y2; - delta = B * B - 4 * A * C; + final double q = (x2 - x1) / (x1 - x0); + final double a = q * (y2 - (1 + q) * y1 + q * y0); + final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0; + final double c = (1 + q) * y2; + final double delta = b * b - 4 * a * c; + double x; + final double denominator; if (delta >= 0.0) { // choose a denominator larger in magnitude - double dplus = B + Math.sqrt(delta); - double dminus = B - Math.sqrt(delta); + double dplus = b + Math.sqrt(delta); + double dminus = b - Math.sqrt(delta); denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus; } else { // take the modulus of (B +/- Math.sqrt(delta)) - denominator = Math.sqrt(B * B - delta); + denominator = Math.sqrt(b * b - delta); } if (denominator != 0) { - x = x2 - 2.0 * C * (x2 - x1) / denominator; + x = x2 - 2.0 * c * (x2 - x1) / denominator; // perturb x if it exactly coincides with x1 or x2 // the equality tests here are intentional while (x == x1 || x == x2) { @@ -321,10 +326,10 @@ x = min + Math.random() * (max - min); oldx = Double.POSITIVE_INFINITY; } - y = f.value(x); + final double y = f.value(x); // check for convergence - tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy); + final double tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy); if (Math.abs(x - oldx) <= tolerance) { setResult(x, i); return result; @@ -335,11 +340,13 @@ } // prepare the next iteration - x0 = x1; y0 = y1; - x1 = x2; y1 = y2; - x2 = x; y2 = y; + x0 = x1; + y0 = y1; + x1 = x2; + y1 = y2; + x2 = x; + y2 = y; oldx = x; - i++; } throw new MaxIterationsExceededException(maximalIterationCount); } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/RiddersSolver.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/RiddersSolver.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/RiddersSolver.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/solvers/RiddersSolver.java Fri Oct 16 14:51:55 2009 @@ -125,34 +125,38 @@ // [x1, x2] is the bracketing interval in each iteration // x3 is the midpoint of [x1, x2] // x is the new root approximation and an endpoint of the new interval - double x1, x2, x3, x, oldx, y1, y2, y3, y, delta, correction, tolerance; - - x1 = min; y1 = f.value(x1); - x2 = max; y2 = f.value(x2); + double x1 = min; + double y1 = f.value(x1); + double x2 = max; + double y2 = f.value(x2); // check for zeros before verifying bracketing - if (y1 == 0.0) { return min; } - if (y2 == 0.0) { return max; } + if (y1 == 0.0) { + return min; + } + if (y2 == 0.0) { + return max; + } verifyBracketing(min, max, f); int i = 1; - oldx = Double.POSITIVE_INFINITY; + double oldx = Double.POSITIVE_INFINITY; while (i <= maximalIterationCount) { // calculate the new root approximation - x3 = 0.5 * (x1 + x2); - y3 = f.value(x3); + final double x3 = 0.5 * (x1 + x2); + final double y3 = f.value(x3); if (Math.abs(y3) <= functionValueAccuracy) { setResult(x3, i); return result; } - delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing - correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) * - (x3 - x1) / Math.sqrt(delta); - x = x3 - correction; // correction != 0 - y = f.value(x); + final double delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing + final double correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) * + (x3 - x1) / Math.sqrt(delta); + final double x = x3 - correction; // correction != 0 + final double y = f.value(x); // check for convergence - tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy); + final double tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy); if (Math.abs(x - oldx) <= tolerance) { setResult(x, i); return result; @@ -166,17 +170,23 @@ // Ridders' method guarantees x1 < x < x2 if (correction > 0.0) { // x1 < x < x3 if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) { - x2 = x; y2 = y; + x2 = x; + y2 = y; } else { - x1 = x; x2 = x3; - y1 = y; y2 = y3; + x1 = x; + x2 = x3; + y1 = y; + y2 = y3; } } else { // x3 < x < x2 if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) { - x1 = x; y1 = y; + x1 = x; + y1 = y; } else { - x1 = x3; x2 = x; - y1 = y3; y2 = y; + x1 = x3; + x2 = x; + y1 = y3; + y2 = y; } } oldx = x; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java Fri Oct 16 14:51:55 2009 @@ -124,7 +124,8 @@ protected void updateJacobian() { incrementJacobianEvaluationsCounter(); Arrays.fill(jacobian, 0); - for (int i = 0, index = 0; i < rows; i++) { + int index = 0; + for (int i = 0; i < rows; i++) { WeightedMeasurement wm = measurements[i]; double factor = -Math.sqrt(wm.getWeight()); for (int j = 0; j < cols; ++j) { @@ -154,7 +155,8 @@ } cost = 0; - for (int i = 0, index = 0; i < rows; i++, index += cols) { + int index = 0; + for (int i = 0; i < rows; i++, index += cols) { WeightedMeasurement wm = measurements[i]; double residual = wm.getResidual(); residuals[i] = Math.sqrt(wm.getWeight()) * residual; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java Fri Oct 16 14:51:55 2009 @@ -256,7 +256,8 @@ lmDir = new double[cols]; // local variables - double delta = 0, xNorm = 0; + double delta = 0; + double xNorm = 0; double[] diag = new double[cols]; double[] oldX = new double[cols]; double[] oldRes = new double[rows]; @@ -315,8 +316,10 @@ double s = jacNorm[pj]; if (s != 0) { double sum = 0; - for (int i = 0, index = pj; i <= j; ++i, index += cols) { + int index = pj; + for (int i = 0; i <= j; ++i) { sum += jacobian[index] * residuals[i]; + index += cols; } maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost)); } @@ -379,8 +382,10 @@ int pj = permutation[j]; double dirJ = lmDir[pj]; work1[j] = 0; - for (int i = 0, index = pj; i <= j; ++i, index += cols) { + int index = pj; + for (int i = 0; i <= j; ++i) { work1[i] += jacobian[index] * dirJ; + index += cols; } } double coeff1 = 0; @@ -500,8 +505,10 @@ for (int k = rank - 1; k >= 0; --k) { int pk = permutation[k]; double ypk = lmDir[pk] / diagR[pk]; - for (int i = 0, index = pk; i < k; ++i, index += cols) { + int index = pk; + for (int i = 0; i < k; ++i) { lmDir[permutation[i]] -= ypk * jacobian[index]; + index += cols; } lmDir[pk] = ypk; } @@ -525,7 +532,8 @@ // if the jacobian is not rank deficient, the Newton step provides // a lower bound, parl, for the zero of the function, // otherwise set this bound to zero - double sum2, parl = 0; + double sum2; + double parl = 0; if (rank == solvedCols) { for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; @@ -535,8 +543,10 @@ for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double sum = 0; - for (int i = 0, index = pj; i < j; ++i, index += cols) { + int index = pj; + for (int i = 0; i < j; ++i) { sum += jacobian[index] * work1[permutation[i]]; + index += cols; } double s = (work1[pj] - sum) / diagR[pj]; work1[pj] = s; @@ -550,8 +560,10 @@ for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double sum = 0; - for (int i = 0, index = pj; i <= j; ++i, index += cols) { + int index = pj; + for (int i = 0; i <= j; ++i) { sum += jacobian[index] * qy[i]; + index += cols; } sum /= diag[pj]; sum2 += sum * sum; @@ -691,14 +703,15 @@ // appropriate element in the current row of d if (lmDiag[k] != 0) { - double sin, cos; + final double sin; + final double cos; double rkk = jacobian[k * cols + pk]; if (Math.abs(rkk) < Math.abs(lmDiag[k])) { - double cotan = rkk / lmDiag[k]; + final double cotan = rkk / lmDiag[k]; sin = 1.0 / Math.sqrt(1.0 + cotan * cotan); cos = sin * cotan; } else { - double tan = lmDiag[k] / rkk; + final double tan = lmDiag[k] / rkk; cos = 1.0 / Math.sqrt(1.0 + tan * tan); sin = cos * tan; } @@ -706,16 +719,16 @@ // compute the modified diagonal element of R and // the modified element of (Qty,0) jacobian[k * cols + pk] = cos * rkk + sin * lmDiag[k]; - double temp = cos * work[k] + sin * qtbpj; + final double temp = cos * work[k] + sin * qtbpj; qtbpj = -sin * work[k] + cos * qtbpj; work[k] = temp; // accumulate the tranformation in the row of s for (int i = k + 1; i < solvedCols; ++i) { double rik = jacobian[i * cols + pk]; - temp = cos * rik + sin * lmDiag[i]; + final double temp2 = cos * rik + sin * lmDiag[i]; lmDiag[i] = -sin * rik + cos * lmDiag[i]; - jacobian[i * cols + pk] = temp; + jacobian[i * cols + pk] = temp2; } } @@ -864,12 +877,16 @@ int pk = permutation[k]; int kDiag = k * cols + pk; double gamma = 0; - for (int i = k, index = kDiag; i < rows; ++i, index += cols) { + int index = kDiag; + for (int i = k; i < rows; ++i) { gamma += jacobian[index] * y[i]; + index += cols; } gamma *= beta[pk]; - for (int i = k, index = kDiag; i < rows; ++i, index += cols) { + index = kDiag; + for (int i = k; i < rows; ++i) { y[i] -= gamma * jacobian[index]; + index += cols; } } } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockFieldMatrix.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockFieldMatrix.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockFieldMatrix.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockFieldMatrix.java Fri Oct 16 14:51:55 2009 @@ -226,11 +226,12 @@ // convert array final Field field = extractField(rawData); final T[][] blocks = buildArray(field, blockRows * blockColumns, -1); - for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { final int pStart = iBlock * BLOCK_SIZE; final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); final int iHeight = pEnd - pStart; - for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final int jWidth = qEnd - qStart; @@ -240,10 +241,14 @@ blocks[blockIndex] = block; // copy data - for (int p = pStart, index = 0; p < pEnd; ++p, index += jWidth) { + int index = 0; + for (int p = pStart; p < pEnd; ++p) { System.arraycopy(rawData[p], qStart, block, index, jWidth); + index += jWidth; } + ++blockIndex; + } } @@ -273,15 +278,17 @@ final int blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE; final T[][] blocks = buildArray(field, blockRows * blockColumns, -1); - for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { final int pStart = iBlock * BLOCK_SIZE; final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); final int iHeight = pEnd - pStart; - for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final int jWidth = qEnd - qStart; blocks[blockIndex] = buildArray(field, iHeight * jWidth); + ++blockIndex; } } @@ -337,9 +344,11 @@ final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); - for (int p = pStart, k = 0; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q, ++k) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q) { outBlock[k] = tBlock[k].add(m.getEntry(p, q)); + ++k; } } @@ -408,9 +417,11 @@ final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); - for (int p = pStart, k = 0; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q, ++k) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q) { outBlock[k] = tBlock[k].subtract(m.getEntry(p, q)); + ++k; } } @@ -528,13 +539,16 @@ final int kWidth = blockWidth(kBlock); final T[] tBlock = blocks[iBlock * blockColumns + kBlock]; final int rStart = kBlock * BLOCK_SIZE; - for (int p = pStart, k = 0; p < pEnd; ++p) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { final int lStart = (p - pStart) * kWidth; final int lEnd = lStart + kWidth; for (int q = qStart; q < qEnd; ++q) { T sum = zero; - for (int l = lStart, r = rStart; l < lEnd; ++l, ++r) { + int r = rStart; + for (int l = lStart; l < lEnd; ++l) { sum = sum.add(tBlock[l].multiply(m.getEntry(r, q))); + ++r; } outBlock[k] = outBlock[k].add(sum); ++k; @@ -590,7 +604,8 @@ final int kWidth = blockWidth(kBlock); final T[] tBlock = blocks[iBlock * blockColumns + kBlock]; final T[] mBlock = m.blocks[kBlock * m.blockColumns + jBlock]; - for (int p = pStart, k = 0; p < pEnd; ++p) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { final int lStart = (p - pStart) * kWidth; final int lEnd = lStart + kWidth; for (int nStart = 0; nStart < jWidth; ++nStart) { @@ -676,9 +691,11 @@ final int columnsShift = startColumn % BLOCK_SIZE; // perform extraction block-wise, to ensure good cache behavior - for (int iBlock = 0, pBlock = blockStartRow; iBlock < out.blockRows; ++iBlock, ++pBlock) { + int pBlock = blockStartRow; + for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) { final int iHeight = out.blockHeight(iBlock); - for (int jBlock = 0, qBlock = blockStartColumn; jBlock < out.blockColumns; ++jBlock, ++qBlock) { + int qBlock = blockStartColumn; + for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) { final int jWidth = out.blockWidth(jBlock); // handle one block of the output matrix @@ -743,7 +760,11 @@ } } + ++qBlock; } + + ++pBlock; + } return out; @@ -1273,10 +1294,14 @@ final int pEnd = Math.min(pStart + BLOCK_SIZE, columns); final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, rows); - for (int p = pStart, k = 0; p < pEnd; ++p) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { final int lInc = pEnd - pStart; - for (int q = qStart, l = p - pStart; q < qEnd; ++q, l+= lInc) { - outBlock[k++] = tBlock[l]; + int l = p - pStart; + for (int q = qStart; q < qEnd; ++q) { + outBlock[k] = tBlock[l]; + ++k; + l+= lInc; } } @@ -1323,7 +1348,8 @@ final T[] block = blocks[iBlock * blockColumns + jBlock]; final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); - for (int p = pStart, k = 0; p < pEnd; ++p) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { T sum = zero; int q = qStart; while (q < qEnd - 3) { @@ -1412,8 +1438,10 @@ final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final T[] block = blocks[iBlock * blockColumns + jBlock]; - for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) { + int k = (p - pStart) * jWidth; + for (int q = qStart; q < qEnd; ++q) { block[k] = visitor.visit(p, q, block[k]); + ++k; } } } @@ -1435,8 +1463,10 @@ final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final T[] block = blocks[iBlock * blockColumns + jBlock]; - for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) { + int k = (p - pStart) * jWidth; + for (int q = qStart; q < qEnd; ++q) { visitor.visit(p, q, block[k]); + ++k; } } } @@ -1463,8 +1493,10 @@ final int qStart = Math.max(startColumn, q0); final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); final T[] block = blocks[iBlock * blockColumns + jBlock]; - for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + int k = (p - p0) * jWidth + qStart - q0; + for (int q = qStart; q < qEnd; ++q) { block[k] = visitor.visit(p, q, block[k]); + ++k; } } } @@ -1491,8 +1523,10 @@ final int qStart = Math.max(startColumn, q0); final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); final T[] block = blocks[iBlock * blockColumns + jBlock]; - for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + int k = (p - p0) * jWidth + qStart - q0; + for (int q = qStart; q < qEnd; ++q) { visitor.visit(p, q, block[k]); + ++k; } } } @@ -1505,18 +1539,22 @@ public T walkInOptimizedOrder(final FieldMatrixChangingVisitor visitor) throws MatrixVisitorException { visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); - for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { final int pStart = iBlock * BLOCK_SIZE; final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); - for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final T[] block = blocks[blockIndex]; - for (int p = pStart, k = 0; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q, ++k) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q) { block[k] = visitor.visit(p, q, block[k]); + ++k; } } + ++blockIndex; } } return visitor.end(); @@ -1527,18 +1565,22 @@ public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor visitor) throws MatrixVisitorException { visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); - for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { final int pStart = iBlock * BLOCK_SIZE; final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); - for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final T[] block = blocks[blockIndex]; - for (int p = pStart, k = 0; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q, ++k) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q) { visitor.visit(p, q, block[k]); + ++k; } } + ++blockIndex; } } return visitor.end(); @@ -1563,8 +1605,10 @@ final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); final T[] block = blocks[iBlock * blockColumns + jBlock]; for (int p = pStart; p < pEnd; ++p) { - for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + int k = (p - p0) * jWidth + qStart - q0; + for (int q = qStart; q < qEnd; ++q) { block[k] = visitor.visit(p, q, block[k]); + ++k; } } } @@ -1591,8 +1635,10 @@ final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); final T[] block = blocks[iBlock * blockColumns + jBlock]; for (int p = pStart; p < pEnd; ++p) { - for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + int k = (p - p0) * jWidth + qStart - q0; + for (int q = qStart; q < qEnd; ++q) { visitor.visit(p, q, block[k]); + ++k; } } } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockRealMatrix.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockRealMatrix.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockRealMatrix.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/BlockRealMatrix.java Fri Oct 16 14:51:55 2009 @@ -220,11 +220,12 @@ // convert array final double[][] blocks = new double[blockRows * blockColumns][]; - for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { final int pStart = iBlock * BLOCK_SIZE; final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); final int iHeight = pEnd - pStart; - for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final int jWidth = qEnd - qStart; @@ -234,10 +235,14 @@ blocks[blockIndex] = block; // copy data - for (int p = pStart, index = 0; p < pEnd; ++p, index += jWidth) { + int index = 0; + for (int p = pStart; p < pEnd; ++p) { System.arraycopy(rawData[p], qStart, block, index, jWidth); + index += jWidth; } + ++blockIndex; + } } @@ -263,15 +268,17 @@ final int blockColumns = (columns + BLOCK_SIZE - 1) / BLOCK_SIZE; final double[][] blocks = new double[blockRows * blockColumns][]; - for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { final int pStart = iBlock * BLOCK_SIZE; final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); final int iHeight = pEnd - pStart; - for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final int jWidth = qEnd - qStart; blocks[blockIndex] = new double[iHeight * jWidth]; + ++blockIndex; } } @@ -327,9 +334,11 @@ final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); - for (int p = pStart, k = 0; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q, ++k) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q) { outBlock[k] = tBlock[k] + m.getEntry(p, q); + ++k; } } @@ -398,9 +407,11 @@ final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); - for (int p = pStart, k = 0; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q, ++k) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q) { outBlock[k] = tBlock[k] - m.getEntry(p, q); + ++k; } } @@ -517,15 +528,19 @@ final int kWidth = blockWidth(kBlock); final double[] tBlock = blocks[iBlock * blockColumns + kBlock]; final int rStart = kBlock * BLOCK_SIZE; - for (int p = pStart, k = 0; p < pEnd; ++p) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { final int lStart = (p - pStart) * kWidth; final int lEnd = lStart + kWidth; for (int q = qStart; q < qEnd; ++q) { double sum = 0; - for (int l = lStart, r = rStart; l < lEnd; ++l, ++r) { + int r = rStart; + for (int l = lStart; l < lEnd; ++l) { sum += tBlock[l] * m.getEntry(r, q); + ++r; } - outBlock[k++] += sum; + outBlock[k] += sum; + ++k; } } } @@ -577,7 +592,8 @@ final int kWidth = blockWidth(kBlock); final double[] tBlock = blocks[iBlock * blockColumns + kBlock]; final double[] mBlock = m.blocks[kBlock * m.blockColumns + jBlock]; - for (int p = pStart, k = 0; p < pEnd; ++p) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { final int lStart = (p - pStart) * kWidth; final int lEnd = lStart + kWidth; for (int nStart = 0; nStart < jWidth; ++nStart) { @@ -596,7 +612,8 @@ sum += tBlock[l++] * mBlock[n]; n += jWidth; } - outBlock[k++] += sum; + outBlock[k] += sum; + ++k; } } } @@ -699,9 +716,11 @@ final int columnsShift = startColumn % BLOCK_SIZE; // perform extraction block-wise, to ensure good cache behavior - for (int iBlock = 0, pBlock = blockStartRow; iBlock < out.blockRows; ++iBlock, ++pBlock) { + int pBlock = blockStartRow; + for (int iBlock = 0; iBlock < out.blockRows; ++iBlock) { final int iHeight = out.blockHeight(iBlock); - for (int jBlock = 0, qBlock = blockStartColumn; jBlock < out.blockColumns; ++jBlock, ++qBlock) { + int qBlock = blockStartColumn; + for (int jBlock = 0; jBlock < out.blockColumns; ++jBlock) { final int jWidth = out.blockWidth(jBlock); // handle one block of the output matrix @@ -766,7 +785,12 @@ } } + ++qBlock; + } + + ++pBlock; + } return out; @@ -1294,10 +1318,14 @@ final int pEnd = Math.min(pStart + BLOCK_SIZE, columns); final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, rows); - for (int p = pStart, k = 0; p < pEnd; ++p) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { final int lInc = pEnd - pStart; - for (int q = qStart, l = p - pStart; q < qEnd; ++q, l+= lInc) { - outBlock[k++] = tBlock[l]; + int l = p - pStart; + for (int q = qStart; q < qEnd; ++q) { + outBlock[k] = tBlock[l]; + ++k; + l+= lInc; } } @@ -1343,7 +1371,8 @@ final double[] block = blocks[iBlock * blockColumns + jBlock]; final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); - for (int p = pStart, k = 0; p < pEnd; ++p) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { double sum = 0; int q = qStart; while (q < qEnd - 3) { @@ -1429,8 +1458,10 @@ final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final double[] block = blocks[iBlock * blockColumns + jBlock]; - for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) { + int k = (p - pStart) * jWidth; + for (int q = qStart; q < qEnd; ++q) { block[k] = visitor.visit(p, q, block[k]); + ++k; } } } @@ -1452,8 +1483,10 @@ final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final double[] block = blocks[iBlock * blockColumns + jBlock]; - for (int q = qStart, k = (p - pStart) * jWidth; q < qEnd; ++q, ++k) { + int k = (p - pStart) * jWidth; + for (int q = qStart; q < qEnd; ++q) { visitor.visit(p, q, block[k]); + ++k; } } } @@ -1480,8 +1513,10 @@ final int qStart = Math.max(startColumn, q0); final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); final double[] block = blocks[iBlock * blockColumns + jBlock]; - for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + int k = (p - p0) * jWidth + qStart - q0; + for (int q = qStart; q < qEnd; ++q) { block[k] = visitor.visit(p, q, block[k]); + ++k; } } } @@ -1508,8 +1543,10 @@ final int qStart = Math.max(startColumn, q0); final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); final double[] block = blocks[iBlock * blockColumns + jBlock]; - for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + int k = (p - p0) * jWidth + qStart - q0; + for (int q = qStart; q < qEnd; ++q) { visitor.visit(p, q, block[k]); + ++k; } } } @@ -1522,18 +1559,22 @@ public double walkInOptimizedOrder(final RealMatrixChangingVisitor visitor) throws MatrixVisitorException { visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); - for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { final int pStart = iBlock * BLOCK_SIZE; final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); - for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final double[] block = blocks[blockIndex]; - for (int p = pStart, k = 0; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q, ++k) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q) { block[k] = visitor.visit(p, q, block[k]); + ++k; } } + ++blockIndex; } } return visitor.end(); @@ -1544,18 +1585,22 @@ public double walkInOptimizedOrder(final RealMatrixPreservingVisitor visitor) throws MatrixVisitorException { visitor.start(rows, columns, 0, rows - 1, 0, columns - 1); - for (int iBlock = 0, blockIndex = 0; iBlock < blockRows; ++iBlock) { + int blockIndex = 0; + for (int iBlock = 0; iBlock < blockRows; ++iBlock) { final int pStart = iBlock * BLOCK_SIZE; final int pEnd = Math.min(pStart + BLOCK_SIZE, rows); - for (int jBlock = 0; jBlock < blockColumns; ++jBlock, ++blockIndex) { + for (int jBlock = 0; jBlock < blockColumns; ++jBlock) { final int qStart = jBlock * BLOCK_SIZE; final int qEnd = Math.min(qStart + BLOCK_SIZE, columns); final double[] block = blocks[blockIndex]; - for (int p = pStart, k = 0; p < pEnd; ++p) { - for (int q = qStart; q < qEnd; ++q, ++k) { + int k = 0; + for (int p = pStart; p < pEnd; ++p) { + for (int q = qStart; q < qEnd; ++q) { visitor.visit(p, q, block[k]); + ++k; } } + ++blockIndex; } } return visitor.end(); @@ -1580,8 +1625,10 @@ final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); final double[] block = blocks[iBlock * blockColumns + jBlock]; for (int p = pStart; p < pEnd; ++p) { - for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + int k = (p - p0) * jWidth + qStart - q0; + for (int q = qStart; q < qEnd; ++q) { block[k] = visitor.visit(p, q, block[k]); + ++k; } } } @@ -1608,8 +1655,10 @@ final int qEnd = Math.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn); final double[] block = blocks[iBlock * blockColumns + jBlock]; for (int p = pStart; p < pEnd; ++p) { - for (int q = qStart, k = (p - p0) * jWidth + qStart - q0; q < qEnd; ++q, ++k) { + int k = (p - p0) * jWidth + qStart - q0; + for (int q = qStart; q < qEnd; ++q) { visitor.visit(p, q, block[k]); + ++k; } } } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java Fri Oct 16 14:51:55 2009 @@ -670,10 +670,12 @@ // sort the realEigenvalues in decreasing order Arrays.sort(realEigenvalues); - for (int i = 0, j = realEigenvalues.length - 1; i < j; ++i, --j) { + int j = realEigenvalues.length - 1; + for (int i = 0; i < j; ++i) { final double tmp = realEigenvalues[i]; realEigenvalues[i] = realEigenvalues[j]; realEigenvalues[j] = tmp; + --j; } } @@ -1126,12 +1128,14 @@ private boolean flipIfWarranted(final int n, final int step) { if (1.5 * work[pingPong] < work[4 * (n - 1) + pingPong]) { // flip array - for (int i = 0, j = 4 * n - 1; i < j; i += 4, j -= 4) { + int j = 4 * n - 1; + for (int i = 0; i < j; i += 4) { for (int k = 0; k < 4; k += step) { final double tmp = work[i + k]; work[i + k] = work[j - k]; work[j - k] = tmp; } + j -= 4; } return true; } @@ -1734,12 +1738,14 @@ // the least diagonal element in the twisted factorization int r = m - 1; double minG = Math.abs(work[6 * r] + work[6 * r + 3] + eigenvalue); - for (int i = 0, sixI = 0; i < m - 1; ++i, sixI += 6) { + int sixI = 0; + for (int i = 0; i < m - 1; ++i) { final double absG = Math.abs(work[sixI] + d[i] * work[sixI + 9] / work[sixI + 10]); if (absG < minG) { r = i; minG = absG; } + sixI += 6; } // solve the singular system by ignoring the equation @@ -1784,7 +1790,8 @@ final double lambda) { final int nM1 = d.length - 1; double si = -lambda; - for (int i = 0, sixI = 0; i < nM1; ++i, sixI += 6) { + int sixI = 0; + for (int i = 0; i < nM1; ++i) { final double di = d[i]; final double li = l[i]; final double diP1 = di + si; @@ -1793,6 +1800,7 @@ work[sixI + 1] = diP1; work[sixI + 2] = liP1; si = li * liP1 * si - lambda; + sixI += 6; } work[6 * nM1 + 1] = d[nM1] + si; work[6 * nM1] = si; @@ -1810,7 +1818,8 @@ final double lambda) { final int nM1 = d.length - 1; double pi = d[nM1] - lambda; - for (int i = nM1 - 1, sixI = 6 * i; i >= 0; --i, sixI -= 6) { + int sixI = 6 * (nM1 - 1); + for (int i = nM1 - 1; i >= 0; --i) { final double di = d[i]; final double li = l[i]; final double diP1 = di * li * li + pi; @@ -1819,6 +1828,7 @@ work[sixI + 10] = diP1; work[sixI + 5] = li * t; pi = pi * t - lambda; + sixI -= 6; } work[3] = pi; work[4] = pi; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java Fri Oct 16 14:51:55 2009 @@ -417,9 +417,10 @@ final double factor = 1.0 / rDiag[j]; final double[] yJ = y[j]; final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock]; - for (int k = 0, index = (j - jStart) * kWidth; k < kWidth; ++k, ++index) { - yJ[k] *= factor; - xBlock[index] = yJ[k]; + int index = (j - jStart) * kWidth; + for (int k = 0; k < kWidth; ++k) { + yJ[k] *= factor; + xBlock[index++] = yJ[k]; } final double[] qrtJ = qrt[j]; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java Fri Oct 16 14:51:55 2009 @@ -208,10 +208,12 @@ objective.length, rows); } cost = 0; - for (int i = 0, index = 0; i < rows; i++, index += cols) { + int index = 0; + for (int i = 0; i < rows; i++) { final double residual = targetValues[i] - objective[i]; residuals[i] = residual; cost += residualsWeights[i] * residual * residual; + index += cols; } cost = Math.sqrt(cost); Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java Fri Oct 16 14:51:55 2009 @@ -220,7 +220,8 @@ lmDir = new double[cols]; // local point - double delta = 0, xNorm = 0; + double delta = 0; + double xNorm = 0; double[] diag = new double[cols]; double[] oldX = new double[cols]; double[] oldRes = new double[rows]; @@ -492,7 +493,8 @@ // if the jacobian is not rank deficient, the Newton step provides // a lower bound, parl, for the zero of the function, // otherwise set this bound to zero - double sum2, parl = 0; + double sum2; + double parl = 0; if (rank == solvedCols) { for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; @@ -658,14 +660,15 @@ // appropriate element in the current row of d if (lmDiag[k] != 0) { - double sin, cos; + final double sin; + final double cos; double rkk = jacobian[k][pk]; if (Math.abs(rkk) < Math.abs(lmDiag[k])) { - double cotan = rkk / lmDiag[k]; + final double cotan = rkk / lmDiag[k]; sin = 1.0 / Math.sqrt(1.0 + cotan * cotan); cos = sin * cotan; } else { - double tan = lmDiag[k] / rkk; + final double tan = lmDiag[k] / rkk; cos = 1.0 / Math.sqrt(1.0 + tan * tan); sin = cos * tan; } @@ -673,16 +676,16 @@ // compute the modified diagonal element of R and // the modified element of (Qty,0) jacobian[k][pk] = cos * rkk + sin * lmDiag[k]; - double temp = cos * work[k] + sin * qtbpj; + final double temp = cos * work[k] + sin * qtbpj; qtbpj = -sin * work[k] + cos * qtbpj; work[k] = temp; // accumulate the tranformation in the row of s for (int i = k + 1; i < solvedCols; ++i) { double rik = jacobian[i][pk]; - temp = cos * rik + sin * lmDiag[i]; + final double temp2 = cos * rik + sin * lmDiag[i]; lmDiag[i] = -sin * rik + cos * lmDiag[i]; - jacobian[i][pk] = temp; + jacobian[i][pk] = temp2; } } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastCosineTransformer.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastCosineTransformer.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastCosineTransformer.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastCosineTransformer.java Fri Oct 16 14:51:55 2009 @@ -218,45 +218,45 @@ protected double[] fct(double f[]) throws IllegalArgumentException { - double A, B, C, F1, x[], F[] = new double[f.length]; + final double transformed[] = new double[f.length]; - int N = f.length - 1; - if (!FastFourierTransformer.isPowerOf2(N)) { + final int n = f.length - 1; + if (!FastFourierTransformer.isPowerOf2(n)) { throw MathRuntimeException.createIllegalArgumentException( "{0} is not a power of 2 plus one", f.length); } - if (N == 1) { // trivial case - F[0] = 0.5 * (f[0] + f[1]); - F[1] = 0.5 * (f[0] - f[1]); - return F; + if (n == 1) { // trivial case + transformed[0] = 0.5 * (f[0] + f[1]); + transformed[1] = 0.5 * (f[0] - f[1]); + return transformed; } // construct a new array and perform FFT on it - x = new double[N]; - x[0] = 0.5 * (f[0] + f[N]); - x[N >> 1] = f[N >> 1]; - F1 = 0.5 * (f[0] - f[N]); // temporary variable for F[1] - for (int i = 1; i < (N >> 1); i++) { - A = 0.5 * (f[i] + f[N-i]); - B = Math.sin(i * Math.PI / N) * (f[i] - f[N-i]); - C = Math.cos(i * Math.PI / N) * (f[i] - f[N-i]); - x[i] = A - B; - x[N-i] = A + B; - F1 += C; + final double[] x = new double[n]; + x[0] = 0.5 * (f[0] + f[n]); + x[n >> 1] = f[n >> 1]; + double t1 = 0.5 * (f[0] - f[n]); // temporary variable for transformed[1] + for (int i = 1; i < (n >> 1); i++) { + final double a = 0.5 * (f[i] + f[n-i]); + final double b = Math.sin(i * Math.PI / n) * (f[i] - f[n-i]); + final double c = Math.cos(i * Math.PI / n) * (f[i] - f[n-i]); + x[i] = a - b; + x[n-i] = a + b; + t1 += c; } FastFourierTransformer transformer = new FastFourierTransformer(); Complex y[] = transformer.transform(x); // reconstruct the FCT result for the original array - F[0] = y[0].getReal(); - F[1] = F1; - for (int i = 1; i < (N >> 1); i++) { - F[2*i] = y[i].getReal(); - F[2*i+1] = F[2*i-1] - y[i].getImaginary(); + transformed[0] = y[0].getReal(); + transformed[1] = t1; + for (int i = 1; i < (n >> 1); i++) { + transformed[2 * i] = y[i].getReal(); + transformed[2 * i + 1] = transformed[2 * i - 1] - y[i].getImaginary(); } - F[N] = y[N >> 1].getReal(); + transformed[n] = y[n >> 1].getReal(); - return F; + return transformed; } } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java Fri Oct 16 14:51:55 2009 @@ -344,58 +344,58 @@ protected Complex[] fft(Complex data[]) throws IllegalArgumentException { - int i, j, k, m, N = data.length; - Complex A, B, C, D, E, F, z, f[] = new Complex[N]; + final int n = data.length; + final Complex f[] = new Complex[n]; // initial simple cases verifyDataSet(data); - if (N == 1) { + if (n == 1) { f[0] = data[0]; return f; } - if (N == 2) { + if (n == 2) { f[0] = data[0].add(data[1]); f[1] = data[0].subtract(data[1]); return f; } // permute original data array in bit-reversal order - j = 0; - for (i = 0; i < N; i++) { - f[i] = data[j]; - k = N >> 1; - while (j >= k && k > 0) { - j -= k; k >>= 1; + int ii = 0; + for (int i = 0; i < n; i++) { + f[i] = data[ii]; + int k = n >> 1; + while (ii >= k && k > 0) { + ii -= k; k >>= 1; } - j += k; + ii += k; } // the bottom base-4 round - for (i = 0; i < N; i += 4) { - A = f[i].add(f[i+1]); - B = f[i+2].add(f[i+3]); - C = f[i].subtract(f[i+1]); - D = f[i+2].subtract(f[i+3]); - E = C.add(D.multiply(Complex.I)); - F = C.subtract(D.multiply(Complex.I)); - f[i] = A.add(B); - f[i+2] = A.subtract(B); + for (int i = 0; i < n; i += 4) { + final Complex a = f[i].add(f[i+1]); + final Complex b = f[i+2].add(f[i+3]); + final Complex c = f[i].subtract(f[i+1]); + final Complex d = f[i+2].subtract(f[i+3]); + final Complex e1 = c.add(d.multiply(Complex.I)); + final Complex e2 = c.subtract(d.multiply(Complex.I)); + f[i] = a.add(b); + f[i+2] = a.subtract(b); // omegaCount indicates forward or inverse transform - f[i+1] = roots.isForward() ? F : E; - f[i+3] = roots.isForward() ? E : F; + f[i+1] = roots.isForward() ? e2 : e1; + f[i+3] = roots.isForward() ? e1 : e2; } // iterations from bottom to top take O(N*logN) time - for (i = 4; i < N; i <<= 1) { - m = N / (i<<1); - for (j = 0; j < N; j += i<<1) { - for (k = 0; k < i; k++) { + for (int i = 4; i < n; i <<= 1) { + final int m = n / (i<<1); + for (int j = 0; j < n; j += i<<1) { + for (int k = 0; k < i; k++) { //z = f[i+j+k].multiply(roots.getOmega(k*m)); final int k_times_m = k*m; final double omega_k_times_m_real = roots.getOmegaReal(k_times_m); final double omega_k_times_m_imaginary = roots.getOmegaImaginary(k_times_m); //z = f[i+j+k].multiply(omega[k*m]); - z = new Complex( + final Complex z = new Complex( f[i+j+k].getReal() * omega_k_times_m_real - f[i+j+k].getImaginary() * omega_k_times_m_imaginary, f[i+j+k].getReal() * omega_k_times_m_imaginary + Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastSineTransformer.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastSineTransformer.java?rev=825919&r1=825918&r2=825919&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastSineTransformer.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/transform/FastSineTransformer.java Fri Oct 16 14:51:55 2009 @@ -212,7 +212,7 @@ */ protected double[] fst(double f[]) throws IllegalArgumentException { - double A, B, x[], F[] = new double[f.length]; + final double transformed[] = new double[f.length]; FastFourierTransformer.verifyDataSet(f); if (f[0] != 0.0) { @@ -220,33 +220,33 @@ "first element is not 0: {0}", f[0]); } - int N = f.length; - if (N == 1) { // trivial case - F[0] = 0.0; - return F; + final int n = f.length; + if (n == 1) { // trivial case + transformed[0] = 0.0; + return transformed; } // construct a new array and perform FFT on it - x = new double[N]; + final double[] x = new double[n]; x[0] = 0.0; - x[N >> 1] = 2.0 * f[N >> 1]; - for (int i = 1; i < (N >> 1); i++) { - A = Math.sin(i * Math.PI / N) * (f[i] + f[N-i]); - B = 0.5 * (f[i] - f[N-i]); - x[i] = A + B; - x[N-i] = A - B; + x[n >> 1] = 2.0 * f[n >> 1]; + for (int i = 1; i < (n >> 1); i++) { + final double a = Math.sin(i * Math.PI / n) * (f[i] + f[n-i]); + final double b = 0.5 * (f[i] - f[n-i]); + x[i] = a + b; + x[n - i] = a - b; } FastFourierTransformer transformer = new FastFourierTransformer(); Complex y[] = transformer.transform(x); // reconstruct the FST result for the original array - F[0] = 0.0; - F[1] = 0.5 * y[0].getReal(); - for (int i = 1; i < (N >> 1); i++) { - F[2*i] = -y[i].getImaginary(); - F[2*i+1] = y[i].getReal() + F[2*i-1]; + transformed[0] = 0.0; + transformed[1] = 0.5 * y[0].getReal(); + for (int i = 1; i < (n >> 1); i++) { + transformed[2 * i] = -y[i].getImaginary(); + transformed[2 * i + 1] = y[i].getReal() + transformed[2 * i - 1]; } - return F; + return transformed; } }