This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
new d8cb52a4ca Resolve more cases about FMA identified by "JDK9" comment.
d8cb52a4ca is described below
commit d8cb52a4ca653c9f32b9ec83181b81e05f5c9c5f
Author: Martin Desruisseaux <[email protected]>
AuthorDate: Thu Dec 29 18:21:41 2022 +0100
Resolve more cases about FMA identified by "JDK9" comment.
https://issues.apache.org/jira/browse/SIS-136
---
.../operation/projection/MeridianArcBased.java | 58 +++++++++++++++-------
.../operation/projection/package-info.java | 2 +-
.../operation/transform/LinearInterpolator1D.java | 2 +-
.../operation/matrix/MatrixTestCase.java | 20 ++++----
.../operation/projection/SinusoidalTest.java | 4 +-
.../org/apache/sis/internal/util/DoubleDouble.java | 56 ++++-----------------
.../org/apache/sis/internal/util/Numerics.java | 8 +++
.../java/org/apache/sis/math/MathFunctions.java | 37 ++++++++------
.../java/org/apache/sis/math/package-info.java | 2 +-
.../org/apache/sis/measure/LinearConverter.java | 9 ++--
.../apache/sis/measure/LinearConverterTest.java | 1 -
11 files changed, 98 insertions(+), 101 deletions(-)
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
index 993181b3e7..b9c784e698 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
@@ -16,6 +16,7 @@
*/
package org.apache.sis.referencing.operation.projection;
+import org.apache.sis.internal.util.Numerics;
import org.apache.sis.internal.util.DoubleDouble;
import static java.lang.Math.*;
@@ -27,7 +28,7 @@ import static java.lang.Math.*;
*
* @author Martin Desruisseaux (MPO, IRD, Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.0
+ * @version 1.4
*
* @see <a href="https://en.wikipedia.org/wiki/Meridian_arc";>Meridian arc on
Wikipedia</a>
*
@@ -113,12 +114,12 @@ abstract class MeridianArcBased extends
NormalizedProjection {
* of 2, which result in exact representations in IEEE 754 double
type. Derivation from Kawase formula can be
* viewed at:
https://svn.apache.org/repos/asf/sis/analysis/Map%20projection%20formulas.ods
*/
- cf0 = -441./0x10000 * e10 + -175./0x4000 * e8 + -5./0x100 *
e6 + -3./0x40 * e4 + -1./4 * e2 + 1;
- cf1 = 1575./0x20000 * e10 + 175./0x4000 * e8 + 5./0x100 *
e6 + 3./0x40 * e4 + -3./4 * e2;
- cf2 = -2625./0x8000 * e10 + 175./0x6000 * e8 + 5120./0x60000 *
e6 + -15./0x20 * e4;
- cf3 = 735./0x800 * e10 + 2240./0x60000 * e8 + -35./0x60 *
e6;
- cf4 = -2205./0x1000 * e10 + -315./0x400 * e8;
- // cf5 = 693./0x20000 * e10 omitted for now (not yet used).
+ cf0 = fma(e10, -441./0x10000, fma(e8, -175./0x4000, fma(e6,
-5./0x100, fma(e4, -3./0x40, fma(e2, -1./4, 1)))));
+ cf1 = fma(e10, 1575./0x20000, fma(e8, 175./0x4000, fma(e6,
5./0x100, fma(e4, 3./0x40, e2*(-3./4)))));
+ cf2 = fma(e10, -2625./0x8000, fma(e8, 175./0x6000, fma(e6,
5120./0x60000, e4*(-15./0x20))));
+ cf3 = fma(e10, 735./0x800, fma(e8, 2240./0x60000, e6*
(-35./0x60)));
+ cf4 = fma(e10, -2205./0x1000, e8*(-315./0x400));
+ // cf5 = e10* (693./0x20000) omitted for now (not yet used).
/*
* Coefficients for inverse transform derived from Snyder 3-26 and
EPSG guidance notes:
*
@@ -134,17 +135,19 @@ abstract class MeridianArcBased extends
NormalizedProjection {
* µ is the rectifying latitude.
* Derivation of coefficients used by this class are provided in the
above-cited spreadsheet.
*/
- rµ = (1 - 1./4*e2 - 3./64*e4 - 5./256*e6); // Part of Snyder
7-19 for computing rectifying latitude.
+ rµ = fma(e6, -5./256,
+ fma(e4, -3./64,
+ fma(e2, -1./4, 1))); // Part of Snyder
7-19 for computing rectifying latitude.
DoubleDouble e1 = initializer.axisLengthRatio();
e1.ratio_1m_1p();
final double ei = e1.doubleValue(); // Equivalent to [1
- √(1 - ℯ²)] / [1 + √(1 - ℯ²)] (Snyder 3-24).
final double ei2 = ei*ei;
final double ei3 = ei*ei2;
final double ei4 = ei2*ei2;
- ci1 = 657./0x40 * ei4 + 31./4 * ei3 + 21./4 * ei2 + 3./1 *
ei;
- ci2 = -5045./0x20 * ei4 + -151./3 * ei3 + -21./2 * ei2;
- ci3 = 3291./0x08 * ei4 + 151./3 * ei3;
- ci4 = -1097./0x04 * ei4;
+ ci1 = fma(ei4, 657./0x40, fma(ei3, 31./4, fma(ei2, 21./4,
ei*(3./1))));
+ ci2 = fma(ei4, -5045./0x20, fma(ei3, -151./3, ei2*(-21./2)));
+ ci3 = fma(ei4, 3291./0x08, ei3* (151./3));
+ ci4 = ei4* (-1097./0x04);
}
/**
@@ -183,7 +186,11 @@ abstract class MeridianArcBased extends
NormalizedProjection {
*/
final double distance(final double φ, final double sinφ, final double
cosφ) {
final double sinφ2 = sinφ * sinφ;
- return cf0*φ + cosφ*sinφ*(cf1 + sinφ2*(cf2 + sinφ2*(cf3 +
sinφ2*cf4))); // TODO: use Math.fma with JDK9.
+ if (Numerics.USE_FMA) {
+ return fma(cosφ*sinφ, fma(sinφ2, fma(sinφ2, fma(sinφ2, cf4, cf3),
cf2), cf1), cf0*φ);
+ } else {
+ return cf0*φ + cosφ*sinφ*(cf1 + sinφ2*(cf2 + sinφ2*(cf3 +
sinφ2*cf4)));
+ }
}
/**
@@ -192,11 +199,19 @@ abstract class MeridianArcBased extends
NormalizedProjection {
* @return the derivative at the specified latitude.
*/
final double dM_dφ(final double sinφ2) {
- return ((((7 - 8*sinφ2)*cf4 - 6*cf3) * sinφ2
- + 5*cf3 - 4*cf2) * sinφ2
- + 3*cf2 - 2*cf1) * sinφ2
- + cf1
- + cf0;
+ if (Numerics.USE_FMA) {
+ return fma(fma(fma(fma(fma(-8, sinφ2, 7),
+ cf4, -6*cf3), sinφ2,
+ 5*cf3 - 4*cf2), sinφ2,
+ 3*cf2 - 2*cf1), sinφ2,
+ cf1) + cf0;
+ } else {
+ return ((((7 + -8*sinφ2)*cf4 - 6*cf3) * sinφ2
+ + 5*cf3 - 4*cf2) * sinφ2
+ + 3*cf2 - 2*cf1) * sinφ2
+ + cf1
+ + cf0;
+ }
}
/**
@@ -207,9 +222,14 @@ abstract class MeridianArcBased extends
NormalizedProjection {
*/
final double latitude(final double distance) {
double φ = distance / rµ; // Rectifying latitude µ
used as first approximation.
+ double cosφ = cos(φ);
double sinφ = sin(φ);
double sinφ2 = sinφ * sinφ;
- φ += cos(φ)*sinφ*(ci1 + sinφ2*(ci2 + sinφ2*(ci3 + sinφ2*ci4)));
// Snyder 3-26.
+ if (Numerics.USE_FMA) {
+ φ = fma(cosφ*sinφ, fma(sinφ2, fma(sinφ2, fma(sinφ2, ci4, ci3),
ci2), ci1), φ);
+ } else {
+ φ += cosφ*sinφ*(ci1 + sinφ2*(ci2 + sinφ2*(ci3 + sinφ2*ci4)));
// Snyder 3-26.
+ }
/*
* We could improve accuracy by continuing from here with Newton's
iterative method
* (see MeridianArcTest.inverse(…) for implementation). However, those
iterations requires
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/package-info.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/package-info.java
index 28a6bcbb26..0fc98bd16b 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/package-info.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/package-info.java
@@ -160,7 +160,7 @@
* @author Rémi Maréchal (Geomatys)
* @author Adrian Custer (Geomatys)
* @author Matthieu Bastianelli (Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see <a href="https://mathworld.wolfram.com/MapProjection.html";>Map
projections on MathWorld</a>
*
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearInterpolator1D.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearInterpolator1D.java
index 642f55ee41..b18873bf7d 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearInterpolator1D.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearInterpolator1D.java
@@ -123,7 +123,7 @@ final class LinearInterpolator1D extends
AbstractMathTransform1D implements Seri
return LinearTransform1D.create(slope, offset);
}
value = values[i];
- } while (Numerics.epsilonEqual(value, offset + slope*i, // TODO:
use Math.fma with JDK9.
+ } while (Numerics.epsilonEqual(value, Math.fma(i, slope, offset),
Math.max(Math.abs(value), as) *
Numerics.COMPARISON_THRESHOLD));
/*
* If the values are in decreasing order, reverse their sign so we get
increasing order.
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java
index 9637f91165..c46f9dd025 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatrixTestCase.java
@@ -45,7 +45,7 @@ import static org.apache.sis.test.Assert.*;
* <p>This class uses <a
href="http://math.nist.gov/javanumerics/jama";>JAMA</a> as the reference
implementation.</p>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.4
* @since 0.4
*/
public abstract class MatrixTestCase extends TestCase {
@@ -184,7 +184,7 @@ public abstract class MatrixTestCase extends TestCase {
{
assertEquals("numRow", numRow, actual.getNumRow());
assertEquals("numCol", numCol, actual.getNumCol());
- assertArrayEquals(expected, actual.getElements(), tolerance); // First
because more informative in case of failure.
+ assertArrayEquals(expected, actual.getElements(), tolerance);
// First because more informative in case of failure.
assertTrue(Matrices.create(numRow, numCol, expected).equals(actual,
tolerance));
}
@@ -525,8 +525,8 @@ public abstract class MatrixTestCase extends TestCase {
if ((n % 10) == 0) {
setRandomValues(at, matrix);
}
- at.translate(vector[0] = random.nextDouble() * 50 - 25,
- vector[1] = random.nextDouble() * 50 - 25);
+ at.translate(vector[0] = Math.fma(random.nextDouble(), 50, -25),
+ vector[1] = Math.fma(random.nextDouble(), 50, -25));
matrix.translate(vector);
assertMatrixEquals("translate", AffineTransforms2D.toMatrix(at),
matrix, TOLERANCE);
}
@@ -538,8 +538,8 @@ public abstract class MatrixTestCase extends TestCase {
private void setRandomValues(final AffineTransform at, final MatrixSIS
matrix) {
at.setToRotation(random.nextDouble() * StrictMath.PI);
at.scale(nextNonZeroRandom(), nextNonZeroRandom());
- at.translate(random.nextDouble() * 100 - 50,
- random.nextDouble() * 100 - 50);
+ at.translate(Math.fma(random.nextDouble(), 100, - 50),
+ Math.fma(random.nextDouble(), 100, - 50));
matrix.setElements(new double[] {
at.getScaleX(), at.getShearX(), at.getTranslateX(),
at.getShearY(), at.getScaleY(), at.getTranslateY(),
@@ -564,8 +564,8 @@ public abstract class MatrixTestCase extends TestCase {
if ((n % 10) == 0) {
setRandomValues(at, matrix);
}
- vector[0] = random.nextDouble() * 50 - 25;
- vector[1] = random.nextDouble() * 50 - 25;
+ vector[0] = Math.fma(random.nextDouble(), 50, -25);
+ vector[1] = Math.fma(random.nextDouble(), 50, -25);
vector[2] = 1;
final double[] result = matrix.multiply(vector); // The
result to verify.
at.transform(vector, 0, vector, 0, 1); // The
expected result.
@@ -596,7 +596,7 @@ public abstract class MatrixTestCase extends TestCase {
final int nx = random.nextInt(8) + 1;
elements = new double[numCol * nx];
for (int k=0; k<elements.length; k++) {
- elements[k] = 8 - random.nextDouble() * 10;
+ elements[k] = Math.fma(random.nextDouble(), -10, 8);
}
final Matrix referenceArg = new Matrix(elements, nx).transpose();
final MatrixSIS matrixArg = Matrices.create(numCol, nx, elements);
@@ -636,7 +636,7 @@ public abstract class MatrixTestCase extends TestCase {
final int nx = random.nextInt(8) + 1;
elements = new double[numCol * nx];
for (int k=0; k<elements.length; k++) {
- elements[k] = 8 - random.nextDouble() * 10;
+ elements[k] = Math.fma(random.nextDouble(), -10, 8);
}
final Matrix referenceArg = new Matrix(elements, nx).transpose();
final MatrixSIS matrixArg = Matrices.create(numCol, nx, elements);
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
index 4a283b6fc4..9d487a3617 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
@@ -28,7 +28,7 @@ import org.junit.Test;
* Tests the {@link Sinusoidal} projection.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.4
* @since 1.0
*/
@DependsOn(MeridianArcTest.class)
@@ -123,7 +123,7 @@ public final class SinusoidalTest extends
MapProjectionTestCase {
createProjection(true);
final double delta = (1.0 / 60) / 1852; //
Approximately 1 metre.
derivativeDeltas = new double[] {delta, delta};
- tolerance = Formulas.LINEAR_TOLERANCE / 10000;
+ tolerance = Formulas.LINEAR_TOLERANCE / 10;
verifyDerivative(105, 30);
verifyDerivative(100, -60);
}
diff --git
a/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
b/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
index ad06dc676b..af78f73076 100644
---
a/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
+++
b/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
@@ -46,17 +46,8 @@ import org.apache.sis.math.DecimalFunctions;
* BigDecimal decimal = new BigDecimal(dd.value).add(new
BigDecimal(dd.error));
* }
*
- * <h2>Impact of availability of FMA instructions</h2>
- * When allowed to use <cite>fused multiply-add</cite> (FMA) instruction added
in JDK9
- * (see <a href="https://issues.apache.org/jira/browse/SIS-136";>SIS-136</a> on
Apache SIS JIRA),
- * then the following methods should be revisited:
- *
- * <ul>
- * <li>{@link #setToProduct(double, double)} - revisit with [Hida & al.]
algorithm 7.</li>
- * </ul>
- *
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see <a
href="https://en.wikipedia.org/wiki/Double-double_%28arithmetic%29#Double-double_arithmetic";>Wikipedia:
Double-double arithmetic</a>
*
@@ -92,29 +83,6 @@ public final class DoubleDouble extends Number {
*/
private static final int ZERO_THRESHOLD = 2;
- /**
- * The split constant used as part of multiplication algorithms. The split
algorithm is as below
- * (we have to inline it in multiplication methods because Java cannot
return multi-values):
- *
- * {@snippet lang="java" :
- * private void split(double a) {
- * double t = SPLIT * a;
- * double ahi = t - (t - a);
- * double alo = a - ahi;
- * }
- * }
- *
- * <p>Source: [Hida & al.] page 4 algorithm 5, itself reproduced from
[Shewchuk] page 325.</p>
- */
- private static final double SPLIT = (1 << 27) + 1;
-
- /**
- * Maximal value that can be handled by {@link #multiply(double, double)}.
- * If a multiplication is using a value greater than {@code MAX_VALUE},
- * then the result will be infinity or NaN.
- */
- public static final double MAX_VALUE = Double.MAX_VALUE / SPLIT;
-
/**
* Pre-defined constants frequently used in SIS, sorted in increasing
order. This table contains only
* constants that cannot be inferred by {@link
DecimalFunctions#deltaForDoubleToDecimal(double)},
@@ -484,20 +452,16 @@ public final class DoubleDouble extends Number {
* Sets this {@code DoubleDouble} to the product of the given numbers.
* The given numbers shall not be greater than {@value #MAX_VALUE} in
magnitude.
*
- * <p>Source: [Hida & al.] page 4 algorithm 6, itself reproduced from
[Shewchuk] page 326.</p>
+ * <p>Source: [Hida & al.] page 5 algorithm 7, itself reproduced from
[Shewchuk] page 326.
+ * This is the algorithm used when FMA instruction is available. For an
algorithm without FMA,
+ * see [Hida & al.] page 4 algorithm 6 (it is more complicated).</p>
*
* @param a the first number to multiply.
* @param b the second number to multiply.
*/
public void setToProduct(final double a, final double b) {
value = a * b;
- double t = SPLIT * a;
- final double ahi = t - (t - a);
- final double alo = a - ahi;
- t = SPLIT * b;
- final double bhi = t - (t - b);
- final double blo = b - bhi;
- error = ((ahi*bhi - value) + ahi*blo + alo*bhi) + alo*blo;
+ error = Math.fma(a, b, -value); // Really needs the `fma(…)`
precision here.
if (DISABLED) error = 0;
}
@@ -656,9 +620,9 @@ public final class DoubleDouble extends Number {
if (value == 0 && error != 0) {
/*
* The two values almost cancelled, only their error terms are
different.
- * The number of significand bits (mantissa) in the IEEE 'double'
representation is 52,
+ * The number of significand bits (mantissa) in the IEEE `double`
representation is 52,
* not counting the hidden bit. So estimate the accuracy of the
double-double number as
- * the accuracy of the 'double' value (which is 1 ULP) scaled as
if we had 52 additional
+ * the accuracy of the `double` value (which is 1 ULP) scaled as
if we had 52 additional
* significand bits (we ignore some more bits if ZERO_THRESHOLD is
greater than 0).
* If the error is not greater than that value, then assume that
it is not significant.
*/
@@ -1073,19 +1037,19 @@ public final class DoubleDouble extends Number {
}
final double denominatorValue = value;
/*
- * The 'b * (a.value / b.value)' part in the method javadoc.
+ * The `b * (a.value / b.value)` part in the method javadoc.
*/
final double quotient = numeratorValue / denominatorValue;
multiply(quotient);
/*
- * Compute 'remainder' as 'a - above_product'.
+ * Compute `remainder` as `a - above_product`.
*/
final double productError = error;
setToSum(numeratorValue, -value);
error -= productError; // Complete the above
subtraction
error += numeratorError;
/*
- * Adds the 'remainder / b' term, using 'remainder / b.value' as an
approximation
+ * Adds the `remainder / b` term, using `remainder / b.value` as an
approximation
* (otherwise we would have to invoke this method recursively). The
approximation
* is assumed okay since the second term is small compared to the
first one.
*/
diff --git
a/core/sis-utility/src/main/java/org/apache/sis/internal/util/Numerics.java
b/core/sis-utility/src/main/java/org/apache/sis/internal/util/Numerics.java
index 7c2e038a4c..f4f2c9d68a 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/internal/util/Numerics.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/internal/util/Numerics.java
@@ -74,6 +74,14 @@ public final class Numerics extends Static {
boxed = -value; CACHE.put(boxed, boxed);
}
+ /**
+ * Whether to use {@link Math#fma(double, double, double)} for performance
reasons.
+ * We do not use this flag when the goal is to get better accuracy rather
than performance.
+ * Use of FMA brings performance benefits on machines having hardware
support,
+ * but come at a high cost on older machines without hardware support.
+ */
+ public static final boolean USE_FMA = true;
+
/**
* Maximum number of rows or columns in Apache SIS matrices. We define a
maximum because SIS is expected to work
* mostly with small matrices, because their sizes are related to the
number of dimensions in coordinate systems.
diff --git
a/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
b/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
index 90b830745a..f3da6bd283 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
@@ -23,9 +23,16 @@ import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.internal.util.DoubleDouble;
+import static java.lang.Math.PI;
+import static java.lang.Math.min;
import static java.lang.Math.abs;
import static java.lang.Math.sqrt;
+import static java.lang.Math.cbrt;
+import static java.lang.Math.fma;
+import static java.lang.Math.cos;
+import static java.lang.Math.copySign;
import static java.lang.Math.multiplyFull;
+import static java.lang.Math.multiplyExact;
import static java.lang.Float.intBitsToFloat;
import static java.lang.Float.floatToIntBits;
import static java.lang.Float.floatToRawIntBits;
@@ -62,7 +69,7 @@ import static
org.apache.sis.internal.util.Numerics.SIGNIFICAND_SIZE;
*
* @author Martin Desruisseaux (MPO, IRD, Geomatys)
* @author Johann Sorel (Geomatys)
- * @version 1.2
+ * @version 1.4
*
* @see DecimalFunctions
* @see org.apache.sis.util.Numbers
@@ -360,9 +367,9 @@ public final class MathFunctions extends Static {
result = base;
}
while ((exponent >>>= 1) != 0) {
- base = Math.multiplyExact(base, base);
+ base = multiplyExact(base, base);
if ((exponent & 1) != 0) {
- result = Math.multiplyExact(result, base);
+ result = multiplyExact(result, base);
}
}
} else if (exponent < 0) {
@@ -801,7 +808,7 @@ public final class MathFunctions extends Static {
int i = primes.length;
int n = Short.toUnsignedInt(primes[i - 1]);
// Compute by block of 16 values, for reducing the amount
of array resize.
- primes = Arrays.copyOf(primes, Math.min((index | 0xF) + 1,
PRIMES_LENGTH_16_BITS));
+ primes = Arrays.copyOf(primes, min((index | 0xF) + 1,
PRIMES_LENGTH_16_BITS));
do {
testNextNumber: while (true) { // Simulate a "goto" statement
(usually not recommanded...)
final int stopAt = (int) sqrt(n += 2);
@@ -839,7 +846,7 @@ testNextNumber: while (true) { // Simulate a
"goto" statement (usua
ArgumentChecks.ensureBetween("number", 2,
HIGHEST_SUPPORTED_PRIME_NUMBER, number);
final short[] primes = MathFunctions.primes;
int lower = 0;
- int upper = Math.min(PRIMES_LENGTH_15_BITS, primes.length);
+ int upper = min(PRIMES_LENGTH_15_BITS, primes.length);
if (number > Short.MAX_VALUE) {
lower = upper;
upper = primes.length;
@@ -1038,7 +1045,7 @@ testNextNumber: while (true) { // Simulate a
"goto" statement (usua
*/
case 2: {
final double b = coefficients[lower + 1];
- final double q = -0.5 * (b + Math.copySign(sqrt(b*b -
4*a*c), b));
+ final double q = -0.5 * (b + copySign(sqrt(b*b - 4*a*c),
b));
final double x1 = q/a;
final double x2 = c/q;
if (Double.isNaN(x1) && Double.isNaN(x2)) break;
@@ -1069,9 +1076,9 @@ testNextNumber: while (true) { // Simulate a
"goto" statement (usua
b = coefficients[lower + 2] / a;
a = coefficients[lower + 3] / a;
final double a2 = a*a;
- final double p = -3./8 * (a2) + b;
- final double q = 1./8 * (a2*a) - 1./2 * (a*b) +
c;
- final double r = -3./256 * (a2*a2) + 1./16 * (a2*b) -
1./4 * (a*c) + d;
+ final double p = fma(-3./8, a2, b);
+ final double q = fma( 1./8* a2 - 1./2 *b, a, c);
+ final double r = fma(-3./256, a2, 1./16*b)*a2 +
fma(-1./4*a, c, d);
final double[] roots = solveCubic(-2*p, p*p-4*r, q*q,
true);
if (roots.length != 4) break;
for (int i=0; i<3; i++) {
@@ -1116,7 +1123,7 @@ testNextNumber: while (true) { // Simulate a
"goto" statement (usua
final double R = (a*(a*a - 4.5*b) + 13.5*c) / 27; // R from
Numerical Recipes 5.6.10.
final double Q3 = Q*Q*Q;
final double R2 = R*R;
- a /= 3; // Last
term of Numerical Recipes 5.6.12, 17 and 18.
+ a /= -3; // Last
term of Numerical Recipes 5.6.12, 17 and 18.
if (R2 < Q3) {
/*
* Numerical Recipes 5.6.11 and 5.6.12 uses acos(R/sqrt(Q³)). It
is possible to rewrite as
@@ -1126,17 +1133,17 @@ testNextNumber: while (true) { // Simulate
a "goto" statement (usua
b = Math.acos(R/sqrt(Q3)) / 3; // θ from
Numerical recipes 5.6.11, then b = θ/3.
c = -2 * sqrt(Q); // First
part of Numerical Recipes 5.6.12.
double[] roots = new double[quartic ? 4 : 3];
- roots[2] = c*Math.cos(b - 2*Math.PI/3) - a; // TODO:
try Math.fma with JDK9.
- roots[1] = c*Math.cos(b + 2*Math.PI/3) - a;
- roots[0] = c*Math.cos(b) - a;
+ roots[2] = fma(c, cos(b - 2*PI/3), a);
+ roots[1] = fma(c, cos(b + 2*PI/3), a);
+ roots[0] = fma(c, cos(b), a);
if (!quartic) {
roots = removeDuplicated(roots);
}
return roots;
}
if (!quartic) {
- b = -Math.copySign(Math.cbrt(abs(R) + sqrt(R2 - Q3)), R); //
A from Numerical Recipes 5.6.15.
- final double x = (b == 0 ? 0 : b + Q/b) - a;
+ b = -copySign(cbrt(abs(R) + sqrt(R2 - Q3)), R); // A from
Numerical Recipes 5.6.15.
+ final double x = (b == 0 ? 0 : b + Q/b) + a;
if (!Double.isNaN(x)) {
return new double[] {x};
}
diff --git
a/core/sis-utility/src/main/java/org/apache/sis/math/package-info.java
b/core/sis-utility/src/main/java/org/apache/sis/math/package-info.java
index e0aef6d448..4d8485cb46 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/math/package-info.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/math/package-info.java
@@ -40,7 +40,7 @@
* </ul>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.2
+ * @version 1.4
* @since 0.3
*/
package org.apache.sis.math;
diff --git
a/core/sis-utility/src/main/java/org/apache/sis/measure/LinearConverter.java
b/core/sis-utility/src/main/java/org/apache/sis/measure/LinearConverter.java
index 11909ba657..d75cc6884a 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/measure/LinearConverter.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/measure/LinearConverter.java
@@ -33,16 +33,16 @@ import org.apache.sis.internal.util.Numerics;
* Note that the "linear" word in this class does not have the same meaning
than the same word
* in the {@link #isLinear()} method inherited from JSR-385.
*
- * <p><b>Implementation note:</b>
+ * <h2>Implementation note</h2>
* for performance reason we should create specialized subtypes for the case
where there is only a scale to apply,
* or only an offset, <i>etc.</i> But we don't do that in Apache SIS
implementation because we will rarely use the
* {@code UnitConverter} for converting a lot of values. We rather use {@code
MathTransform} for operations on
* <var>n</var>-dimensional tuples, and unit conversions are only a special
case of those more generic operations.
* The {@code sis-referencing} module provided the specialized implementations
needed for efficient operations
- * and know how to copy the {@code UnitConverter} coefficients into an affine
transform matrix.</p>
+ * and know how to copy the {@code UnitConverter} coefficients into an affine
transform matrix.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.4
* @since 0.8
*/
final class LinearConverter extends AbstractConverter implements
LenientComparable {
@@ -273,8 +273,7 @@ final class LinearConverter extends AbstractConverter
implements LenientComparab
*/
@Override
public double convert(final double value) {
- // TODO: use JDK9' Math.fma(…) and verify if it solve the accuracy
issue in LinearConverterTest.inverse().
- return (value * scale + offset) / divisor;
+ return Math.fma(value, scale, offset) / divisor;
}
/**
diff --git
a/core/sis-utility/src/test/java/org/apache/sis/measure/LinearConverterTest.java
b/core/sis-utility/src/test/java/org/apache/sis/measure/LinearConverterTest.java
index b0277ac213..5c8fce4e25 100644
---
a/core/sis-utility/src/test/java/org/apache/sis/measure/LinearConverterTest.java
+++
b/core/sis-utility/src/test/java/org/apache/sis/measure/LinearConverterTest.java
@@ -176,7 +176,6 @@ public final class LinearConverterTest extends TestCase {
c = LinearConverter.offset(27315, 100);
inv = (LinearConverter) c.inverse();
assertEquals(12.3, c.convert(inv.convert(12.3)), 1E-13);
- // TODO: use JDK9' Math.fma(…) in LinearConverter.convert(double) and
verify if it solve the accuracy issue.
}
/**
