Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

tighten checkstyle rules: declaring multiple variables in one statement is now forbidden 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@825919 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-10-16 14:51:55 +00:00
parent 88885e7a1c
commit 0cb01403a1
24 changed files with 493 additions and 336 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
checkstyle.xml
View File

@ -122,15 +122,24 @@
<!-- Don't add up parentheses when they are not required -->
<module name="UnnecessaryParentheses" />
<!-- Don't use too widespread catch (Exception, Throwable, RuntimeException) -->
<!-- Don't use too widespread catch (Exception, Throwable, RuntimeException) -->
<module name="IllegalCatch" />
<!--
<!-- Don't use = or != for string comparisons -->
<module name="StringLiteralEquality" />
<module name="MultipleStringLiterals" />
<!-- Don't declare multiple variables in the same statement -->
<module name="MultipleVariableDeclarations" />
<module name="TodoComment" />
<!-- String literals more than one character long should not be repeated several times -->
<!--
<module name="MultipleStringLiterals" >
<property name="ignoreStringsRegexp" value='^(("")|("."))$'/>
</module>
-->
<!-- <module name="TodoComment" /> -->
</module>
<!-- Require files to end with newline characters -->

16
src/main/java/org/apache/commons/math/analysis/integration/SimpsonIntegrator.java
View File

@ -66,25 +66,22 @@ public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
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 @@ public class SimpsonIntegrator extends UnivariateRealIntegratorImpl {
}
olds = s;
oldt = t;
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}

22
src/main/java/org/apache/commons/math/analysis/integration/TrapezoidIntegrator.java
View File

@ -77,17 +77,15 @@ public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl {
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 @@ public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl {
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 @@ public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl {
}
}
oldt = t;
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}

30
src/main/java/org/apache/commons/math/analysis/interpolation/DividedDifferenceInterpolator.java
View File

@ -57,8 +57,6 @@ public class DividedDifferenceInterpolator implements UnivariateRealInterpolator
* 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 @@ public class DividedDifferenceInterpolator implements UnivariateRealInterpolator
* <p>
* Note x[], y[], a[] have the same length but c[]'s size is one less.</p>
*/
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 @@ public class DividedDifferenceInterpolator implements UnivariateRealInterpolator
* @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);

5
src/main/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java
View File

@ -277,7 +277,10 @@ public class LoessInterpolator
// 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];

9
src/main/java/org/apache/commons/math/analysis/interpolation/MicrosphereInterpolatingFunction.java
View File

@ -162,12 +162,13 @@ public class MicrosphereInterpolatingFunction
// Copy data samples.
samples = new HashMap<RealVector, Double>(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<MicrosphereSurfaceElement>(microsphereElements);

63
src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionLagrangeForm.java
View File

@ -43,9 +43,14 @@ public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
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 class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
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 @@ public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
}
// 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 @@ public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
* @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 @@ public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
}
}
}
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];
}

16
src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionNewtonForm.java
View File

@ -43,11 +43,15 @@ public class PolynomialFunctionNewtonForm implements UnivariateRealFunction {
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 @@ public class PolynomialFunctionNewtonForm implements UnivariateRealFunction {
* 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];

95
src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver.java
View File

@ -139,39 +139,43 @@ public class MullerSolver extends UnivariateRealSolverImpl {
// 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 @@ public class MullerSolver extends UnivariateRealSolverImpl {
y1 = f.value(x1);
oldx = Double.POSITIVE_INFINITY;
}
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
@ -279,38 +282,40 @@ public class MullerSolver extends UnivariateRealSolverImpl {
// 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 @@ public class MullerSolver extends UnivariateRealSolverImpl {
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 @@ public class MullerSolver extends UnivariateRealSolverImpl {
}
// 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);
}

52
src/main/java/org/apache/commons/math/analysis/solvers/RiddersSolver.java
View File

@ -125,34 +125,38 @@ public class RiddersSolver extends UnivariateRealSolverImpl {
// [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 @@ public class RiddersSolver extends UnivariateRealSolverImpl {
// 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;

6
src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java
View File

@ -124,7 +124,8 @@ public abstract class AbstractEstimator implements Estimator {
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 @@ public abstract class AbstractEstimator implements Estimator {
}
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;

47
src/main/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java
View File

@ -256,7 +256,8 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
// 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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
// 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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
// 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 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
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;
}
}
}

110
src/main/java/org/apache/commons/math/linear/BlockFieldMatrix.java
View File

@ -226,11 +226,12 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
// convert array
final Field<T> 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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
}
}
++qBlock;
}
++pBlock;
}
return out;
@ -1273,10 +1294,14 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
public T walkInOptimizedOrder(final FieldMatrixChangingVisitor<T> 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 class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
public T walkInOptimizedOrder(final FieldMatrixPreservingVisitor<T> 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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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 @@ public class BlockFieldMatrix<T extends FieldElement<T>> extends AbstractFieldMa
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;
}
}
}

117
src/main/java/org/apache/commons/math/linear/BlockRealMatrix.java
View File

@ -220,11 +220,12 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
// 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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
sum += tBlock[l++] * mBlock[n];
n += jWidth;
}
outBlock[k++] += sum;
outBlock[k] += sum;
++k;
}
}
}
@ -699,9 +716,11 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
}
}
++qBlock;
}
++pBlock;
}
return out;
@ -1294,10 +1318,14 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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 @@ public class BlockRealMatrix extends AbstractRealMatrix implements Serializable
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;
}
}
}

20
src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
View File

@ -670,10 +670,12 @@ public class EigenDecompositionImpl implements EigenDecomposition {
// 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 @@ public class EigenDecompositionImpl implements EigenDecomposition {
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 @@ public class EigenDecompositionImpl implements EigenDecomposition {
// 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 @@ public class EigenDecompositionImpl implements EigenDecomposition {
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 @@ public class EigenDecompositionImpl implements EigenDecomposition {
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 @@ public class EigenDecompositionImpl implements EigenDecomposition {
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 @@ public class EigenDecompositionImpl implements EigenDecomposition {
work[sixI + 10] = diP1;
work[sixI + 5] = li * t;
pi = pi * t - lambda;
sixI -= 6;
}
work[3] = pi;
work[4] = pi;

7
src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java
View File

@ -417,9 +417,10 @@ public class QRDecompositionImpl implements QRDecomposition {
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];

4
src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
View File

@ -208,10 +208,12 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
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);

19
src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
View File

@ -220,7 +220,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
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 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
// 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 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
// 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 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
// 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;
}
}

50
src/main/java/org/apache/commons/math/transform/FastCosineTransformer.java
View File

@ -218,45 +218,45 @@ public class FastCosineTransformer implements RealTransformer {
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;
}
}

54
src/main/java/org/apache/commons/math/transform/FastFourierTransformer.java
View File

@ -344,58 +344,58 @@ public class FastFourierTransformer implements Serializable {
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 +

36
src/main/java/org/apache/commons/math/transform/FastSineTransformer.java
View File

@ -212,7 +212,7 @@ public class FastSineTransformer implements RealTransformer {
*/
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 @@ public class FastSineTransformer implements RealTransformer {
"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;
}
}

16
src/main/java/org/apache/commons/math/util/MathUtils.java
View File

@ -205,29 +205,35 @@ public final class MathUtils {
long result = 1;
if (n <= 61) {
// For n <= 61, the naive implementation cannot overflow.
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
int i = n - k + 1;
for (int j = 1; j <= k; j++) {
result = result * i / j;
i++;
}
} else if (n <= 66) {
// For n > 61 but n <= 66, the result cannot overflow,
// but we must take care not to overflow intermediate values.
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
int i = n - k + 1;
for (int j = 1; j <= k; j++) {
// We know that (result * i) is divisible by j,
// but (result * i) may overflow, so we split j:
// Filter out the gcd, d, so j/d and i/d are integer.
// result is divisible by (j/d) because (j/d)
// is relative prime to (i/d) and is a divisor of
// result * (i/d).
long d = gcd(i, j);
final long d = gcd(i, j);
result = (result / (j / d)) * (i / d);
i++;
}
} else {
// For n > 66, a result overflow might occur, so we check
// the multiplication, taking care to not overflow
// unnecessary.
for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
long d = gcd(i, j);
int i = n - k + 1;
for (int j = 1; j <= k; j++) {
final long d = gcd(i, j);
result = mulAndCheck(result / (j / d), i / d);
i++;
}
}
return result;

9
src/main/java/org/apache/commons/math/util/OpenIntToDoubleHashMap.java
View File

@ -186,7 +186,8 @@ public class OpenIntToDoubleHashMap implements Serializable {
return missingEntries;
}
for (int perturb = perturb(hash), j = index; states[index] != FREE; perturb >>= PERTURB_SHIFT) {
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {
@ -215,7 +216,8 @@ public class OpenIntToDoubleHashMap implements Serializable {
return false;
}
for (int perturb = perturb(hash), j = index; states[index] != FREE; perturb >>= PERTURB_SHIFT) {
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {
@ -358,7 +360,8 @@ public class OpenIntToDoubleHashMap implements Serializable {
return missingEntries;
}
for (int perturb = perturb(hash), j = index; states[index] != FREE; perturb >>= PERTURB_SHIFT) {
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {

9
src/main/java/org/apache/commons/math/util/OpenIntToFieldHashMap.java
View File

@ -198,7 +198,8 @@ public class OpenIntToFieldHashMap<T extends FieldElement<T>> implements Seriali
return missingEntries;
}
for (int perturb = perturb(hash), j = index; states[index] != FREE; perturb >>= PERTURB_SHIFT) {
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {
@ -227,7 +228,8 @@ public class OpenIntToFieldHashMap<T extends FieldElement<T>> implements Seriali
return false;
}
for (int perturb = perturb(hash), j = index; states[index] != FREE; perturb >>= PERTURB_SHIFT) {
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {
@ -370,7 +372,8 @@ public class OpenIntToFieldHashMap<T extends FieldElement<T>> implements Seriali
return missingEntries;
}
for (int perturb = perturb(hash), j = index; states[index] != FREE; perturb >>= PERTURB_SHIFT) {
int j = index;
for (int perturb = perturb(hash); states[index] != FREE; perturb >>= PERTURB_SHIFT) {
j = probe(perturb, j);
index = j & mask;
if (containsKey(key, index)) {