You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2022/12/30 17:35:31 UTC
[sis] branch geoapi-4.0 updated (d8cb52a4ca -> 17ac6d41c2)
This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a change to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
from d8cb52a4ca Resolve more cases about FMA identified by "JDK9" comment.
new 4f11c0244d Add more verifications of `GeneralMatrix` validity. Fix a bug in `GeneralMatrix.setElements(double[])` identified by above-cited verifications.
new 17ac6d41c2 More attemps to fix accuracy problems, in particular in `LinearTransform1D`. Include deprecation of `InterpolatedMolodenskyTransform`.
The 2 revisions listed above as "new" are entirely new to this
repository and will be described in separate emails. The revisions
listed as "add" were already present in the repository and have only
been added to this reference.
Summary of changes:
.../apache/sis/internal/referencing/Formulas.java | 10 +-
.../sis/internal/referencing/package-info.java | 2 +-
.../provider/MolodenskyInterpolation.java | 3 +
.../operation/CoordinateOperationFinder.java | 19 ++-
.../operation/matrix/GeneralMatrix.java | 67 ++++++--
.../sis/referencing/operation/matrix/Matrices.java | 46 ++---
.../referencing/operation/matrix/MatrixSIS.java | 10 +-
.../referencing/operation/matrix/package-info.java | 2 +-
.../sis/referencing/operation/package-info.java | 2 +-
.../operation/projection/MeridianArcBased.java | 8 +-
.../transform/AbstractLinearTransform.java | 14 +-
.../operation/transform/ConstantTransform1D.java | 13 +-
.../operation/transform/ContextualParameters.java | 18 +-
.../transform/EllipsoidToCentricTransform.java | 1 +
.../transform/ExponentialTransform1D.java | 1 +
.../operation/transform/IdentityTransform1D.java | 4 +-
.../transform/InterpolatedMolodenskyTransform.java | 4 +
.../InterpolatedMolodenskyTransform2D.java | 3 +
.../operation/transform/LinearInterpolator1D.java | 8 +-
.../operation/transform/LinearTransform1D.java | 185 ++++++++++++++++-----
.../operation/transform/MathTransforms.java | 8 +-
.../operation/transform/PoleRotation.java | 1 +
.../operation/transform/ProjectiveTransform.java | 27 ++-
.../operation/CoordinateOperationFinderTest.java | 4 +-
.../operation/matrix/GeneralMatrixTest.java | 5 +-
.../referencing/operation/matrix/MatricesTest.java | 19 ++-
.../operation/matrix/MatrixTestCase.java | 27 +--
.../org/apache/sis/internal/util/DoubleDouble.java | 14 +-
.../org/apache/sis/internal/util/Numerics.java | 8 -
.../org/apache/sis/measure/AbstractConverter.java | 2 +-
.../org/apache/sis/measure/LinearConverter.java | 27 ++-
.../main/java/org/apache/sis/measure/Units.java | 1 +
.../apache/sis/internal/util/DoubleDoubleTest.java | 8 +-
33 files changed, 398 insertions(+), 173 deletions(-)
[sis] 02/02: More attemps to fix accuracy problems, in particular in `LinearTransform1D`. Include deprecation of `InterpolatedMolodenskyTransform`.
Posted by de...@apache.org.
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 17ac6d41c2a3aa68b7c63fb0b8aa5fb42c1f6484
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Fri Dec 30 18:34:44 2022 +0100
More attemps to fix accuracy problems, in particular in `LinearTransform1D`.
Include deprecation of `InterpolatedMolodenskyTransform`.
---
.../apache/sis/internal/referencing/Formulas.java | 10 +-
.../sis/internal/referencing/package-info.java | 2 +-
.../provider/MolodenskyInterpolation.java | 3 +
.../operation/CoordinateOperationFinder.java | 19 ++-
.../operation/matrix/GeneralMatrix.java | 5 +
.../sis/referencing/operation/matrix/Matrices.java | 41 ++---
.../sis/referencing/operation/package-info.java | 2 +-
.../operation/projection/MeridianArcBased.java | 8 +-
.../transform/AbstractLinearTransform.java | 14 +-
.../operation/transform/ConstantTransform1D.java | 13 +-
.../operation/transform/ContextualParameters.java | 18 +-
.../transform/EllipsoidToCentricTransform.java | 1 +
.../transform/ExponentialTransform1D.java | 1 +
.../operation/transform/IdentityTransform1D.java | 4 +-
.../transform/InterpolatedMolodenskyTransform.java | 4 +
.../InterpolatedMolodenskyTransform2D.java | 3 +
.../operation/transform/LinearInterpolator1D.java | 8 +-
.../operation/transform/LinearTransform1D.java | 185 ++++++++++++++++-----
.../operation/transform/MathTransforms.java | 8 +-
.../operation/transform/PoleRotation.java | 1 +
.../operation/transform/ProjectiveTransform.java | 27 ++-
.../operation/CoordinateOperationFinderTest.java | 4 +-
.../referencing/operation/matrix/MatricesTest.java | 19 ++-
.../org/apache/sis/internal/util/Numerics.java | 8 -
.../org/apache/sis/measure/AbstractConverter.java | 2 +-
.../org/apache/sis/measure/LinearConverter.java | 27 ++-
.../main/java/org/apache/sis/measure/Units.java | 1 +
27 files changed, 305 insertions(+), 133 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
index d6c41e07f1..2face4aff3 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
@@ -33,7 +33,7 @@ import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE
* do not want to expose publicly those arbitrary values (or at least not in a too direct way).
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
* @since 0.4
*/
public final class Formulas extends Static {
@@ -98,6 +98,14 @@ public final class Formulas extends Static {
*/
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.
+ */
+ public static final boolean USE_FMA = true;
+
/**
* Do not allow instantiation of this class.
*/
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/package-info.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/package-info.java
index 30f5f420b2..7d28e5d2fd 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/package-info.java
@@ -24,7 +24,7 @@
* may change in incompatible ways in any future version without notice.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
* @since 0.3
*/
package org.apache.sis.internal.referencing;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/MolodenskyInterpolation.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/MolodenskyInterpolation.java
index 3f3a264f57..33f3eb118f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/MolodenskyInterpolation.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/MolodenskyInterpolation.java
@@ -42,8 +42,11 @@ import org.apache.sis.referencing.operation.transform.InterpolatedMolodenskyTran
* @author Martin Desruisseaux (Geomatys)
* @version 0.7
* @since 0.7
+ *
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-500">Deprecate (for removal) InterpolatedMolodenskyTransform</a>
*/
@XmlTransient
+@Deprecated(since="1.4", forRemoval=true)
public final class MolodenskyInterpolation extends FranceGeocentricInterpolation {
/**
* Serial number for inter-operability with different versions.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
index c0e620c9ab..436735e4d3 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
@@ -49,6 +49,7 @@ import org.apache.sis.internal.referencing.provider.GeographicToGeocentric;
import org.apache.sis.internal.referencing.provider.GeocentricToGeographic;
import org.apache.sis.internal.referencing.provider.GeocentricAffine;
import org.apache.sis.internal.referencing.Resources;
+import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.internal.util.Constants;
import org.apache.sis.measure.Units;
import org.apache.sis.metadata.iso.citation.Citations;
@@ -62,6 +63,7 @@ import org.apache.sis.referencing.cs.CoordinateSystems;
import org.apache.sis.referencing.datum.BursaWolfParameters;
import org.apache.sis.referencing.datum.DefaultGeodeticDatum;
import org.apache.sis.referencing.operation.matrix.Matrices;
+import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.resources.Vocabulary;
@@ -113,7 +115,7 @@ import static org.apache.sis.util.Utilities.equalsIgnoreMetadata;
* </ul>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see DefaultCoordinateOperationFactory#createOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, CoordinateOperationContext)
*
@@ -871,9 +873,9 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
* This "epoch shift" is in units of `targetCRS`.
*/
final Unit<Time> targetUnit = targetCS.getAxis(0).getUnit().asType(Time.class);
- double epochShift = sourceDatum.getOrigin().getTime() -
- targetDatum.getOrigin().getTime();
- epochShift = Units.MILLISECOND.getConverterTo(targetUnit).convert(epochShift);
+ DoubleDouble epochShift = new DoubleDouble(sourceDatum.getOrigin().getTime());
+ epochShift.subtract(new DoubleDouble(targetDatum.getOrigin().getTime()));
+ epochShift = DoubleDouble.castOrCopy(Units.MILLISECOND.getConverterTo(targetUnit).convert(epochShift));
/*
* Check axis directions. The method `swapAndScaleAxes` should returns a matrix of size 2×2.
* The element at index (0,0) may be +1 if source and target axes are in the same direction,
@@ -883,15 +885,16 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
*
* The "epoch shift" previously computed is a translation. Consequently, it is added to element (0,1).
*/
- final Matrix matrix;
+ final MatrixSIS matrix;
try {
- matrix = CoordinateSystems.swapAndScaleAxes(sourceCS, targetCS);
+ matrix = MatrixSIS.castOrCopy(CoordinateSystems.swapAndScaleAxes(sourceCS, targetCS));
} catch (IllegalArgumentException | IncommensurableException exception) {
throw new OperationNotFoundException(notFoundMessage(sourceCRS, targetCRS), exception);
}
final int translationColumn = matrix.getNumCol() - 1; // Paranoiac check: should always be 1.
- final double translation = matrix.getElement(0, translationColumn);
- matrix.setElement(0, translationColumn, translation + epochShift);
+ final DoubleDouble translation = DoubleDouble.castOrCopy(matrix.getNumber(0, translationColumn));
+ translation.add(epochShift);
+ matrix.setNumber(0, translationColumn, translation);
return asList(createFromAffineTransform(AXIS_CHANGES, sourceCRS, targetCRS, matrix));
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
index 0fba129a0c..b603a1e18d 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
@@ -90,12 +90,17 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
* If {@code setToIdentity} is {@code true}, then the elements
* on the diagonal (<var>j</var> == <var>i</var>) are set to 1.
*
+ * <p>Do not invoke this constructor directly (except by {@link NonSquareMatrix} constructor) unless
+ * the matrix is known to be square. If this is not the case, invoke a factory method instead.</p>
+ *
* @param numRow number of rows.
* @param numCol number of columns.
* @param setToIdentity {@code true} for initializing the matrix to the identity matrix,
* or {@code false} for leaving it initialized to zero.
* @param precision 1 for normal precision, or 2 for extended precision.
* No other value is allowed (this is not verified).
+ *
+ * @see #createExtendedPrecision(int, int, boolean)
*/
GeneralMatrix(final int numRow, final int numCol, final boolean setToIdentity, final int precision) {
ensureValidSize(numRow, numCol);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
index c1f50e882c..b20281f0aa 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
@@ -253,7 +253,7 @@ public final class Matrices extends Static {
* enough because callers in other package may perform additional arithmetic operations
* on it (for example org.apache.sis.referencing.cs.CoordinateSystems.swapAndScaleAxes).
*/
- final MatrixSIS matrix = new GeneralMatrix(dstAxes.length+1, srcAxes.length+1, false, 2);
+ final MatrixSIS matrix = GeneralMatrix.createExtendedPrecision(dstAxes.length+1, srcAxes.length+1, false);
/*
* Maps source axes to destination axes. If no axis is moved (for example if the user
* want to transform (NORTH,EAST) to (SOUTH,EAST)), then source and destination index
@@ -404,10 +404,9 @@ public final class Matrices extends Static {
* then an exception will be thrown.</li>
* </ul>
*
- * <div class="note"><b>Example:</b>
- * it is legal to transform from (<i>easting</i>, <i>northing</i>, <i>up</i>) to
- * (<i>easting</i>, <i>northing</i>) — this is the first above case — but illegal
- * to transform (<i>easting</i>, <i>northing</i>) to (<i>easting</i>, <i>up</i>).</div>
+ * For example, it is legal to transform from (<i>easting</i>, <i>northing</i>, <i>up</i>)
+ * to (<i>easting</i>, <i>northing</i>) — this is the first above case — but illegal
+ * to transform (<i>easting</i>, <i>northing</i>) to (<i>easting</i>, <i>up</i>).
*
* <h4>Example</h4>
* The following method call:
@@ -444,8 +443,10 @@ public final class Matrices extends Static {
/*
* createTransform(…) may fail if the arrays contain two axes with the same direction, for example
* AxisDirection.OTHER. This check prevents that failure for the common case of an identity transform.
+ * The returned matrix must use extended precision for reason documented in `createTransform(…)`.
*/
- return Matrices.createIdentity(srcAxes.length + 1);
+ final int n = srcAxes.length + 1;
+ return new GeneralMatrix(n, n, true, 2);
}
return createTransform(null, srcAxes, null, dstAxes, false);
}
@@ -647,13 +648,13 @@ public final class Matrices extends Static {
* The 'stride' and 'length' values will be used for computing indices in that array.
* The DoubleDouble temporary object is used only if the array contains error terms.
*/
- final int stride = sourceDimensions;
- final int length = sourceDimensions * targetDimensions;
- final double[] sources = getExtendedElements(subMatrix);
- final DoubleDouble transfer = (sources.length > length) ? new DoubleDouble() : null;
- final MatrixSIS matrix = createZero(targetDimensions-- + expansion,
- sourceDimensions-- + expansion,
- transfer != null);
+ final int stride = sourceDimensions;
+ final int length = sourceDimensions * targetDimensions;
+ final double[] sources = getExtendedElements(subMatrix);
+ final var transfer = (sources.length > length) ? new DoubleDouble() : null;
+ final MatrixSIS matrix = createZero(targetDimensions-- + expansion,
+ sourceDimensions-- + expansion,
+ transfer != null);
/*
* Following code processes from upper row to lower row.
* First, set the diagonal elements on leading new dimensions.
@@ -788,13 +789,13 @@ public final class Matrices extends Static {
if (numRow == srcRow && numCol == srcCol) {
return matrix;
}
- final int stride = srcCol;
- final int length = srcCol * srcRow;
- final double[] sources = getExtendedElements(matrix);
- final DoubleDouble transfer = (sources.length > length) ? new DoubleDouble() : null;
- final MatrixSIS resized = createZero(numRow, numCol, transfer != null);
- final int copyRow = Math.min(--numRow, --srcRow);
- final int copyCol = Math.min(--numCol, --srcCol);
+ final int stride = srcCol;
+ final int length = srcCol * srcRow;
+ final double[] sources = getExtendedElements(matrix);
+ final var transfer = (sources.length > length) ? new DoubleDouble() : null;
+ final MatrixSIS resized = createZero(numRow, numCol, transfer != null);
+ final int copyRow = Math.min(--numRow, --srcRow);
+ final int copyCol = Math.min(--numCol, --srcCol);
for (int j=copyRow; j<numRow; j++) {
resized.setElement(j, j, 1);
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/package-info.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/package-info.java
index 23a4dbf177..6732e798ed 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/package-info.java
@@ -91,7 +91,7 @@
* for example by specifying the area of interest.
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 1.3
+ * @version 1.4
* @since 0.6
*/
@XmlSchema(location = "http://schemas.opengis.net/gml/3.2.1/coordinateOperations.xsd",
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 b9c784e698..4c253979c1 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,8 +16,8 @@
*/
package org.apache.sis.referencing.operation.projection;
-import org.apache.sis.internal.util.Numerics;
import org.apache.sis.internal.util.DoubleDouble;
+import org.apache.sis.internal.referencing.Formulas;
import static java.lang.Math.*;
@@ -186,7 +186,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
*/
final double distance(final double φ, final double sinφ, final double cosφ) {
final double sinφ2 = sinφ * sinφ;
- if (Numerics.USE_FMA) {
+ 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)));
@@ -199,7 +199,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
* @return the derivative at the specified latitude.
*/
final double dM_dφ(final double sinφ2) {
- if (Numerics.USE_FMA) {
+ 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,
@@ -225,7 +225,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
double cosφ = cos(φ);
double sinφ = sin(φ);
double sinφ2 = sinφ * sinφ;
- if (Numerics.USE_FMA) {
+ 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.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
index 3533809e00..e41b38a184 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
@@ -27,6 +27,7 @@ import org.opengis.referencing.operation.NoninvertibleTransformException;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.internal.referencing.provider.Affine;
import org.apache.sis.internal.referencing.Resources;
+import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.internal.util.Numerics;
import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.resources.Errors;
@@ -34,8 +35,8 @@ import org.opengis.util.FactoryException;
/**
- * Base class of linear transforms. For efficiency reasons, this transform implements itself the matrix
- * to be returned by {@link #getMatrix()}.
+ * Base class of linear transforms.
+ * For efficiency reasons, this transform implements itself the matrix to be returned by {@link #getMatrix()}.
*
* <p>Subclasses need to implement the following methods:</p>
* <ul>
@@ -44,7 +45,7 @@ import org.opengis.util.FactoryException;
* </ul>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.4
* @since 0.6
*/
abstract class AbstractLinearTransform extends AbstractMathTransform implements LinearTransform, Matrix, Serializable {
@@ -60,6 +61,7 @@ abstract class AbstractLinearTransform extends AbstractMathTransform implements
*
* @see #inverse()
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
volatile LinearTransform inverse;
/**
@@ -239,7 +241,11 @@ 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.
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ConstantTransform1D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ConstantTransform1D.java
index 47388b935a..9b72b07479 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ConstantTransform1D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ConstantTransform1D.java
@@ -30,7 +30,7 @@ import java.util.Arrays;
* those values to NaN. Overwriting NaN by the constant value is required by {@link org.apache.sis.coverage}.</p>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.5
+ * @version 1.4
* @since 0.5
*/
final class ConstantTransform1D extends LinearTransform1D {
@@ -42,20 +42,21 @@ final class ConstantTransform1D extends LinearTransform1D {
/**
* A transform for the positive zero constant.
*/
- static final ConstantTransform1D ZERO = new ConstantTransform1D(0);
+ static final ConstantTransform1D ZERO = new ConstantTransform1D(0, 0);
/**
* A transform for the one constant.
*/
- static final ConstantTransform1D ONE = new ConstantTransform1D(1);
+ static final ConstantTransform1D ONE = new ConstantTransform1D(1, 0);
/**
* Constructs a new constant transform.
*
- * @param offset the {@code offset} term in the linear equation.
+ * @param offset the {@code offset} term in the linear equation.
+ * @param offsetError error term of {@code offset} for double-double arithmetic.
*/
- ConstantTransform1D(final double offset) {
- super(0, offset);
+ ConstantTransform1D(final double offset, final double offsetError) {
+ super(0, offset, 0, offsetError);
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ContextualParameters.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ContextualParameters.java
index 22891e1f70..9c6c8ecc71 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ContextualParameters.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ContextualParameters.java
@@ -124,7 +124,7 @@ import static java.util.logging.Logger.getLogger;
* Serialization should be used only for short term storage or RMI between applications running the same SIS version.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see org.apache.sis.referencing.operation.projection.NormalizedProjection
* @see AbstractMathTransform#getContextualParameters()
@@ -187,6 +187,7 @@ public class ContextualParameters extends Parameters implements Serializable {
*
* @see #getDescriptor()
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private final ParameterDescriptorGroup descriptor;
/**
@@ -198,6 +199,7 @@ public class ContextualParameters extends Parameters implements Serializable {
*
* @see #getMatrix(MatrixRole)
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private Matrix normalize, denormalize;
/**
@@ -210,6 +212,7 @@ public class ContextualParameters extends Parameters implements Serializable {
* @see #parameter(String)
* @see #freeze()
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private ParameterValue<?>[] values;
/**
@@ -300,19 +303,6 @@ public class ContextualParameters extends Parameters implements Serializable {
isFrozen = true;
}
- /**
- * Creates a matrix for a linear step of the transforms chain.
- * It is important that the matrices created here are instances of {@link MatrixSIS}, in order
- * to allow {@link #getMatrix(MatrixRole)} to return the reference to the (de)normalize matrices.
- */
- private static MatrixSIS linear(final String name, final Integer size) {
- if (size == null) {
- throw new IllegalArgumentException(Errors.format(Errors.Keys.MissingValueForProperty_1, name));
- }
- final int n = size + 1;
- return Matrices.create(n, n, ExtendedPrecisionMatrix.IDENTITY);
- }
-
/**
* Creates a new and frozen {@code ContextualParameters} for the inverse of this operation.
* An optional {@code mapper} can be specified for setting the parameter values.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/EllipsoidToCentricTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/EllipsoidToCentricTransform.java
index af458d6e1f..8ce552ab11 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/EllipsoidToCentricTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/EllipsoidToCentricTransform.java
@@ -214,6 +214,7 @@ public class EllipsoidToCentricTransform extends AbstractMathTransform implement
* construction time). In addition this field is part of serialization form in order to preserve the
* references graph.</div>
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private final AbstractMathTransform inverse;
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ExponentialTransform1D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ExponentialTransform1D.java
index 247acb0f25..03512753b2 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ExponentialTransform1D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ExponentialTransform1D.java
@@ -83,6 +83,7 @@ final class ExponentialTransform1D extends AbstractMathTransform1D implements Se
* The inverse of this transform. Created only when first needed. Serialized in order to avoid
* rounding error if this transform is actually the one which was created from the inverse.
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private MathTransform1D inverse;
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/IdentityTransform1D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/IdentityTransform1D.java
index 44872e4d40..6dbc23b56e 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/IdentityTransform1D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/IdentityTransform1D.java
@@ -22,7 +22,7 @@ package org.apache.sis.referencing.operation.transform;
* This class is a special case of {@link LinearTransform1D} optimized for speed.
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.5
+ * @version 1.4
* @since 0.5
*/
final class IdentityTransform1D extends LinearTransform1D {
@@ -40,7 +40,7 @@ final class IdentityTransform1D extends LinearTransform1D {
* Constructs a new identity transform.
*/
private IdentityTransform1D() {
- super(1, 0);
+ super(1, 0, 0, 0);
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform.java
index 0a308e9876..1008bce1df 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform.java
@@ -59,13 +59,17 @@ import org.apache.sis.util.Debug;
* But in the {@code InterpolatedMolodenskyTransform} case, the interpolated translations are rather the
* ({@linkplain #tX}, {@linkplain #tY}, {@linkplain #tZ}) parameters of a Molodensky transformation.</p>
*
+ * @deprecated This operation method is non-standard, of little use and has greater errors than intended.
+ *
* @author Martin Desruisseaux (Geomatys)
* @version 1.0
*
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-500">Deprecate (for removal) InterpolatedMolodenskyTransform</a>
* @see InterpolatedGeocentricTransform
*
* @since 0.7
*/
+@Deprecated(since="1.4", forRemoval=true)
public class InterpolatedMolodenskyTransform extends MolodenskyFormula {
/**
* Serial number for inter-operability with different versions.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform2D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform2D.java
index 8b341bc1ce..5cc19d8ce7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform2D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform2D.java
@@ -33,7 +33,10 @@ import org.apache.sis.referencing.datum.DatumShiftGrid;
* @author Martin Desruisseaux (Geomatys)
* @version 0.7
* @since 0.7
+ *
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-500">Deprecate (for removal) InterpolatedMolodenskyTransform</a>
*/
+@Deprecated(since="1.4", forRemoval=true)
final class InterpolatedMolodenskyTransform2D extends InterpolatedMolodenskyTransform implements MathTransform2D {
/**
* Serial number for compatibility with different versions.
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 b18873bf7d..c3faf5bdd7 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
@@ -54,7 +54,7 @@ import org.apache.sis.util.resources.Errors;
* @author Johann Sorel (Geomatys)
* @author Rémi Maréchal (Geomatys)
* @author Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.4
*
* @see MathTransforms#interpolate(double[], double[])
*
@@ -116,15 +116,15 @@ final class LinearInterpolator1D extends AbstractMathTransform1D implements Seri
* return a one-dimensional affine transform instead of an interpolator.
* We need to perform this check before the sign reversal applied after this loop.
*/
- double value;
+ double value, tolerance;
int i = 0;
do {
if (++i >= n) {
return LinearTransform1D.create(slope, offset);
}
value = values[i];
- } while (Numerics.epsilonEqual(value, Math.fma(i, slope, offset),
- Math.max(Math.abs(value), as) * Numerics.COMPARISON_THRESHOLD));
+ tolerance = Math.max(Math.abs(value), as) * Numerics.COMPARISON_THRESHOLD;
+ } while (Numerics.epsilonEqual(value, Math.fma(i, slope, offset), tolerance));
/*
* If the values are in decreasing order, reverse their sign so we get increasing order.
* We will multiply the results by -1 after the transformation.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform1D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform1D.java
index 72e8e9db50..18414ac790 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform1D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform1D.java
@@ -26,8 +26,11 @@ import org.opengis.referencing.operation.TransformException;
import org.opengis.referencing.operation.NoninvertibleTransformException;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.matrix.Matrix1;
-import org.apache.sis.referencing.operation.matrix.Matrix2;
+import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
import org.apache.sis.internal.referencing.provider.Affine;
+import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.internal.util.DoubleDouble;
+import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.ComparisonMode;
/*
* We really want to use doubleToRawLongBits, not doubleToLongBits, because the
@@ -46,14 +49,15 @@ import static java.lang.Double.doubleToRawLongBits;
* is faster. This kind of transform is extensively used by {@code org.apache.sis.coverage.grid.GridCoverage2D}.
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.7
+ * @version 1.4
*
* @see LogarithmicTransform1D
* @see ExponentialTransform1D
*
* @since 0.5
*/
-class LinearTransform1D extends AbstractMathTransform1D implements LinearTransform, Serializable {
+@SuppressWarnings("CloneInNonCloneableClass")
+class LinearTransform1D extends AbstractMathTransform1D implements LinearTransform, ExtendedPrecisionMatrix, Serializable {
/**
* Serial number for inter-operability with different versions.
*/
@@ -62,7 +66,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
/**
* A transform that just reverse the sign of input values.
*/
- static final LinearTransform1D NEGATE = new LinearTransform1D(-1, 0);
+ static final LinearTransform1D NEGATE = new LinearTransform1D(-1, 0, 0, 0);
/**
* The value which is multiplied to input values.
@@ -74,6 +78,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
*/
final double offset;
+ /**
+ * Error terms of {@link #scale} and {@link #offset} used in double-double arithmetic.
+ */
+ private final double scaleError, offsetError;
+
/**
* The inverse of this transform. Created only when first needed.
*/
@@ -84,14 +93,36 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* Instances should be created using the {@linkplain #create(double, double) factory method},
* which may returns optimized implementations for some particular argument values.
*
- * @param scale the {@code scale} term in the linear equation.
- * @param offset the {@code offset} term in the linear equation.
+ * @param scale the {@code scale} term in the linear equation.
+ * @param offset the {@code offset} term in the linear equation.
+ * @param scaleError error term of {@code scale} for double-double arithmetic.
+ * @param offsetError error term of {@code offset} for double-double arithmetic.
*
* @see #create(double, double)
*/
- protected LinearTransform1D(final double scale, final double offset) {
- this.scale = scale;
- this.offset = offset;
+ LinearTransform1D(final double scale, final double offset, final double scaleError, final double offsetError) {
+ this.scale = scale;
+ this.offset = offset;
+ this.scaleError = scaleError;
+ this.offsetError = offsetError;
+ }
+
+ /**
+ * Constructs a new linear transform and remember error terms.
+ *
+ * @param scale the {@code scale} term in the linear equation.
+ * @param offset the {@code offset} term in the linear equation.
+ * @return the linear transform for the given scale and offset.
+ */
+ static LinearTransform1D create(final DoubleDouble scale, final DoubleDouble offset) {
+ if (scale.error == 0) {
+ if (offset.error == 0) {
+ return create(scale.value, offset.value);
+ } else if (scale.value == 0) {
+ return new ConstantTransform1D(offset.value, offset.error);
+ }
+ }
+ return new LinearTransform1D(scale.value, offset.value, scale.error, offset.error);
}
/**
@@ -111,9 +142,9 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
if (scale == 0) {
if (offset == 0) return ConstantTransform1D.ZERO;
if (offset == 1) return ConstantTransform1D.ONE;
- return new ConstantTransform1D(offset);
+ return new ConstantTransform1D(offset, 0);
}
- return new LinearTransform1D(scale, offset);
+ return new LinearTransform1D(scale, offset, 0, 0);
}
/**
@@ -129,6 +160,43 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
return tr;
}
+ /**
+ * Implementation of Matrix API. No need for clone because this matrix is immutable.
+ */
+ @SuppressWarnings("CloneDoesntCallSuperClone")
+ @Override public final Matrix clone() {return this;}
+ @Override public final Matrix getMatrix() {return this;}
+ @Override public final int getNumRow() {return 2;}
+ @Override public final int getNumCol() {return 2;}
+
+ /** Unsupported operation because this matrix shall be immutable. */
+ @Override public void setElement(int row, int column, double value) {
+ throw new UnsupportedOperationException(Errors.format(Errors.Keys.UnmodifiableObject_1, Matrix.class));
+ }
+
+ /**
+ * Retrieves the value at the specified row and column of the matrix.
+ * Row and column indices can be only 0 or 1.
+ */
+ @Override
+ public final double getElement(final int row, final int column) {
+ if (((row | column) & ~1) != 0) {
+ throw new IndexOutOfBoundsException();
+ }
+ return (row == 0) ? ((column == 0) ? scale : offset) : column;
+ }
+
+ /**
+ * Returns a copy of all matrix elements followed by the error terms for extended-precision arithmetic.
+ * Matrix elements are returned in a flat, row-major (column indices vary fastest) array.
+ *
+ * @return a copy of matrix elements followed by error terms.
+ */
+ @Override
+ public final double[] getExtendedElements() {
+ return new double[] {scale, offset, 0, 1, scaleError, offsetError, 0, 0};
+ }
+
/**
* Returns the parameter descriptors for this math transform.
*/
@@ -149,19 +217,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
return Affine.parameters(getMatrix());
}
- /**
- * Returns this transform as an affine transform matrix.
- */
- @Override
- public Matrix getMatrix() {
- return new Matrix2(scale, offset, 0, 1);
- }
-
/**
* Creates the inverse transform of this object.
*/
@Override
- public LinearTransform1D inverse() throws NoninvertibleTransformException {
+ public synchronized LinearTransform1D inverse() throws NoninvertibleTransformException {
if (inverse == null) {
/*
* Note: we do not perform the following optimization, because MathTransforms.linear(…)
@@ -173,7 +233,15 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
*/
if (scale != 0) {
final LinearTransform1D inverse;
- inverse = create(1/scale, -offset/scale);
+ if (DoubleDouble.DISABLED) {
+ inverse = create(1/scale, -offset/scale);
+ } else {
+ final var sd = new DoubleDouble(scale, scaleError);
+ final var od = new DoubleDouble(sd);
+ od.inverseDivide(-offset, -offsetError);
+ sd.inverseDivide(1, 0);
+ inverse = create(sd, od);
+ }
inverse.inverse = this;
this.inverse = inverse;
} else {
@@ -187,20 +255,19 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* Returns {@code true} since this transform is affine.
*/
@Override
- public boolean isAffine() {
+ public final boolean isAffine() {
return true;
}
/**
* Tests whether this transform does not move any points.
- *
- * <div class="note"><b>Note:</b> this method should always returns {@code false}, since
+ * This method should always returns {@code false}, because
* {@code MathTransforms.linear(…)} should have created specialized implementations for identity cases.
* Nevertheless we perform the full check as a safety, in case someone instantiated this class directly
- * instead of using a factory method.</div>
+ * instead of using a factory method.
*/
@Override
- public boolean isIdentity() {
+ public final boolean isIdentity() {
return offset == 0 && scale == 1;
}
@@ -211,7 +278,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* @return the derivative at the given point.
*/
@Override
- public Matrix derivative(final DirectPosition point) throws TransformException {
+ public final Matrix derivative(final DirectPosition point) throws TransformException {
return new Matrix1(scale);
}
@@ -222,7 +289,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* @return the derivative at the given point.
*/
@Override
- public double derivative(final double value) {
+ public final double derivative(final double value) {
return scale;
}
@@ -230,8 +297,12 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* Transforms the specified value.
*/
@Override
- public double transform(double value) {
- return offset + scale * value;
+ public double transform(final double value) {
+ if (Formulas.USE_FMA) {
+ return Math.fma(value, scale, offset);
+ } else {
+ return offset + scale * value;
+ }
}
/**
@@ -245,7 +316,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
final boolean derivate)
{
if (dstPts != null) {
- dstPts[dstOff] = offset + scale*srcPts[srcOff];
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff] = Math.fma(srcPts[srcOff], scale, offset);
+ } else {
+ dstPts[dstOff] = offset + scale*srcPts[srcOff];
+ }
}
return derivate ? new Matrix1(scale) : null;
}
@@ -260,13 +335,21 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
{
if (srcPts != dstPts || srcOff >= dstOff) {
while (--numPts >= 0) {
- dstPts[dstOff++] = offset + scale * srcPts[srcOff++];
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff++] = Math.fma(srcPts[srcOff++], scale, offset);
+ } else {
+ dstPts[dstOff++] = offset + scale * srcPts[srcOff++];
+ }
}
} else {
srcOff += numPts;
dstOff += numPts;
while (--numPts >= 0) {
- dstPts[--dstOff] = offset + scale * srcPts[--srcOff];
+ if (Formulas.USE_FMA) {
+ dstPts[--dstOff] = Math.fma(srcPts[--srcOff], scale, offset);
+ } else {
+ dstPts[--dstOff] = offset + scale * srcPts[--srcOff];
+ }
}
}
}
@@ -282,13 +365,21 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
{
if (srcPts != dstPts || srcOff >= dstOff) {
while (--numPts >= 0) {
- dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]);
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff++] = (float) Math.fma(srcPts[srcOff++], scale, offset);
+ } else {
+ dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]);
+ }
}
} else {
srcOff += numPts;
dstOff += numPts;
while (--numPts >= 0) {
- dstPts[--dstOff] = (float) (offset + scale * srcPts[--srcOff]);
+ if (Formulas.USE_FMA) {
+ dstPts[--dstOff] = (float) Math.fma(srcPts[--srcOff], scale, offset);
+ } else {
+ dstPts[--dstOff] = (float) (offset + scale * srcPts[--srcOff]);
+ }
}
}
}
@@ -303,7 +394,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
final float [] dstPts, int dstOff, int numPts)
{
while (--numPts >= 0) {
- dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]);
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff++] = (float) Math.fma(srcPts[srcOff++], scale, offset);
+ } else {
+ dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]);
+ }
}
}
@@ -316,7 +411,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
final double[] dstPts, int dstOff, int numPts)
{
while (--numPts >= 0) {
- dstPts[dstOff++] = offset + scale * srcPts[srcOff++];
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff++] = Math.fma(srcPts[srcOff++], scale, offset);
+ } else {
+ dstPts[dstOff++] = offset + scale * srcPts[srcOff++];
+ }
}
}
@@ -348,7 +447,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
*/
@Override
protected int computeHashCode() {
- return Long.hashCode(doubleToRawLongBits(offset) + 31*doubleToRawLongBits(scale)) ^ super.computeHashCode();
+ return Long.hashCode(super.computeHashCode() ^
+ (doubleToRawLongBits(offset) +
+ doubleToRawLongBits(scale) * 31 +
+ doubleToRawLongBits(offsetError) * 37 +
+ doubleToRawLongBits(scaleError) * 7));
}
/**
@@ -365,8 +468,10 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
}
} else if (super.equals(object, mode)) {
final LinearTransform1D that = (LinearTransform1D) object;
- return doubleToRawLongBits(this.scale) == doubleToRawLongBits(that.scale) &&
- doubleToRawLongBits(this.offset) == doubleToRawLongBits(that.offset);
+ return doubleToRawLongBits(scale) == doubleToRawLongBits(that.scale) &&
+ doubleToRawLongBits(offset) == doubleToRawLongBits(that.offset) &&
+ doubleToRawLongBits(scaleError) == doubleToRawLongBits(that.scaleError) &&
+ doubleToRawLongBits(offsetError) == doubleToRawLongBits(that.offsetError);
/*
* NOTE: 'LinearTransform1D' and 'ConstantTransform1D' are heavily used by 'Category'
* from 'org.apache.sis.coverage' package. It is essential for Cateory to differenciate
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
index e9a8bae8a2..ac484de9ef 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
@@ -36,6 +36,7 @@ import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
import org.apache.sis.referencing.operation.matrix.AffineTransforms2D;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.matrix.Matrices;
+import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.util.OptionalCandidate;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.ArraysExt;
@@ -57,7 +58,7 @@ import org.apache.sis.util.Static;
* GeoAPI factory interfaces instead.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see MathTransformFactory
*
@@ -198,7 +199,10 @@ public final class MathTransforms extends Static {
if (Matrices.isAffine(matrix)) {
switch (sourceDimension) {
case 1: {
- return linear(matrix.getElement(0,0), matrix.getElement(0,1));
+ final MatrixSIS m = MatrixSIS.castOrCopy(matrix);
+ return LinearTransform1D.create(
+ DoubleDouble.castOrCopy(m.getNumber(0,0)),
+ DoubleDouble.castOrCopy(m.getNumber(0,1)));
}
case 2: {
if (matrix instanceof ExtendedPrecisionMatrix) {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/PoleRotation.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/PoleRotation.java
index 3fd77bdec3..a445b93e44 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/PoleRotation.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/PoleRotation.java
@@ -126,6 +126,7 @@ public class PoleRotation extends AbstractMathTransform2D implements Serializabl
*
* @see #inverse()
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private MathTransform2D inverse;
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
index 1ae2d841df..ed8506e906 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
@@ -23,6 +23,7 @@ import org.apache.sis.internal.referencing.DirectPositionView;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
+import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.util.ArgumentChecks;
@@ -35,7 +36,7 @@ import org.apache.sis.util.ArgumentChecks;
* lines in the source is preserved in the output.</p>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 1.2
+ * @version 1.4
*
* @see java.awt.geom.AffineTransform
*
@@ -281,7 +282,11 @@ 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.
*/
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
@@ -345,7 +350,11 @@ 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[], ...)
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
@@ -381,7 +390,11 @@ 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[], ...)
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
@@ -416,7 +429,11 @@ 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[], ...)
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java
index 0fe369f82f..e897d9feab 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java
@@ -663,7 +663,7 @@ public final class CoordinateOperationFinderTest extends MathTransformTestCase {
assertInstanceOf("operation", Conversion.class, operation);
transform = operation.getMathTransform();
- tolerance = 1E-12;
+ tolerance = 2E-12;
verifyTransform(new double[] {
// December 31, 1899 at 12:00 UTC in seconds.
CommonCRS.Temporal.DUBLIN_JULIAN.datum().getOrigin().getTime() / 1000
@@ -990,7 +990,7 @@ public final class CoordinateOperationFinderTest extends MathTransformTestCase {
0, 0, 0, 1
}), ((LinearTransform) transform).getMatrix(), 1E-12);
- tolerance = 1E-12;
+ tolerance = 2E-12;
verifyTransform(new double[] {
-5, -8, CommonCRS.Temporal.DUBLIN_JULIAN.datum().getOrigin().getTime() / 1000
}, new double[] {
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatricesTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatricesTest.java
index 64bc647065..58d455e271 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatricesTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatricesTest.java
@@ -41,7 +41,7 @@ import static org.opengis.referencing.cs.AxisDirection.*;
*
* @author Martin Desruisseaux (Geomatys)
* @author Johann Sorel (Geomatys)
- * @version 1.1
+ * @version 1.4
* @since 0.4
*/
@DependsOn({
@@ -106,6 +106,7 @@ public final class MatricesTest extends TestCase {
new AxisDirection[] {NORTH, EAST, UP},
new AxisDirection[] {NORTH, EAST, UP});
+ assertExtendedPrecision(matrix);
assertTrue ("isAffine", matrix.isAffine());
assertTrue ("isIdentity", matrix.isIdentity());
assertEquals("numRow", 4, matrix.getNumRow());
@@ -128,6 +129,7 @@ public final class MatricesTest extends TestCase {
new AxisDirection[] {NORTH, EAST, UP},
new AxisDirection[] {WEST, UP, SOUTH});
+ assertExtendedPrecision(matrix);
assertTrue ("isAffine", matrix.isAffine());
assertFalse("isIdentity", matrix.isIdentity());
assertEquals("numRow", 4, matrix.getNumRow());
@@ -155,6 +157,7 @@ public final class MatricesTest extends TestCase {
new AxisDirection[] {NORTH, EAST, UP},
new AxisDirection[] {DOWN, NORTH});
+ assertExtendedPrecision(matrix);
assertFalse("isIdentity", matrix.isIdentity());
assertEquals("numRow", 3, matrix.getNumRow());
assertEquals("numCol", 4, matrix.getNumCol());
@@ -180,6 +183,7 @@ public final class MatricesTest extends TestCase {
new AxisDirection[] {NORTH, EAST, UP},
new AxisDirection[] {DOWN, DOWN});
+ assertExtendedPrecision(matrix);
assertFalse("isIdentity", matrix.isIdentity());
assertEquals("numRow", 3, matrix.getNumRow());
assertEquals("numCol", 4, matrix.getNumCol());
@@ -242,6 +246,19 @@ public final class MatricesTest extends TestCase {
}
}
+ /**
+ * Asserts that the given matrix uses extended precision. This is mandatory for all matrices
+ * returned by {@link Matrices#createTransform(AxisDirection[], AxisDirection[])} because
+ * {@code CoordinateSystems.swapAndScaleAxes(…)} will modify those matrices in-place with
+ * the assumption that they accept extended precision.
+ *
+ * @param matrix the matrix to test.
+ */
+ private static void assertExtendedPrecision(final Matrix matrix) {
+ assertTrue(matrix instanceof GeneralMatrix);
+ assertTrue(((GeneralMatrix) matrix).isExtendedPrecision());
+ }
+
/**
* Tests {@link Matrices#createTransform(Envelope, Envelope)}.
* This method tests the example given in {@code Matrices.createTransform(…)} javadoc.
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 f4f2c9d68a..7c2e038a4c 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,14 +74,6 @@ 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/measure/AbstractConverter.java b/core/sis-utility/src/main/java/org/apache/sis/measure/AbstractConverter.java
index 6005d0aef9..e719e227ca 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/measure/AbstractConverter.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/measure/AbstractConverter.java
@@ -64,7 +64,7 @@ abstract class AbstractConverter implements UnitConverter, Serializable {
/**
* If the conversion can be represented by a polynomial equation, returns the coefficients of that equation.
- * Otherwise returns {@code null}.
+ * Otherwise returns {@code null}. This is the implementation of {@link Units#coefficients(UnitConverter)}.
*/
Number[] coefficients() {
return isIdentity() ? new Number[0] : null;
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 d75cc6884a..355b2b495a 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
@@ -26,6 +26,7 @@ import org.apache.sis.util.LenientComparable;
import org.apache.sis.math.MathFunctions;
import org.apache.sis.math.Fraction;
import org.apache.sis.internal.util.Numerics;
+import org.apache.sis.internal.util.DoubleDouble;
/**
@@ -239,6 +240,7 @@ final class LinearConverter extends AbstractConverter implements LenientComparab
/**
* Returns the coefficient of this linear conversion.
+ * Coefficients are offset and scale, in that order.
*/
@Override
@SuppressWarnings("fallthrough")
@@ -258,11 +260,13 @@ final class LinearConverter extends AbstractConverter implements LenientComparab
* package to perform more accurate calculations.
*/
private static Number ratio(final double value, final double divisor) {
- final int numerator = (int) value;
- if (numerator == value) {
- final int denominator = (int) divisor;
- if (denominator == divisor) {
- return (denominator == 1) ? numerator : new Fraction(numerator, denominator);
+ if (value != 0) {
+ final int numerator = (int) value;
+ if (numerator == value) {
+ final int denominator = (int) divisor;
+ if (denominator == divisor) {
+ return (denominator == 1) ? numerator : new Fraction(numerator, denominator);
+ }
}
}
return value / divisor;
@@ -278,15 +282,20 @@ final class LinearConverter extends AbstractConverter implements LenientComparab
/**
* Applies the linear conversion on the given value. This method uses {@link BigDecimal} arithmetic if
- * the given value is an instance of {@code BigDecimal}, or IEEE 754 floating-point arithmetic otherwise.
- *
- * <p>Apache SIS rarely uses {@link BigDecimal} arithmetic, so providing an efficient implementation of
- * this method is not a goal.</p>
+ * the given value is an instance of {@code BigDecimal}, or double-double arithmetic if the value is an
+ * instance of {@link DoubleDouble}, or IEEE 754 floating-point arithmetic otherwise.
*/
@Override
public Number convert(Number value) {
ArgumentChecks.ensureNonNull("value", value);
if (!isIdentity()) {
+ if (value instanceof DoubleDouble) {
+ var dd = new DoubleDouble((DoubleDouble) value);
+ dd.multiply(scale);
+ dd.add(offset);
+ dd.divide(divisor);
+ return dd;
+ }
if (value instanceof BigInteger) {
value = new BigDecimal((BigInteger) value);
}
diff --git a/core/sis-utility/src/main/java/org/apache/sis/measure/Units.java b/core/sis-utility/src/main/java/org/apache/sis/measure/Units.java
index 7327047740..1bec091089 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/measure/Units.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/measure/Units.java
@@ -1654,6 +1654,7 @@ public final class Units extends Static {
*
* @since 0.8
*/
+ @OptionalCandidate
@SuppressWarnings("fallthrough")
public static Number[] coefficients(final UnitConverter converter) {
if (converter != null) {
[sis] 01/02: Add more verifications of `GeneralMatrix` validity. Fix a bug in `GeneralMatrix.setElements(double[])` identified by above-cited verifications.
Posted by de...@apache.org.
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 4f11c0244de664dce82d95aebebd3eac72b9489a
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Fri Dec 30 14:14:38 2022 +0100
Add more verifications of `GeneralMatrix` validity.
Fix a bug in `GeneralMatrix.setElements(double[])`
identified by above-cited verifications.
---
.../operation/matrix/GeneralMatrix.java | 62 ++++++++++++++++++----
.../sis/referencing/operation/matrix/Matrices.java | 5 +-
.../referencing/operation/matrix/MatrixSIS.java | 10 +++-
.../referencing/operation/matrix/package-info.java | 2 +-
.../operation/matrix/GeneralMatrixTest.java | 5 +-
.../operation/matrix/MatrixTestCase.java | 27 +++++-----
.../org/apache/sis/internal/util/DoubleDouble.java | 14 +++--
.../apache/sis/internal/util/DoubleDoubleTest.java | 8 +--
8 files changed, 93 insertions(+), 40 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
index 90a4fcc596..0fba129a0c 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
@@ -36,7 +36,7 @@ import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
* the {@link DoubleDouble#error}.
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see Matrices#createDiagonal(int, int)
*
@@ -125,6 +125,7 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
this.numRow = (short) numRow;
this.numCol = (short) numCol;
this.elements = elements.clone();
+ assert isValid();
}
/**
@@ -145,6 +146,7 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
elements = new double[numRow * numCol];
getElements(matrix, numRow, numCol, elements);
}
+ assert isValid();
}
/**
@@ -154,6 +156,34 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
numRow = matrix.numRow;
numCol = matrix.numCol;
elements = matrix.elements.clone();
+ assert isValid();
+ }
+
+ /**
+ * Verifies that this matrix is well-formed. This method verifies that the matrix size is valid,
+ * that the {@link #elements} array is non-null and has a valid length, and that all error terms
+ * are smaller than 1 ULP of the corresponding matrix element value.
+ *
+ * @return whether this matrix is well-formed.
+ * @throws NullPointerException if the {@link #elements} array is null.
+ * @throws AssertionError (thrown by {@link DoubleDouble#DoubleDouble(double, double)} constructor)
+ * if an error term is not smaller than the corresponding matrix element value.
+ */
+ private boolean isValid() {
+ final int numRow = this.numRow;
+ final int numCol = this.numCol;
+ final int length = elements.length;
+ int i = numRow * numCol; // Cannot overflow.
+ if ((numRow | numCol) < 0 || (length != i && length != i*2) ||
+ (numRow != numCol) != (this instanceof NonSquareMatrix))
+ {
+ return false;
+ }
+ boolean isValid = true;
+ while (--i >= 0) {
+ isValid &= Numerics.equals(getNumber(i / numCol, i % numCol).doubleValue(), elements[i]);
+ }
+ return isValid;
}
/**
@@ -267,6 +297,7 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
* @param row the row index, from 0 inclusive to {@link #getNumRow()} exclusive.
* @param column the column index, from 0 inclusive to {@link #getNumCol()} exclusive.
* @return the current value at the given row and column, rounded to nearest integer.
+ * @throws ArithmeticException if the value is NaN or overflows integer capacity.
*
* @see DoubleDouble#longValue()
*
@@ -276,15 +307,18 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
public long getInteger(int row, int column) {
if (row >= 0 && row < numRow && column >= 0 && column < numCol) {
int i = row * numCol + column;
- long value = Math.round(elements[i]);
+ final double value = elements[i];
+ long r = Math.round(value);
i += numRow * numCol;
if (i < elements.length) {
- value += (long) elements[i]; // Really want rounding toward zero.
+ r += (long) elements[i]; // Really want rounding toward zero.
}
- return value;
- } else {
- throw indexOutOfBounds(row, column);
+ if (Math.abs(r - value) <= 0.5) {
+ return r;
+ }
+ throw new ArithmeticException(Errors.format(Errors.Keys.CanNotConvertValue_2, value, Long.TYPE));
}
+ throw indexOutOfBounds(row, column);
}
/**
@@ -395,7 +429,8 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
}
/**
- * {@inheritDoc}
+ * Returns a copy of all matrix elements in a flat, row-major (column indices vary fastest) array.
+ * The returned array does <em>not</em> include error terms used in double-double arithmetic.
*/
@Override
public final double[] getElements() {
@@ -412,15 +447,17 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
}
/**
- * {@inheritDoc}
+ * Sets all matrix elements from a flat, row-major (column indices vary fastest) array.
+ * The given array shall not contain error terms. The error terms will be set to default values.
*/
@Override
public final void setElements(final double[] newValues) {
ensureLengthMatch(numRow*numCol, newValues);
System.arraycopy(newValues, 0, elements, 0, newValues.length);
if (elements.length != newValues.length) {
- inferErrors(newValues);
+ inferErrors(elements);
}
+ assert isValid();
}
/**
@@ -436,7 +473,7 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
* </ul>
*
* @param newValues the new matrix elements in a row-major array.
- * @return {@code true} if at leat one {@link DoubleDouble} instance has been found, in which case all
+ * @return {@code true} if at least one {@link DoubleDouble} instance has been found, in which case all
* errors terms have been initialized, or {@code false} otherwise, in which case no error term
* has been initialized (this is a <cite>all or nothing</cite> operation).
* @throws IllegalArgumentException if the given array does not have the expected length.
@@ -486,6 +523,9 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
}
elements[i + length] = error;
}
+ if (isExtended) {
+ assert isValid();
+ }
return isExtended;
}
@@ -513,6 +553,7 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
} else {
super.setMatrix(matrix);
}
+ assert isValid();
}
/**
@@ -662,6 +703,7 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
sum.storeTo(elements, k++, errorOffset);
}
}
+ assert isValid();
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
index 6a74205416..c1f50e882c 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
@@ -212,8 +212,9 @@ public final class Matrices extends Static {
* Below is an intantionally undocumented feature. We use those sentinel values as a way to create
* matrices with extended precision without exposing our double-double arithmetic in public API.
*/
- if (elements == ExtendedPrecisionMatrix.IDENTITY || elements == ExtendedPrecisionMatrix.ZERO) {
- return GeneralMatrix.createExtendedPrecision(numRow, numCol, elements == ExtendedPrecisionMatrix.IDENTITY);
+ final boolean setToIdentity = (elements == ExtendedPrecisionMatrix.IDENTITY);
+ if (setToIdentity || elements == ExtendedPrecisionMatrix.ZERO) {
+ return GeneralMatrix.createExtendedPrecision(numRow, numCol, setToIdentity);
}
final GeneralMatrix matrix = GeneralMatrix.createExtendedPrecision(numRow, numCol, false);
if (matrix.setElements(elements)) {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
index c3d9f19e49..ab9afb41ec 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
@@ -55,7 +55,7 @@ import org.apache.sis.util.resources.Errors;
* </ul>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see Matrices
*
@@ -170,11 +170,17 @@ public abstract class MatrixSIS implements Matrix, LenientComparable, Cloneable,
* @param row the row index, from 0 inclusive to {@link #getNumRow()} exclusive.
* @param column the column index, from 0 inclusive to {@link #getNumCol()} exclusive.
* @return the current value at the given row and column, rounded to nearest integer.
+ * @throws ArithmeticException if the value is NaN or overflows integer capacity.
*
* @since 1.3
*/
public long getInteger(int row, int column) {
- return Math.round(getElement(row, column));
+ final double value = getElement(row, column);
+ final long r = Math.round(value);
+ if (Math.abs(r - value) <= 0.5) {
+ return r;
+ }
+ throw new ArithmeticException(Errors.format(Errors.Keys.CanNotConvertValue_2, value, Long.TYPE));
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/package-info.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/package-info.java
index b424dae044..ed8b445699 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/package-info.java
@@ -69,7 +69,7 @@
* Like SIS, Vecmath is optimized for small matrices of interest for 2D and 3D graphics.</p>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 1.3
+ * @version 1.4
* @since 0.4
*/
package org.apache.sis.referencing.operation.matrix;
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/GeneralMatrixTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/GeneralMatrixTest.java
index 7321fae51e..2febade1b6 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/GeneralMatrixTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/GeneralMatrixTest.java
@@ -66,7 +66,7 @@ public final class GeneralMatrixTest extends MatrixTestCase {
@Test
public void testExtendedPrecision() {
final long value = 1000000000000010000L;
- assertNotEquals(value, Math.round((double) value));
+ assertNotEquals(value, StrictMath.round((double) value));
final GeneralMatrix m = new GeneralMatrix(1, 1, false, 2);
final DoubleDouble ddval = new DoubleDouble((Long) value);
m.setNumber(0, 0, ddval);
@@ -155,7 +155,6 @@ public final class GeneralMatrixTest extends MatrixTestCase {
@Test
public void testTranslateVector() {
testTranslateVector(new GeneralMatrix(3, 3, true, 1)); // Double precision
-// testTranslateVector(new GeneralMatrix(3, 3, true, 2)); // Double-double precision
- // TODO: revisit commented-out test after using Math.fma with JDK9.
+ testTranslateVector(new GeneralMatrix(3, 3, true, 2)); // Double-double precision
}
}
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 c46f9dd025..120667f6c3 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
@@ -26,6 +26,7 @@ import org.apache.sis.test.TestUtilities;
import org.apache.sis.test.DependsOnMethod;
import org.junit.Test;
+import static java.lang.StrictMath.*;
import static org.apache.sis.test.Assert.*;
@@ -168,7 +169,7 @@ public abstract class MatrixTestCase extends TestCase {
assertEquals(name, e, actual.getNumber(j,i).doubleValue(), tolerance);
if (tolerance != STRICT && statistics != null) {
synchronized (statistics) {
- statistics.accept(StrictMath.abs(e - a));
+ statistics.accept(abs(e - a));
}
}
}
@@ -194,7 +195,7 @@ public abstract class MatrixTestCase extends TestCase {
private static void assertEqualsRelative(final String message, final double expected,
final MatrixSIS matrix, final int row, final int column)
{
- assertEquals(message, expected, matrix.getElement(row, column), StrictMath.abs(expected) * 1E-12);
+ assertEquals(message, expected, matrix.getElement(row, column), abs(expected) * 1E-12);
}
/**
@@ -204,7 +205,7 @@ public abstract class MatrixTestCase extends TestCase {
*/
private double nextNonZeroRandom() {
double value = random.nextDouble() * 200 - 100;
- value += StrictMath.copySign(0.001, value);
+ value += copySign(0.001, value);
if (random.nextBoolean()) {
value = 1 / value;
}
@@ -399,7 +400,7 @@ public abstract class MatrixTestCase extends TestCase {
final double e = matrix.getElement(j, i);
sum += e*e;
}
- return StrictMath.sqrt(sum);
+ return sqrt(sum);
}
/**
@@ -525,8 +526,8 @@ public abstract class MatrixTestCase extends TestCase {
if ((n % 10) == 0) {
setRandomValues(at, matrix);
}
- at.translate(vector[0] = Math.fma(random.nextDouble(), 50, -25),
- vector[1] = Math.fma(random.nextDouble(), 50, -25));
+ at.translate(vector[0] = fma(random.nextDouble(), 50, -25),
+ vector[1] = fma(random.nextDouble(), 50, -25));
matrix.translate(vector);
assertMatrixEquals("translate", AffineTransforms2D.toMatrix(at), matrix, TOLERANCE);
}
@@ -536,10 +537,10 @@ public abstract class MatrixTestCase extends TestCase {
* Sets random values in the given affine transform and a copy of those values in the given matrix.
*/
private void setRandomValues(final AffineTransform at, final MatrixSIS matrix) {
- at.setToRotation(random.nextDouble() * StrictMath.PI);
+ at.setToRotation(random.nextDouble() * PI);
at.scale(nextNonZeroRandom(), nextNonZeroRandom());
- at.translate(Math.fma(random.nextDouble(), 100, - 50),
- Math.fma(random.nextDouble(), 100, - 50));
+ at.translate(fma(random.nextDouble(), 100, - 50),
+ fma(random.nextDouble(), 100, - 50));
matrix.setElements(new double[] {
at.getScaleX(), at.getShearX(), at.getTranslateX(),
at.getShearY(), at.getScaleY(), at.getTranslateY(),
@@ -564,8 +565,8 @@ public abstract class MatrixTestCase extends TestCase {
if ((n % 10) == 0) {
setRandomValues(at, matrix);
}
- vector[0] = Math.fma(random.nextDouble(), 50, -25);
- vector[1] = Math.fma(random.nextDouble(), 50, -25);
+ vector[0] = fma(random.nextDouble(), 50, -25);
+ vector[1] = 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 +597,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] = Math.fma(random.nextDouble(), -10, 8);
+ elements[k] = fma(random.nextDouble(), -10, 8);
}
final Matrix referenceArg = new Matrix(elements, nx).transpose();
final MatrixSIS matrixArg = Matrices.create(numCol, nx, elements);
@@ -636,7 +637,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] = Math.fma(random.nextDouble(), -10, 8);
+ elements[k] = 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-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 af78f73076..6cf925f1ba 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
@@ -347,11 +347,15 @@ public final class DoubleDouble extends Number {
public static double errorForWellKnownValue(final double value) {
if (DISABLED) return 0;
final int i = Arrays.binarySearch(VALUES, Math.abs(value));
+ final double error;
if (i >= 0) {
- return MathFunctions.xorSign(ERRORS[i], value);
+ error = MathFunctions.xorSign(ERRORS[i], value);
+ } else {
+ final double delta = DecimalFunctions.deltaForDoubleToDecimal(value);
+ error = Double.isNaN(delta) ? 0 : delta;
}
- final double delta = DecimalFunctions.deltaForDoubleToDecimal(value);
- return Double.isNaN(delta) ? 0 : delta;
+ assert !(Math.abs(error) >= Math.ulp(value)) : value; // Use ! for being tolerant to NaN.
+ return error;
}
/**
@@ -834,8 +838,8 @@ public final class DoubleDouble extends Number {
final double thisValue = this.value;
final double thisError = this.error;
setToProduct(thisValue, otherValue);
- error += otherError * thisValue;
- error += otherValue * thisError;
+ error = Math.fma(otherError, thisValue, error);
+ error = Math.fma(otherValue, thisError, error);
normalize();
}
diff --git a/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java b/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
index 9b15b57586..f54ab22923 100644
--- a/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
+++ b/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
@@ -71,7 +71,7 @@ public final class DoubleDoubleTest extends TestCase {
* in order to change only the exponent part of the IEEE representation.
*/
private double nextRandom() {
- return random.nextDouble() * 2048 - 1024;
+ return fma(random.nextDouble(), 2048, -1024);
}
/**
@@ -224,7 +224,7 @@ public final class DoubleDoubleTest extends TestCase {
nextRandom(dd);
nextRandom(op);
final BigDecimal expected = toBigDecimal(dd).add(toBigDecimal(op));
- dd.add(op); // Must be after 'expected' computation.
+ dd.add(op); // Must be after `expected` computation.
assertExtendedEquals(expected, dd, ADD_TOLERANCE_FACTOR);
}
}
@@ -241,7 +241,7 @@ public final class DoubleDoubleTest extends TestCase {
nextRandom(dd);
nextRandom(op);
final BigDecimal expected = toBigDecimal(dd).multiply(toBigDecimal(op), MathContext.DECIMAL128);
- dd.multiply(op); // Must be after 'expected' computation.
+ dd.multiply(op); // Must be after `expected` computation.
assertExtendedEquals(expected, dd, PRODUCT_TOLERANCE_FACTOR);
}
}
@@ -258,7 +258,7 @@ public final class DoubleDoubleTest extends TestCase {
nextRandom(dd);
nextRandom(op);
final BigDecimal expected = toBigDecimal(dd).divide(toBigDecimal(op), MathContext.DECIMAL128);
- dd.divide(op); // Must be after 'expected' computation.
+ dd.divide(op); // Must be after `expected` computation.
assertExtendedEquals(expected, dd, PRODUCT_TOLERANCE_FACTOR);
}
}