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
commit e9c3ce60ea17cd871e3a3390f435cb99d3555b49 Author: Martin Desruisseaux <[email protected]> AuthorDate: Tue Jan 6 11:23:59 2026 +0100 Remove the `USE_FMA` internal flag. Was always on for last years. --- .../sis/referencing/internal/shared/Formulas.java | 9 ---- .../operation/projection/AuthalicConversion.java | 14 +++---- .../operation/projection/ConformalProjection.java | 8 +--- .../operation/projection/MeridianArcBased.java | 33 ++++----------- .../transform/AbstractLinearTransform.java | 7 +--- .../operation/transform/LinearTransform1D.java | 49 ++++------------------ .../operation/transform/ProjectiveTransform.java | 25 ++--------- .../transform/ProjectiveTransformTest.java | 13 ++---- 8 files changed, 32 insertions(+), 126 deletions(-) diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/shared/Formulas.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/shared/Formulas.java index 6b902b6a73..bf49dd3e57 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/shared/Formulas.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/shared/Formulas.java @@ -91,15 +91,6 @@ public final class Formulas { @Configuration public static final int MAXIMUM_ITERATIONS = 18; - /** - * 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. - */ - @Configuration - public static final boolean USE_FMA = true; - /** * Do not allow instantiation of this class. */ diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/AuthalicConversion.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/AuthalicConversion.java index 9623175770..a6c8df6df5 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/AuthalicConversion.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/AuthalicConversion.java @@ -18,7 +18,6 @@ package org.apache.sis.referencing.operation.projection; import static java.lang.Math.*; import org.apache.sis.referencing.internal.Resources; -import org.apache.sis.referencing.internal.shared.Formulas; import static org.apache.sis.math.MathFunctions.atanh; @@ -247,15 +246,12 @@ abstract class AuthalicConversion extends NormalizedProjection { final double sinβ2 = sinβ * sinβ; final double β = asin(sinβ); /* - * Snyder 3-18, but rewriten using trigonometric identities in order to avoid - * multiple calls to sin(double) method. + * Snyder 3-18, but rewriten using trigonometric identities + * in order to avoid multiple calls to sin(double) method. + * + * φ = β + cos(β)*sinβ*(c2β + sinβ2*(c4β + sinβ2*c6β)) */ - double φ; - if (Formulas.USE_FMA) { - φ = fma(fma(fma(sinβ2, c6β, c4β), sinβ2, c2β), cos(β)*sinβ, β); - } else { - φ = β + cos(β)*sinβ*(c2β + sinβ2*(c4β + sinβ2*c6β)); - } + double φ = fma(fma(fma(sinβ2, c6β, c4β), sinβ2, c2β), cos(β)*sinβ, β); if (useIterations) { /* * Mathematical note: Snyder 3-16 gives q/(1-ℯ²) instead of y in the calculation of Δφ below. diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/ConformalProjection.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/ConformalProjection.java index eda2f01ba1..ca35ec0854 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/ConformalProjection.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/ConformalProjection.java @@ -18,7 +18,6 @@ package org.apache.sis.referencing.operation.projection; import static java.lang.Math.*; import org.apache.sis.referencing.internal.Resources; -import org.apache.sis.referencing.internal.shared.Formulas; /** @@ -207,11 +206,8 @@ abstract class ConformalProjection extends NormalizedProjection { */ final double sin_2φ = sin(2*φ); final double cos_2φ = cos(2*φ); - if (Formulas.USE_FMA) { - φ = fma(sin_2φ, fma(cos_2φ, fma(cos_2φ, fma(cos_2φ, c8χ, c6χ), c4χ), c2χ), φ); - } else { - φ += sin_2φ * (c2χ + cos_2φ * (c4χ + cos_2φ * (c6χ + cos_2φ * c8χ))); - } + φ = fma(sin_2φ, fma(cos_2φ, fma(cos_2φ, fma(cos_2φ, c8χ, c6χ), c4χ), c2χ), φ); + // Equivalent to: φ += sin_2φ * (c2χ + cos_2φ * (c4χ + cos_2φ * (c6χ + cos_2φ * c8χ))) /* * Note: a previous version checked if the value of the smallest term c8χ⋅sin(8φ) was smaller than * the iteration tolerance. But this was not reliable enough. We use now a hard coded threshold diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/MeridianArcBased.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/MeridianArcBased.java index 5c8021b57e..2e398cb081 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/MeridianArcBased.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/MeridianArcBased.java @@ -17,7 +17,6 @@ package org.apache.sis.referencing.operation.projection; import static java.lang.Math.*; -import org.apache.sis.referencing.internal.shared.Formulas; /** @@ -181,11 +180,8 @@ abstract class MeridianArcBased extends NormalizedProjection { */ final double distance(final double φ, final double sinφ, final double cosφ) { final double sinφ2 = sinφ * sinφ; - if (Formulas.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))); - } + return fma(cosφ*sinφ, fma(sinφ2, fma(sinφ2, fma(sinφ2, cf4, cf3), cf2), cf1), cf0*φ); + // Equivalent to: cf0*φ + cosφ*sinφ*(cf1 + sinφ2*(cf2 + sinφ2*(cf3 + sinφ2*cf4))) } /** @@ -194,19 +190,11 @@ abstract class MeridianArcBased extends NormalizedProjection { * @return the derivative at the specified latitude. */ final double dM_dφ(final double sinφ2) { - if (Formulas.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; - } + 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; } /** @@ -220,11 +208,8 @@ abstract class MeridianArcBased extends NormalizedProjection { double cosφ = cos(φ); double sinφ = sin(φ); double sinφ2 = sinφ * sinφ; - if (Formulas.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. - } + φ = fma(cosφ*sinφ, fma(sinφ2, fma(sinφ2, fma(sinφ2, ci4, ci3), ci2), ci1), φ); + // Equivalent to: φ += 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/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java index 59ef1b3145..03f4c38c44 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java @@ -25,7 +25,6 @@ import org.opengis.referencing.operation.NoninvertibleTransformException; import org.apache.sis.referencing.operation.matrix.Matrices; import org.apache.sis.referencing.operation.provider.Affine; import org.apache.sis.referencing.internal.Resources; -import org.apache.sis.referencing.internal.shared.Formulas; import org.apache.sis.util.ComparisonMode; import org.apache.sis.util.internal.shared.Numerics; import org.apache.sis.util.resources.Errors; @@ -237,11 +236,7 @@ abstract class AbstractLinearTransform extends AbstractMathTransform implements for (int i=0; i<srcDim; i++) { final double e = getElement(j, i); if (e != 0) { // See the comment in ProjectiveTransform for the purpose of this test. - if (Formulas.USE_FMA) { - sum = Math.fma(srcPts[srcOff + i], e, sum); - } else { - sum += srcPts[srcOff + i] * e; - } + sum = Math.fma(srcPts[srcOff + i], e, sum); } } buffer[j] = sum; diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/LinearTransform1D.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/LinearTransform1D.java index 97bc219999..af601d99e5 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/LinearTransform1D.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/LinearTransform1D.java @@ -28,7 +28,6 @@ import org.apache.sis.referencing.operation.matrix.Matrices; import org.apache.sis.referencing.operation.matrix.Matrix1; import org.apache.sis.referencing.internal.Arithmetic; import org.apache.sis.referencing.internal.shared.ExtendedPrecisionMatrix; -import org.apache.sis.referencing.internal.shared.Formulas; import org.apache.sis.referencing.operation.provider.Affine; import org.apache.sis.util.ComparisonMode; import org.apache.sis.util.internal.shared.DoubleDouble; @@ -282,11 +281,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo */ @Override public double transform(final double value) { - if (Formulas.USE_FMA) { - return Math.fma(value, scale, offset); - } else { - return offset + scale * value; - } + return Math.fma(value, scale, offset); } /** @@ -300,11 +295,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo final boolean derivate) { if (dstPts != null) { - if (Formulas.USE_FMA) { - dstPts[dstOff] = Math.fma(srcPts[srcOff], scale, offset); - } else { - dstPts[dstOff] = offset + scale*srcPts[srcOff]; - } + dstPts[dstOff] = Math.fma(srcPts[srcOff], scale, offset); } return derivate ? new Matrix1(scale) : null; } @@ -319,21 +310,13 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo { if (srcPts != dstPts || srcOff >= dstOff) { while (--numPts >= 0) { - if (Formulas.USE_FMA) { - dstPts[dstOff++] = Math.fma(srcPts[srcOff++], scale, offset); - } else { - dstPts[dstOff++] = offset + scale * srcPts[srcOff++]; - } + dstPts[dstOff++] = Math.fma(srcPts[srcOff++], scale, offset); } } else { srcOff += numPts; dstOff += numPts; while (--numPts >= 0) { - if (Formulas.USE_FMA) { - dstPts[--dstOff] = Math.fma(srcPts[--srcOff], scale, offset); - } else { - dstPts[--dstOff] = offset + scale * srcPts[--srcOff]; - } + dstPts[--dstOff] = Math.fma(srcPts[--srcOff], scale, offset); } } } @@ -349,21 +332,13 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo { if (srcPts != dstPts || srcOff >= dstOff) { while (--numPts >= 0) { - if (Formulas.USE_FMA) { - dstPts[dstOff++] = (float) Math.fma(srcPts[srcOff++], scale, offset); - } else { - dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]); - } + dstPts[dstOff++] = (float) Math.fma(srcPts[srcOff++], scale, offset); } } else { srcOff += numPts; dstOff += numPts; while (--numPts >= 0) { - if (Formulas.USE_FMA) { - dstPts[--dstOff] = (float) Math.fma(srcPts[--srcOff], scale, offset); - } else { - dstPts[--dstOff] = (float) (offset + scale * srcPts[--srcOff]); - } + dstPts[--dstOff] = (float) Math.fma(srcPts[--srcOff], scale, offset); } } } @@ -378,11 +353,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo final float [] dstPts, int dstOff, int numPts) { while (--numPts >= 0) { - if (Formulas.USE_FMA) { - dstPts[dstOff++] = (float) Math.fma(srcPts[srcOff++], scale, offset); - } else { - dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]); - } + dstPts[dstOff++] = (float) Math.fma(srcPts[srcOff++], scale, offset); } } @@ -395,11 +366,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo final double[] dstPts, int dstOff, int numPts) { while (--numPts >= 0) { - if (Formulas.USE_FMA) { - dstPts[dstOff++] = Math.fma(srcPts[srcOff++], scale, offset); - } else { - dstPts[dstOff++] = offset + scale * srcPts[srcOff++]; - } + dstPts[dstOff++] = Math.fma(srcPts[srcOff++], scale, offset); } } diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java index 9a0a95b463..ab934631d7 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java @@ -22,7 +22,6 @@ import org.opengis.referencing.operation.Matrix; import org.apache.sis.referencing.internal.Arithmetic; import org.apache.sis.referencing.internal.shared.DirectPositionView; import org.apache.sis.referencing.internal.shared.ExtendedPrecisionMatrix; -import org.apache.sis.referencing.internal.shared.Formulas; import org.apache.sis.referencing.operation.matrix.Matrices; import org.apache.sis.referencing.operation.matrix.MatrixSIS; import org.apache.sis.util.ArgumentChecks; @@ -416,11 +415,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre * getting 2D points from 3D points. In such case, the fact that the excluded dimensions had * NaN values should not force the retained dimensions to get NaN values. */ - if (Formulas.USE_FMA) { - sum = Math.fma(srcPts[srcOff + i], e, sum); - } else { - sum += srcPts[srcOff + i] * e; - } + sum = Math.fma(srcPts[srcOff + i], e, sum); } } buffer[j] = sum; @@ -487,11 +482,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre for (int i=0; i<srcDim; i++) { final double e = elt[mix++]; if (e != 0) { // See comment in transform(double[], ...) - if (Formulas.USE_FMA) { - sum = Math.fma(srcPts[srcOff + i], e, sum); - } else { - sum += srcPts[srcOff + i] * e; - } + sum = Math.fma(srcPts[srcOff + i], e, sum); } } buffer[j] = sum; @@ -530,11 +521,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre for (int i=0; i<srcDim; i++) { final double e = elt[mix++]; if (e != 0) { // See comment in transform(double[], ...) - if (Formulas.USE_FMA) { - sum = Math.fma(srcPts[srcOff + i], e, sum); - } else { - sum += srcPts[srcOff + i] * e; - } + sum = Math.fma(srcPts[srcOff + i], e, sum); } } buffer[j] = sum; @@ -572,11 +559,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre for (int i=0; i<srcDim; i++) { final double e = elt[mix++]; if (e != 0) { // See comment in transform(double[], ...) - if (Formulas.USE_FMA) { - sum = Math.fma(srcPts[srcOff + i], e, sum); - } else { - sum += srcPts[srcOff + i] * e; - } + sum = Math.fma(srcPts[srcOff + i], e, sum); } } buffer[j] = sum; diff --git a/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/operation/transform/ProjectiveTransformTest.java b/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/operation/transform/ProjectiveTransformTest.java index f9443e90d2..b435531bf8 100644 --- a/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/operation/transform/ProjectiveTransformTest.java +++ b/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/operation/transform/ProjectiveTransformTest.java @@ -175,22 +175,15 @@ public class ProjectiveTransformTest extends AffineTransformTest { /** * Tests the concatenation of transforms that would result in rounding errors - * in extended-precision matrix operations were not used. - * - * Actually there are two sources of rounding errors tested by this method. - * The first source is rounding errors caused by matrix multiplications. - * The other source is rounding errors inside the {@code transform(…)} methods, - * which is reduced by a denominator column in {@link ProjectiveTransform#elt}. - * For demonstrating the latter rounding errors, it may be necessary to set the - * {@link org.apache.sis.referencing.internal.shared.Formulas#USE_FMA} flag to {@code false}. + * if extended-precision matrix operations were not used. * * @throws FactoryException if the transform cannot be created. * @throws TransformException if a coordinate conversion failed. */ @Test public void testRoundingErrors() throws FactoryException, TransformException { - final Matrix4 num = new Matrix4(); num.m00 = 2; num.m11 = 3.25; num.m22 = -17; - final Matrix4 den = new Matrix4(); den.m00 = 37; den.m11 = 1000; den.m22 = 127; + final var num = new Matrix4(); num.m00 = 2; num.m11 = 3.25; num.m22 = -17; + final var den = new Matrix4(); den.m00 = 37; den.m11 = 1000; den.m22 = 127; // Add translation terms. num.m03 = 4*37; num.m13 = 17; num.m23 = -2*127;
