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 2023/01/05 09:49:20 UTC
[sis] branch geoapi-4.0 updated: Reduce the use of double-double arithmetic where we do not expect significant improvement. Make explicit whether or not the numbers where intended to be exact in base 10.
This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
new d208d976b7 Reduce the use of double-double arithmetic where we do not expect significant improvement. Make explicit whether or not the numbers where intended to be exact in base 10.
d208d976b7 is described below
commit d208d976b7168a5c0e77e4295e714c4cddb23ef8
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Wed Jan 4 20:28:19 2023 +0100
Reduce the use of double-double arithmetic where we do not expect significant improvement.
Make explicit whether or not the numbers where intended to be exact in base 10.
---
.../org/apache/sis/coverage/grid/GridExtent.java | 4 +-
.../org/apache/sis/coverage/grid/GridGeometry.java | 10 +-
.../sis/coverage/grid/ResampledGridCoverage.java | 4 +-
.../main/java/org/apache/sis/portrayal/Canvas.java | 2 +-
.../referencing/provider/Equirectangular.java | 6 +-
.../sis/referencing/cs/CoordinateSystems.java | 9 +-
.../sis/referencing/datum/BursaWolfParameters.java | 16 +-
.../sis/referencing/datum/DatumShiftGrid.java | 2 +-
.../sis/referencing/datum/DefaultEllipsoid.java | 14 +-
.../sis/referencing/datum/TimeDependentBWP.java | 2 +-
.../operation/CoordinateOperationFinder.java | 4 +-
.../operation/matrix/GeneralMatrix.java | 4 +-
.../sis/referencing/operation/matrix/Matrices.java | 9 +-
.../referencing/operation/matrix/MatrixSIS.java | 18 +-
.../operation/projection/AlbersEqualArea.java | 24 +--
.../operation/projection/CylindricalEqualArea.java | 6 +-
.../operation/projection/Initializer.java | 70 +++----
.../projection/LambertConicConformal.java | 35 ++--
.../operation/projection/LongitudeWraparound.java | 4 +-
.../referencing/operation/projection/Mercator.java | 9 +-
.../operation/projection/MeridianArcBased.java | 8 +-
.../projection/ModifiedAzimuthalEquidistant.java | 10 +-
.../operation/projection/ObliqueStereographic.java | 4 +-
.../operation/projection/Orthographic.java | 6 +-
.../operation/projection/PolarStereographic.java | 18 +-
.../operation/projection/TransverseMercator.java | 4 +-
.../operation/transform/ContextualParameters.java | 2 +-
.../transform/EllipsoidToCentricTransform.java | 2 +-
.../operation/transform/MathTransforms.java | 4 +-
.../operation/projection/InitializerTest.java | 4 +-
.../org/apache/sis/internal/util/DoubleDouble.java | 201 +++++++++------------
.../src/main/java/org/apache/sis/math/Line.java | 8 +-
.../src/main/java/org/apache/sis/math/Plane.java | 18 +-
.../main/java/org/apache/sis/math/Statistics.java | 4 +-
.../org/apache/sis/measure/LinearConverter.java | 2 +-
.../apache/sis/internal/util/DoubleDoubleTest.java | 2 +-
.../sis/storage/geotiff/GridGeometryBuilder.java | 3 +-
37 files changed, 264 insertions(+), 288 deletions(-)
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
index 38b6b3b259..5d0ef798c0 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
@@ -1627,12 +1627,12 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable
final boolean flip = (flippedAxes & Numerics.bitmask(j)) != 0;
DoubleDouble offset = DoubleDouble.of(coordinates[i]);
DoubleDouble size = DoubleDouble.of(coordinates[i+srcDim]).subtract(offset).add(1);
- scale = DoubleDouble.of(env.getSpan(j)).divide(size);
+ scale = DoubleDouble.of(env.getSpan(j), true).divide(size);
if (flip) scale = scale.negate();
if (!offset.isZero()) { // Use `if` for keeping the value if scale is NaN.
offset = offset.multiply(scale).negate();
}
- offset = offset.add(flip ? env.getMaximum(j) : env.getMinimum(j));
+ offset = offset.add(flip ? env.getMaximum(j) : env.getMinimum(j), true);
affine.setNumber(j, srcDim, offset);
} else {
scale = DoubleDouble.NaN;
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
index 02d3fc6866..1e32bf653f 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
@@ -692,10 +692,10 @@ public class GridGeometry implements LenientComparable, Serializable {
resolution = new double[tgtDim];
for (int j=0; j<tgtDim; j++) {
final int i = (sourceDimensions != null) ? sourceDimensions[j] : j;
- DoubleDouble scale = DoubleDouble.of(affine.getNumber(j, i));
- DoubleDouble offset = DoubleDouble.of(affine.getNumber(j, srcDim));
+ DoubleDouble scale = DoubleDouble.of(affine.getNumber(j, i), true);
+ DoubleDouble offset = DoubleDouble.of(affine.getNumber(j, srcDim), true);
resolution[j] = Math.abs(scale.doubleValue());
- offset = offset.add(scale.multiply0(0.5));
+ offset = offset.add(scale.scalb(-1));
affine.setNumber(j, srcDim, offset);
}
gridToCRS = MathTransforms.linear(affine);
@@ -1422,9 +1422,9 @@ public class GridGeometry implements LenientComparable, Serializable {
final double p = periods[i];
if (p != 1) {
newResolution[i] /= p;
- final DoubleDouble pd = DoubleDouble.of(p);
+ final DoubleDouble pd = DoubleDouble.of(p, true);
for (int j=0; j<tgtDim; j++) {
- DoubleDouble e = DoubleDouble.of(matrix.getNumber(j, i));
+ DoubleDouble e = DoubleDouble.of(matrix.getNumber(j, i), true);
matrix.setNumber(j, i, e.divide(pd));
changed = true;
}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
index 611c4d360e..cb99e29e8c 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
@@ -401,14 +401,14 @@ final class ResampledGridCoverage extends DerivedGridCoverage {
vectors.setElement(tcDim, i, Double.NaN); // For preventing this row to be selected again.
}
DoubleDouble m = DoubleDouble.of(sign);
- m = m.divide(magnitudes.getNumber(0, tgDim));
+ m = m.divide(magnitudes.getNumber(0, tgDim), false);
crsToGrid.setNumber(tgDim, tcDim, m); // Scale factor from CRS coordinates to grid coordinates.
/*
* Move the point of interest in a place where conversion to source grid coordinates
* will be close to integer. The exact location does not matter; an additional shift
* will be applied later for translating to target grid extent.
*/
- m = m.multiply(originToPOI[tcDim] - targetPOI[tcDim]);
+ m = m.multiply(DoubleDouble.sum(originToPOI[tcDim], -targetPOI[tcDim]));
crsToGrid.setNumber(tgDim, crsDim, m);
}
crsToGrid.setElement(gridDim, crsDim, 1);
diff --git a/core/sis-portrayal/src/main/java/org/apache/sis/portrayal/Canvas.java b/core/sis-portrayal/src/main/java/org/apache/sis/portrayal/Canvas.java
index 23312c1e2f..e9188ab0ca 100644
--- a/core/sis-portrayal/src/main/java/org/apache/sis/portrayal/Canvas.java
+++ b/core/sis-portrayal/src/main/java/org/apache/sis/portrayal/Canvas.java
@@ -616,7 +616,7 @@ public class Canvas extends Observable implements Localized {
final int srcDim = magnitudes.getNumCol();
DoubleDouble scale = DoubleDouble.ZERO; // Will be set to average magnitude value.
for (int i=0; i<srcDim; i++) {
- scale = scale.add(magnitudes.getNumber(0, i));
+ scale = scale.add(magnitudes.getNumber(0, i), false);
}
scale = scale.divide(srcDim);
/*
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Equirectangular.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Equirectangular.java
index 56c076540a..4f1f4f35a1 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Equirectangular.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Equirectangular.java
@@ -351,10 +351,10 @@ public final class Equirectangular extends AbstractProvider {
final double sinφ1 = sin(φ1);
a = b / (1 - (1 - rs*rs) * (sinφ1*sinφ1));
}
- final DoubleDouble k = DoubleDouble.of(a);
+ final DoubleDouble k = DoubleDouble.of(a, true);
final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
- denormalize.convertAfter(0, k, DoubleDouble.of(fe));
- denormalize.convertAfter(1, k, DoubleDouble.of(fn));
+ denormalize.convertAfter(0, k, DoubleDouble.of(fe, true));
+ denormalize.convertAfter(1, k, DoubleDouble.of(fn, true));
/*
* Creates the ConcatenatedTransform, letting the factory returns the cached instance
* if the caller already invoked this method previously (which usually do not happen).
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/cs/CoordinateSystems.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/cs/CoordinateSystems.java
index 18cb2c22dd..3ab2ddb261 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/cs/CoordinateSystems.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/cs/CoordinateSystems.java
@@ -419,10 +419,11 @@ next: for (final CoordinateSystem cs : targets) {
default: throw new IncommensurableException(Resources.format(
Resources.Keys.NonLinearUnitConversion_2, sourceUnit, targetUnit));
}
- final var shift = matrix.getNumber(j, sourceDim);
- final var factor = DoubleDouble.of(matrix.getNumber(j,i));
- matrix.setNumber(j, i, factor.multiply(scale));
- matrix.setNumber(j, sourceDim, factor.multiply(offset).add(shift));
+ final boolean decimal = true; // Whether values were intended to be exact in base 10.
+ final var shift = DoubleDouble.of(matrix.getNumber(j, sourceDim), decimal);
+ final var factor = DoubleDouble.of(matrix.getNumber(j, i), decimal);
+ matrix.setNumber(j, i, factor.multiply(scale, decimal));
+ matrix.setNumber(j, sourceDim, factor.multiply(offset, decimal).add(shift));
}
}
return matrix;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/BursaWolfParameters.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/BursaWolfParameters.java
index fba3af963f..c5ef1a854a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/BursaWolfParameters.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/BursaWolfParameters.java
@@ -433,7 +433,7 @@ public class BursaWolfParameters extends FormattableObject implements Cloneable,
case 6: p = dS; break;
default: throw new AssertionError(index);
}
- return DoubleDouble.of(p);
+ return DoubleDouble.of(p, true);
}
/**
@@ -551,9 +551,9 @@ public class BursaWolfParameters extends FormattableObject implements Cloneable,
* elements should have the same value, but we tolerate slight deviation
* (this will be verified later).
*/
- DoubleDouble RS = DoubleDouble.of(getNumber(matrix, 0, 0))
- .add(getNumber(matrix, 1, 1))
- .add(getNumber(matrix, 2, 2)).divide(3);
+ DoubleDouble RS = DoubleDouble.of(getNumber(matrix, 0, 0), true)
+ .add(getNumber(matrix, 1, 1), true)
+ .add(getNumber(matrix, 2, 2), true).divide(3);
/*
* Computes: RS = S * toRadians(1″)
* dS = (S-1) * PPM
@@ -570,12 +570,12 @@ public class BursaWolfParameters extends FormattableObject implements Cloneable,
throw new IllegalArgumentException(Resources.format(Resources.Keys.NonUniformScale));
}
for (int i = j+1; i < SIZE-1; i++) {
- final DoubleDouble mr = DoubleDouble.of(getNumber(matrix, j, i)).divide(RS); // Negative rotation term.
- final DoubleDouble pr = DoubleDouble.of(getNumber(matrix, i, j)).divide(RS); // Positive rotation term.
- if (!(abs(pr.value + mr.value) <= 2*tolerance)) { // We expect mr ≈ -pr.
+ final DoubleDouble mr = DoubleDouble.of(getNumber(matrix, j, i), true).divide(RS); // Negative rotation term.
+ final DoubleDouble pr = DoubleDouble.of(getNumber(matrix, i, j), true).divide(RS); // Positive rotation term.
+ if (!(abs(pr.value + mr.value) <= 2*tolerance)) { // We expect mr ≈ -pr.
throw new IllegalArgumentException(Resources.format(Resources.Keys.NotASkewSymmetricMatrix));
}
- final double value = pr.subtract(mr).multiply0(0.5).doubleValue(); // Average of the two rotation terms.
+ final double value = pr.subtract(mr).scalb(-1).doubleValue(); // Average of the two rotation terms.
switch (j*SIZE + i) {
case 1: rZ = value; break;
case 2: rY = -value; break;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
index a12f47b114..f5b5f96883 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
@@ -707,7 +707,7 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
final int ny = gridSize[1];
for (int gridY=0; gridY<ny; gridY++) {
for (int gridX=0; gridX<nx; gridX++) {
- sum = sum.add0(getCellValue(dim, gridX, gridY));
+ sum = sum.add(getCellValue(dim, gridX, gridY), false);
}
}
return sum.doubleValue() / (nx * ny);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultEllipsoid.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultEllipsoid.java
index 5eab357cde..31ff7a650d 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultEllipsoid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultEllipsoid.java
@@ -436,7 +436,7 @@ public class DefaultEllipsoid extends AbstractIdentifiedObject implements Ellips
*/
private DoubleDouble eccentricitySquared() {
final DoubleDouble f = flattening(this);
- return f.multiply(2).subtract(f.square());
+ return f.scalb(1).subtract(f.square());
}
/**
@@ -452,11 +452,12 @@ public class DefaultEllipsoid extends AbstractIdentifiedObject implements Ellips
* </div>
*/
private static DoubleDouble flattening(final Ellipsoid e) {
+ final boolean decimal = true; // Parameters are presumed accurate in base 10 (not 2) by definition.
if (e.isIvfDefinitive()) {
- return DoubleDouble.ONE.divide(e.getInverseFlattening());
+ return DoubleDouble.ONE.divide(e.getInverseFlattening(), decimal);
} else {
- var a = DoubleDouble.of(e.getSemiMajorAxis()); // Presumed accurate in base 10 (not 2) by definition.
- return a.subtract(e.getSemiMinorAxis()).divide(a);
+ var a = DoubleDouble.of(e.getSemiMajorAxis(), decimal);
+ return a.subtract(e.getSemiMinorAxis(), decimal).divide(a);
}
}
@@ -522,8 +523,9 @@ public class DefaultEllipsoid extends AbstractIdentifiedObject implements Ellips
public double semiMajorAxisDifference(final Ellipsoid other) {
double semiMajor = other.getSemiMajorAxis();
semiMajor = other.getAxisUnit().getConverterTo(getAxisUnit()).convert(semiMajor); // Often a no-op.
- DoubleDouble a = DoubleDouble.of(semiMajor); // Presumed accurate in base 10 if no unit conversion.
- return a.subtract(getSemiMajorAxis()).doubleValue(); // Presumed accurate in base 10 (not 2) by definition.
+ // Presumed accurate in base 10 (not 2) by definition.
+ final DoubleDouble a = DoubleDouble.of(semiMajor, true);
+ return a.subtract(getSemiMajorAxis(), true).doubleValue();
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
index 3ae2c3bee1..c7350ac85d 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
@@ -190,7 +190,7 @@ public class TimeDependentBWP extends BursaWolfParameters {
case 6: d = ddS; break;
default: throw new AssertionError(index);
}
- p = p.add(factor.multiply(d));
+ p = p.add(factor.multiply(d, true));
}
return p;
}
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 c3b776b210..db15de88ac 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
@@ -874,7 +874,7 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
final Unit<Time> targetUnit = targetCS.getAxis(0).getUnit().asType(Time.class);
DoubleDouble epochShift = DoubleDouble.of(sourceDatum.getOrigin().getTime());
epochShift = epochShift.subtract(targetDatum.getOrigin().getTime());
- epochShift = DoubleDouble.of(Units.MILLISECOND.getConverterTo(targetUnit).convert(epochShift));
+ epochShift = DoubleDouble.of(Units.MILLISECOND.getConverterTo(targetUnit).convert(epochShift), true);
/*
* 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,
@@ -891,7 +891,7 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
throw new OperationNotFoundException(notFoundMessage(sourceCRS, targetCRS), exception);
}
final int translationColumn = matrix.getNumCol() - 1; // Paranoiac check: should always be 1.
- var translation = DoubleDouble.of(matrix.getNumber(0, translationColumn));
+ final var translation = DoubleDouble.of(matrix.getNumber(0, translationColumn), true);
matrix.setNumber(0, translationColumn, translation.add(epochShift));
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 d1ecadc85b..bcd6a6bc63 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
@@ -278,7 +278,7 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
assert dd.equals(getNumber(row, column));
return dd;
} else {
- return DoubleDouble.of(value);
+ return DoubleDouble.of(value, true);
}
}
@@ -498,7 +498,7 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
}
boolean isExtended = false;
for (int i=0; i<length; i++) {
- final DoubleDouble value = DoubleDouble.of(newValues[i]);
+ final DoubleDouble value = DoubleDouble.of(newValues[i], true);
final double element = value.doubleValue();
elements[i] = element;
final double error;
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 32c497b027..f62086b7a6 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
@@ -278,18 +278,19 @@ public final class Matrices extends Static {
*/
final boolean same = srcDir.equals(dstDir);
if (useEnvelopes) {
+ final boolean decimal = true; // Whether values are assumed exact in base 10.
/*
* See the comment in transform(Envelope, Envelope) for an explanation about why
* we use the lower/upper corners instead of getMinimum()/getMaximum() methods.
*/
DoubleDouble scale, translate;
- scale = DoubleDouble.of(dstEnvelope.getSpan(dstIndex))
- .divide(srcEnvelope.getSpan(srcIndex));
+ scale = DoubleDouble.of(dstEnvelope.getSpan(dstIndex), decimal)
+ .divide(srcEnvelope.getSpan(srcIndex), decimal);
if (!same) {
scale = scale.negate();
}
- translate = scale.multiply((same ? srcCorner : srcOppositeCorner).getOrdinate(srcIndex));
- translate = DoubleDouble.of(dstCorner.getOrdinate(dstIndex)).subtract(translate);
+ translate = scale.multiply((same ? srcCorner : srcOppositeCorner).getOrdinate(srcIndex), decimal);
+ translate = DoubleDouble.of(dstCorner.getOrdinate(dstIndex), decimal).subtract(translate);
matrix.setNumber(dstIndex, srcIndex, scale);
matrix.setNumber(dstIndex, srcAxes.length, translate);
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 d7e479d967..9f7a990c83 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
@@ -148,7 +148,7 @@ public abstract class MatrixSIS implements Matrix, LenientComparable, Cloneable,
* This method does not need to verify argument validity.
*/
DoubleDouble getDD(final int row, final int column) {
- return DoubleDouble.of(getElement(row, column));
+ return DoubleDouble.of(getElement(row, column), true);
}
/**
@@ -207,7 +207,7 @@ public abstract class MatrixSIS implements Matrix, LenientComparable, Cloneable,
* @since 0.8
*/
public void setNumber(int row, int column, final Number value) {
- final DoubleDouble n = DoubleDouble.of(value);
+ final DoubleDouble n = DoubleDouble.of(value, true);
if (n.error != 0) {
set(row, column, n); // Package-private method when public method can not be invoked.
} else {
@@ -509,17 +509,18 @@ public abstract class MatrixSIS implements Matrix, LenientComparable, Cloneable,
* @since 0.6
*/
public void convertBefore(final int srcDim, final Number scale, final Number offset) {
+ final boolean decimal = false; // Whether given values were intended to be exact in base 10.
final int lastCol = getNumCol() - 1;
ArgumentChecks.ensureValidIndex(lastCol, srcDim);
for (int j = getNumRow(); --j >= 0;) {
if (offset != null) {
DoubleDouble s = getDD(j, srcDim); // Scale factor
DoubleDouble t = getDD(j, lastCol); // Translation factor
- set(j, lastCol, t.add(s.multiply(offset)));
+ set(j, lastCol, t.add(s.multiply(offset, decimal)));
}
if (scale != null) {
DoubleDouble s = getDD(j, srcDim); // Scale factor
- set(j, srcDim, s.multiply(scale));
+ set(j, srcDim, s.multiply(scale, decimal));
}
}
}
@@ -539,18 +540,19 @@ public abstract class MatrixSIS implements Matrix, LenientComparable, Cloneable,
* @since 0.6
*/
public void convertAfter(final int tgtDim, final Number scale, final Number offset) {
+ final boolean decimal = false; // Whether given values were intended to be exact in base 10.
final int lastRow = getNumRow() - 1;
final int lastCol = getNumCol() - 1;
ArgumentChecks.ensureValidIndex(lastRow, tgtDim);
if (scale != null) {
for (int i=lastCol; i>=0; i--) {
DoubleDouble s = getDD(tgtDim, i);
- set(tgtDim, i, s.multiply(scale));
+ set(tgtDim, i, s.multiply(scale, decimal));
}
}
if (offset != null) {
DoubleDouble t = getDD(tgtDim, lastCol);
- set(tgtDim, lastCol, t.add(offset));
+ set(tgtDim, lastCol, t.add(offset, decimal));
}
}
@@ -605,7 +607,7 @@ public abstract class MatrixSIS implements Matrix, LenientComparable, Cloneable,
for (int j=0; j<target.length; j++) {
var sum = DoubleDouble.ZERO;
for (int i=0; i<numCol; i++) {
- sum = sum.add(getDD(j, i).multiply(vector[i]));
+ sum = sum.add(getDD(j, i).multiply(vector[i], true));
}
target[j] = sum.doubleValue();
}
@@ -648,7 +650,7 @@ public abstract class MatrixSIS implements Matrix, LenientComparable, Cloneable,
for (int j=0; j<numRow; j++) {
var sum = DoubleDouble.ZERO;
for (int i=0; i<numCol; i++) {
- sum = sum.add(getDD(j, i).multiply(vector[i]));
+ sum = sum.add(getDD(j, i).multiply(vector[i], true));
}
set(j, numCol-1, sum);
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
index 986a27e6e9..8133f8bc5d 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
@@ -127,29 +127,29 @@ public class AlbersEqualArea extends AuthalicConversion {
final double cosφ1 = cos(φ1);
final double sinφ2 = sin(φ2);
final double cosφ2 = cos(φ2);
- final var m1sq = initializer.scaleAtφ(sinφ1, cosφ1).square(); // (cos(φ₁) / √(1 – ℯ²sin²φ₁))²
- final double α1 = qm(sinφ1); // Omitted ×(1-ℯ²) (see below)
+ final double m1 = initializer.scaleAtφ(sinφ1, cosφ1); // (cos(φ₁) / √(1 – ℯ²sin²φ₁))²
+ final double α1 = qm(sinφ1); // Omitted ×(1-ℯ²) (see below)
if (secant) {
- var m2sq = initializer.scaleAtφ(sinφ2, cosφ2).square(); // (cos(φ₂) / √(1 – ℯ²sin²φ₂))²
- double α2 = qm(sinφ2); // Omitted ×(1-ℯ²) (see below)
- nm = m1sq.subtract(m2sq).divide0(α2 - α1).doubleValue(); // n = nm / (1-ℯ²)
+ final double m2 = initializer.scaleAtφ(sinφ2, cosφ2); // (cos(φ₂) / √(1 – ℯ²sin²φ₂))²
+ final double α2 = qm(sinφ2); // Omitted ×(1-ℯ²) (see below)
+ nm = (m1*m1 - m2*m2) / (α2 - α1); // n = nm / (1-ℯ²)
} else {
nm = sinφ1;
}
- C = m1sq.add0(nm*α1).doubleValue(); // Omitted (1-ℯ²) term in nm cancels with omitted (1-ℯ²) term in α₁.
+ C = m1*m1 + nm*α1; // Omitted (1-ℯ²) term in nm cancels with omitted (1-ℯ²) term in α₁.
/*
* Compute rn = (1-ℯ²)/nm, which is the reciprocal of the "real" n used in Snyder and EPSG guidance note.
* We opportunistically use double-double arithmetic since the MatrixSIS operations use them anyway, but
- * we do not really have that accuracy because of the limited precision of 'nm'. The intent is rather to
- * increase the chances term cancellations happen during concatenation of coordinate operations.
+ * we do not really have that accuracy because of the limited precision of `nm`. The intent is rather to
+ * increase the chances that term cancellations happen during concatenation of coordinate operations.
*/
- final DoubleDouble rn = DoubleDouble.ONE.subtract(initializer.eccentricitySquared).divide0(nm);
+ final DoubleDouble rn = DoubleDouble.ONE.subtract(initializer.eccentricitySquared).divide(nm, false);
/*
- * Compute ρ₀ = √(C - n⋅q(sinφ₀))/n with multiplication by a omitted because already taken in account
- * by the denormalization matrix. Omitted (1-ℯ²) term in nm cancels with omitted (1-ℯ²) term in qm(…).
+ * Compute ρ₀ = √(C - n⋅q(sinφ₀))/n with multiplication by `a` omitted because already taken in account
+ * by the denormalization matrix. Omitted (1-ℯ²) term in `nm` cancels with omitted (1-ℯ²) term in qm(…).
* See above note about double-double arithmetic usage.
*/
- final var ρ0 = DoubleDouble.of0(C - nm*qm(sinφ0)).sqrt().multiply(rn);
+ final var ρ0 = DoubleDouble.sum(C, -nm*qm(sinφ0)).sqrt().multiply(rn);
/*
* At this point, all parameters have been processed. Now process to their
* validation and the initialization of (de)normalize affine transforms.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
index 5dc53165ea..58128e97ed 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
@@ -166,8 +166,8 @@ public class CylindricalEqualArea extends AuthalicConversion {
* but we nevertheless support it.
*/
final double φ1 = toRadians(initializer.getAndStore(STANDARD_PARALLEL));
- final DoubleDouble k0 = initializer.scaleAtφ(sin(φ1), cos(φ1))
- .multiply(initializer.getAndStore(Mercator1SP.SCALE_FACTOR));
+ final DoubleDouble k0 = DoubleDouble.of(initializer.scaleAtφ(sin(φ1), cos(φ1)), false)
+ .multiply(initializer.getAndStore(Mercator1SP.SCALE_FACTOR), true);
/*
* In most Apache SIS map projection implementations, the scale factor is handled by the super-class by
* specifying a ParameterRole.SCALE_FACTOR. However, in the case of this CylindricalEqualArea we rather
@@ -181,7 +181,7 @@ public class CylindricalEqualArea extends AuthalicConversion {
*/
DoubleDouble ik;
ik = DoubleDouble.ONE.subtract(initializer.eccentricitySquared);
- ik = ik.multiply0(0.5); // This line need to be cancelled when using spherical formulas.
+ ik = ik.scalb(-1); // This line need to be cancelled when using spherical formulas.
ik = ik.divide(k0);
final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
denormalize.convertAfter(0, k0, null);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
index 29269e30f1..6fe0baa6b1 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
@@ -128,11 +128,11 @@ final class Initializer {
final double fn = getAndStore(roles.get(ParameterRole.FALSE_NORTHING))
- getAndStore(roles.get(ParameterRole.FALSE_SOUTHING));
- DoubleDouble k = DoubleDouble.of(a); // The value by which to multiply all results of normalized projection.
+ DoubleDouble k = DoubleDouble.of(a, true); // The value by which to multiply all results of normalized projection.
if (a == b) {
eccentricitySquared = DoubleDouble.ZERO;
} else if (variant != null && variant.useAuthalicRadius()) {
- k = DoubleDouble.of(Formulas.getAuthalicRadius(a, b));
+ k = DoubleDouble.of(Formulas.getAuthalicRadius(a, b), false);
eccentricitySquared = DoubleDouble.ZERO;
} else {
/*
@@ -161,10 +161,10 @@ final class Initializer {
* constructor applies corrections for making those values more accurate in base 10 rather than 2.
*/
if (isIvfDefinitive) {
- var f = DoubleDouble.of(parameters.parameter(Constants.INVERSE_FLATTENING).doubleValue()).inverse();
- eccentricitySquared = f.multiply(2).subtract(f.square());
+ var f = DoubleDouble.of(parameters.parameter(Constants.INVERSE_FLATTENING).doubleValue(), true).inverse();
+ eccentricitySquared = f.scalb(1).subtract(f.square());
} else {
- eccentricitySquared = DoubleDouble.ONE.subtract(DoubleDouble.of(b).divide(k).square());
+ eccentricitySquared = DoubleDouble.ONE.subtract(DoubleDouble.of(b, true).divide(k).square());
}
final ParameterDescriptor<? extends Number> radius = roles.get(ParameterRole.LATITUDE_OF_CONFORMAL_SPHERE_RADIUS);
if (radius != null) {
@@ -185,7 +185,7 @@ final class Initializer {
* k = b / (1 - eccentricitySquared * (sinφ*sinφ));
*/
k = rν2(sin(toRadians(parameters.doubleValue(radius))));
- k = DoubleDouble.of(b).divide(k);
+ k = DoubleDouble.of(b, true).divide(k);
}
}
/*
@@ -195,7 +195,7 @@ final class Initializer {
*/
final ParameterDescriptor<? extends Number> scaleFactor = roles.get(ParameterRole.SCALE_FACTOR);
if (scaleFactor != null) {
- k = k.multiply(getAndStore(scaleFactor));
+ k = k.multiply(getAndStore(scaleFactor), true);
}
/*
* Set meridian rotation, scale factor, false easting and false northing parameter values
@@ -203,8 +203,8 @@ final class Initializer {
*/
context.normalizeGeographicInputs(λ0);
final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
- denormalize.convertAfter(0, k, DoubleDouble.of(fe));
- denormalize.convertAfter(1, k, DoubleDouble.of(fn));
+ denormalize.convertAfter(0, k, DoubleDouble.of(fe, true));
+ denormalize.convertAfter(1, k, DoubleDouble.of(fn, true));
}
/**
@@ -276,6 +276,13 @@ final class Initializer {
return DoubleDouble.ONE.subtract(eccentricitySquared).sqrt();
}
+ /**
+ * Returns the radius of a hypothetical sphere having the same surface than the ellipsoid.
+ */
+ final double authalicRadius() {
+ return Formulas.getAuthalicRadius(1, axisLengthRatio().doubleValue());
+ }
+
/**
* Computes the square of the reciprocal of the radius of curvature of the ellipsoid
* perpendicular to the meridian at latitude φ. That radius of curvature is:
@@ -283,37 +290,32 @@ final class Initializer {
* <blockquote>ν = 1 / √(1 - ℯ²⋅sin²φ)</blockquote>
*
* This method returns 1/ν², which is the (1 - ℯ²⋅sin²φ) part of above equation.
- * Special cases:
+ * In other words, the radius of curvature is 1 / √[rν²(sinφ)].
+ *
+ * <h4>Special cases</h4>
* <ul>
* <li>If φ is 0°, then <var>m</var> is 1.</li>
* <li>If φ is ±90°, then <var>m</var> is 0 provided that we are not in the spherical case
* (otherwise we get {@link Double#NaN}).</li>
* </ul>
*
+ * <h4>Accuracy</h4>
+ * The accuracy is only a few digits more than {@code double} value, not the full double-double precision.
+ * This extra accuracy come from the (1 - <var>small value</var>) part of the equation.
+ *
* @param sinφ the sine of the φ latitude.
* @return reciprocal squared of the radius of curvature of the ellipsoid
* perpendicular to the meridian at latitude φ.
*/
final DoubleDouble rν2(final double sinφ) {
if (DoubleDouble.DISABLED) {
- return DoubleDouble.of0(1 - eccentricitySquared.doubleValue() * (sinφ*sinφ));
+ return DoubleDouble.of(1 - eccentricitySquared.doubleValue() * (sinφ*sinφ), false);
}
/*
- * Compute 1 - ℯ²⋅sin²φ. Since ℯ²⋅sin²φ may be small,
+ * Compute 1 - ℯ²⋅sin²φ. Because ℯ²⋅sin²φ may be small,
* this is where double-double arithmetic has more value.
*/
- return DoubleDouble.ONE.subtract(eccentricitySquared.multiply(DoubleDouble.of0(sinφ).square()));
- }
-
- /**
- * Returns the radius of curvature of the ellipsoid perpendicular to the meridian at latitude φ.
- * This is {@code 1/sqrt(rν2(sinφ))}.
- *
- * @param sinφ the sine of the φ latitude.
- * @return radius of curvature of the ellipsoid perpendicular to the meridian at latitude φ.
- */
- final double radiusOfCurvature(final double sinφ) {
- return rν2(sinφ).sqrt().inverse().doubleValue();
+ return DoubleDouble.ONE.subtract(eccentricitySquared.multiply(DoubleDouble.product(sinφ, sinφ)));
}
/**
@@ -327,27 +329,31 @@ final class Initializer {
* the use of φ₀ (or φ₁ as relevant to method) for φ is suggested, except if the projection is
* equal area when the radius of authalic sphere should be used.
*
+ * <h4>Accuracy</h4>
+ * The accuracy is only a few digits more than {@code double} value, not the full double-double precision.
+ * This extra accuracy come from the (1 - <var>small value</var>) parts in both numerator and denominator.
+ *
* @param sinφ the sine of the φ latitude.
* @return radius of the conformal sphere at latitude φ.
*/
- final double radiusOfConformalSphere(final double sinφ) {
- return DoubleDouble.ONE.subtract(eccentricitySquared).sqrt().divide(rν2(sinφ)).doubleValue();
+ final DoubleDouble radiusOfConformalSphere(final double sinφ) {
+ return DoubleDouble.ONE.subtract(eccentricitySquared).sqrt().divide(rν2(sinφ));
}
/**
* Returns the scale factor at latitude φ (Snyder 14-15). This is computed as:
*
- * <blockquote>cosφ / sqrt(rν2(sinφ))</blockquote>
+ * <blockquote>cosφ / √[rν²(sinφ)]</blockquote>
*
- * The result is returned as a {@code DoubleDouble} for allowing callers to perform some additional
- * operations on it, but the {@link DoubleDouble#error} term should be considered mostly meaningless
- * because of the limited precision of {@code sinφ} and {@code cosφ} terms.
+ * <h4>Accuracy</h4>
+ * The result is returned as a {@code double} because the limited precision of {@code sinφ} and {@code cosφ}
+ * makes the error term meaningless. We use double-double arithmetic only for intermediate calculation.
*
* @param sinφ the sine of the φ latitude.
* @param cosφ the cosine of the φ latitude.
* @return scale factor at latitude φ.
*/
- final DoubleDouble scaleAtφ(final double sinφ, final double cosφ) {
- return DoubleDouble.of0(cosφ).divide(rν2(sinφ).sqrt());
+ final double scaleAtφ(final double sinφ, final double cosφ) {
+ return DoubleDouble.of(cosφ, false).divide(rν2(sinφ).sqrt()).doubleValue();
}
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
index c8a1b90af6..52ebbce184 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
@@ -292,7 +292,7 @@ public class LambertConicConformal extends ConformalProjection {
* a 0/0 indetermination.
*/
final double sinφ1 = sin(φ1);
- final var m1 = initializer.scaleAtφ(sinφ1, cos(φ1));
+ final double m1 = initializer.scaleAtφ(sinφ1, cos(φ1));
final double t1 = expΨ(φ1, eccentricity*sinφ1);
/*
* Compute n = (ln m₁ – ln m₂) / (ln t₁ – ln t₂), which we rewrite as ln(m₁/m₂) / ln(t₁/t₂)
@@ -301,9 +301,9 @@ public class LambertConicConformal extends ConformalProjection {
*/
if (abs(φ1 - φ2) >= ANGULAR_TOLERANCE) { // Should be `true` for 2SP case.
final double sinφ2 = sin(φ2);
- final var m2 = initializer.scaleAtφ(sinφ2, cos(φ2));
+ final double m2 = initializer.scaleAtφ(sinφ2, cos(φ2));
final double t2 = expΨ(φ2, eccentricity*sinφ2);
- n = log(m1.divide(m2).doubleValue()) / log(t1/t2); // Tend toward 0/0 if φ1 ≈ φ2.
+ n = log(m1/m2) / log(t1/t2); // Tend toward 0/0 if φ1 ≈ φ2.
} else {
n = -sinφ1;
}
@@ -311,12 +311,10 @@ public class LambertConicConformal extends ConformalProjection {
* Compute F = m₁/(n⋅t₁ⁿ) from Geomatics Guidance Note number 7.
* Following constants will be stored in the denormalization matrix, to be applied
* after the non-linear formulas implemented by this LambertConicConformal class.
- * Opportunistically use double-double arithmetic since the matrix coefficients will
- * be stored in that format anyway. This makes a change in the 2 or 3 last digits.
*/
- DoubleDouble F = m1.divide(DoubleDouble.of0(n).multiply0(pow(t1, n)));
+ double F = m1 / (n * pow(t1, n));
if (!isNorth) {
- F = F.negate();
+ F = -F;
}
/*
* Compute the radius of the parallel of latitude of the false origin.
@@ -326,9 +324,9 @@ public class LambertConicConformal extends ConformalProjection {
*
* EPSG uses this term in the computation of y = FN + rF – r⋅cos(θ).
*/
- DoubleDouble rF = null;
+ Number rF = null;
if (φ0 != copySign(PI/2, -n)) { // For reducing the rounding error documented in expΨ(+π/2).
- rF = F.multiply0(pow(expΨ(φ0, eccentricity*sin(φ0)), n));
+ rF = F * pow(expΨ(φ0, eccentricity*sin(φ0)), n);
}
/*
* At this point, all parameters have been processed. Now store
@@ -347,20 +345,19 @@ public class LambertConicConformal extends ConformalProjection {
* - Multiply by the scale factor (done by the super-class constructor).
* - Add false easting and false northing (done by the super-class constructor).
*/
- DoubleDouble sλ = DoubleDouble.of0(n);
- DoubleDouble sφ = null;
- if (isNorth) {
- // Reverse the sign of either longitude or latitude values before map projection.
- sφ = DoubleDouble.of(-1);
+ double sλ = n;
+ Number sφ = null;
+ if (isNorth) { // Reverse the sign of either longitude or latitude values before map projection.
+ sφ = -1;
} else {
- sλ = sλ.negate();
+ sλ = -sλ;
}
final MatrixSIS normalize = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
- normalize .convertAfter(0, sλ, (variant == Variant.BELGIUM) ? BELGE_A : null);
- normalize .convertAfter(1, sφ, null);
- denormalize.convertBefore(0, F, null);
- denormalize.convertBefore(1, F.negate(), rF);
+ normalize .convertAfter (0, sλ, (variant == Variant.BELGIUM) ? BELGE_A : null);
+ normalize .convertAfter (1, sφ, null);
+ denormalize.convertBefore(0, F, null);
+ denormalize.convertBefore(1, -F, rF);
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LongitudeWraparound.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LongitudeWraparound.java
index fff0f4e4e4..e28be4ffee 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LongitudeWraparound.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LongitudeWraparound.java
@@ -130,8 +130,8 @@ final class LongitudeWraparound extends AbstractMathTransform2D implements Seria
* @return a bound of the [−n⋅π … n⋅π] range.
*/
static double boundOfScaledLongitude(final MatrixSIS normalize, final boolean negative) {
- DoubleDouble bound = DoubleDouble.of(normalize.getNumber(0, 0));
- bound = bound.multiply0(negative ? Longitude.MIN_VALUE : Longitude.MAX_VALUE);
+ DoubleDouble bound = DoubleDouble.of(normalize.getNumber(0, 0), true);
+ bound = bound.multiply(negative ? Longitude.MIN_VALUE : Longitude.MAX_VALUE, false);
return bound.doubleValue();
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
index 551824e591..3e0a60247f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
@@ -34,7 +34,6 @@ import org.apache.sis.internal.referencing.provider.MercatorSpherical;
import org.apache.sis.internal.referencing.provider.MercatorAuxiliarySphere;
import org.apache.sis.internal.referencing.provider.RegionalMercator;
import org.apache.sis.internal.referencing.provider.PseudoMercator;
-import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.referencing.operation.matrix.Matrix2;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
@@ -268,7 +267,7 @@ public class Mercator extends ConformalProjection {
denormalize.convertBefore(0, k0, null);
denormalize.convertBefore(1, k0, null);
if (φ0 != 0) {
- denormalize.convertBefore(1, null, DoubleDouble.of0(-log(expΨ(φ0, eccentricity * sin(φ0)))));
+ denormalize.convertBefore(1, null, -log(expΨ(φ0, eccentricity * sin(φ0))));
}
/*
* Variants of the Mercator projection which can be handled by scale factors.
@@ -284,7 +283,7 @@ public class Mercator extends ConformalProjection {
normalize .convertBefore(1, 0.80, null);
denormalize.convertBefore(1, 1.25, null);
} else if (variant == Variant.AUXILIARY) {
- DoubleDouble ratio = null;
+ final Number ratio;
final int type = initializer.getAndStore(MercatorAuxiliarySphere.AUXILIARY_SPHERE_TYPE, 0);
switch (type) {
default: {
@@ -292,9 +291,9 @@ public class Mercator extends ConformalProjection {
MercatorAuxiliarySphere.AUXILIARY_SPHERE_TYPE.getName().getCode(), type));
}
case AuthalicMercator.TYPE:
- case 2: ratio = DoubleDouble.of0(Formulas.getAuthalicRadius(1, initializer.axisLengthRatio().doubleValue())); break;
+ case 2: ratio = initializer.authalicRadius(); break;
case 1: ratio = initializer.axisLengthRatio(); break;
- case 0: break; // Same as "Popular Visualisation Pseudo Mercator".
+ case 0: ratio = null; break; // Same as "Popular Visualisation Pseudo Mercator".
}
denormalize.convertAfter(0, ratio, null);
denormalize.convertAfter(1, ratio, null);
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 a74b0be97e..65c18a19da 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,7 +16,6 @@
*/
package org.apache.sis.referencing.operation.projection;
-import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.internal.referencing.Formulas;
import static java.lang.Math.*;
@@ -137,9 +136,10 @@ abstract class MeridianArcBased extends NormalizedProjection {
*/
rµ = fma(e6, -5./256,
fma(e4, -3./64,
- fma(e2, -1./4, 1))); // Part of Snyder 7-19 for computing rectifying latitude.
- DoubleDouble e1 = initializer.axisLengthRatio().ratio_1m_1p();
- final double ei = e1.doubleValue(); // Equivalent to [1 - √(1 - ℯ²)] / [1 + √(1 - ℯ²)] (Snyder 3-24).
+ fma(e2, -1./4, 1))); // Part of Snyder 7-19 for computing rectifying latitude.
+
+ // Equivalent to [1 - √(1 - ℯ²)] / [1 + √(1 - ℯ²)] (Snyder 3-24).
+ final double ei = initializer.axisLengthRatio().ratio_1m_1p().doubleValue();
final double ei2 = ei*ei;
final double ei3 = ei*ei2;
final double ei4 = ei2*ei2;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
index 83aabd1529..a2199f0c06 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
@@ -45,7 +45,7 @@ import static org.apache.sis.internal.referencing.provider.ModifiedAzimuthalEqui
*
* @author Martin Desruisseaux (Geomatys)
* @author Maxime Gavens (Geomatys)
- * @version 1.1
+ * @version 1.4
*
* @see AzimuthalEquidistant
*
@@ -118,11 +118,9 @@ public class ModifiedAzimuthalEquidistant extends AzimuthalEquidistant {
@Workaround(library="JDK", version="1.8")
private ModifiedAzimuthalEquidistant(final Initializer initializer) {
super(initializer);
- final double axisRatio, ν0, f;
- axisRatio = initializer.axisLengthRatio().doubleValue();
- ν0 = initializer.radiusOfCurvature(sinφ0);
- ℯ2_ν0_sinφ0 = eccentricitySquared * ν0 * sinφ0;
- f = eccentricity / axisRatio; // √(1 - ℯ²) = b/a
+ var ν0 = initializer.rν2(sinφ0).sqrt().inverse();
+ ℯ2_ν0_sinφ0 = initializer.eccentricitySquared.multiply(ν0).doubleValue() * sinφ0;
+ double f = eccentricity / initializer.axisLengthRatio().doubleValue(); // √(1 - ℯ²) = b/a
G = f * sinφ0;
Hp = f * cosφ0;
Bp = 3*eccentricitySquared * (sinφ0*cosφ0) / (1 - eccentricitySquared);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
index 21f718433c..7fb8643462 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
@@ -63,7 +63,7 @@ import static org.apache.sis.internal.referencing.provider.ObliqueStereographic.
* @author Rémi Maréchal (Geomatys)
* @author Martin Desruisseaux (Geomatys)
* @author Emmanuel Giasson (Thales)
- * @version 1.3
+ * @version 1.4
* @since 0.7
*/
public class ObliqueStereographic extends NormalizedProjection {
@@ -171,7 +171,7 @@ public class ObliqueStereographic extends NormalizedProjection {
* by 2 times the radius of the conformal sphere. Since this is a linear operation, we combine it with other
* linear operations performed by the denormalization matrix.
*/
- final double R2 = 2 * initializer.radiusOfConformalSphere(sinφ0);
+ final var R2 = initializer.radiusOfConformalSphere(sinφ0).scalb(1);
denormalize.convertBefore(0, R2, null);
denormalize.convertBefore(1, R2, null);
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Orthographic.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Orthographic.java
index 0b28ba9298..459fb5bc87 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Orthographic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Orthographic.java
@@ -47,7 +47,7 @@ import static org.apache.sis.internal.referencing.provider.Orthographic.*;
* @author Rueben Schulz (UBC)
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.1
+ * @version 1.4
* @since 1.1
*/
public class Orthographic extends NormalizedProjection {
@@ -115,9 +115,9 @@ public class Orthographic extends NormalizedProjection {
* since the computed value below is zero in the spherical case.
*/
if (eccentricity != 0) {
- final double ν0 = initializer.radiusOfCurvature(sinφ0);
+ final double ν0_cosφ0 = initializer.scaleAtφ(sinφ0, cosφ0);
final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
- denormalize.convertBefore(1, null, eccentricitySquared * ν0 * (sinφ0 * cosφ0));
+ denormalize.convertBefore(1, null, eccentricitySquared * ν0_cosφ0 * sinφ0);
}
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
index ed9ab81178..f63414f30a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
@@ -218,8 +218,8 @@ public class PolarStereographic extends ConformalProjection {
*/
φ1 = -φ1;
}
- φ1 = toRadians(φ1); // May be anything in [-π/2 … 0] range.
- final Number ρ, ρF; // This ρF is actually -ρF in EPSG guide.
+ φ1 = toRadians(φ1); // May be anything in [-π/2 … 0] range.
+ final Number ρ, ρF; // This ρF is actually -ρF in EPSG guide.
if (abs(φ1 + PI/2) < ANGULAR_TOLERANCE) {
/*
* Polar Stereographic (variant A)
@@ -228,13 +228,13 @@ public class PolarStereographic extends ConformalProjection {
* ρ = 2⋅a⋅k₀⋅t / √[(1+ℯ)^(1+ℯ) ⋅ (1–ℯ)^(1–ℯ)]
*
* In this implementation, we omit:
- * - the 'a' and 'k₀' factors, because they are handled outside this class,
- * - the 't' factor, because it needs to be computed in the transform(…) method.
+ * - the `a` and `k₀` factors, because they are handled outside this class,
+ * - the `t` factor, because it needs to be computed in the transform(…) method.
*
* In the spherical case, should give ρ == 2.
*/
- ρ = DoubleDouble.of0(2 / sqrt(pow(1+eccentricity, 1+eccentricity)
- * pow(1-eccentricity, 1-eccentricity)));
+ ρ = 2 / sqrt(pow(1+eccentricity, 1+eccentricity)
+ * pow(1-eccentricity, 1-eccentricity));
ρF = null;
} else {
/*
@@ -255,9 +255,9 @@ public class PolarStereographic extends ConformalProjection {
* In the spherical case, should give ρ = 1 + sinφ1 (Snyder 21-7 and 21-11).
*/
final double sinφ1 = sin(φ1);
- final DoubleDouble mF = initializer.scaleAtφ(sinφ1, cos(φ1));
- ρ = mF.divide0(expΨ(φ1, eccentricity*sinφ1));
- ρF = (variant == Variant.C) ? mF.negate() : null;
+ final double mF = initializer.scaleAtφ(sinφ1, cos(φ1));
+ ρ = mF / expΨ(φ1, eccentricity*sinφ1);
+ ρF = (variant == Variant.C) ? -mF : null;
}
/*
* At this point, all parameters have been processed. Now process to their
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
index 1ec07c0d89..860e57643f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
@@ -286,11 +286,11 @@ public class TransverseMercator extends NormalizedProjection {
*/
final double Q = asinh(tan(φ0)) - eccentricity * atanh(eccentricity * sin(φ0));
final double β = atan(sinh(Q));
- DoubleDouble M0 = B.negate().multiply0(
+ DoubleDouble M0 = B.negate().multiply(
β + fma(cf2, sin(2*β),
fma(cf4, sin(4*β),
fma(cf6, sin(6*β),
- cf8* sin(8*β)))));
+ cf8* sin(8*β)))), false);
/*
* At this point, all parameters have been processed. Now store
* the linear operations in the (de)normalize affine transforms:
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 a979b1cebe..823375db7d 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
@@ -467,7 +467,7 @@ public class ContextualParameters extends Parameters implements Serializable {
final DoubleDouble toRadians = DoubleDouble.DEGREES_TO_RADIANS;
DoubleDouble offset = null;
if (λ0 != 0) {
- offset = DoubleDouble.of(-λ0).multiply(toRadians);
+ offset = DoubleDouble.of(-λ0, true).multiply(toRadians);
}
final MatrixSIS normalize = (MatrixSIS) this.normalize; // Must be the same instance, not a copy.
normalize.convertBefore(0, toRadians, offset);
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 caaa88cb0c..e38b769da0 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
@@ -287,7 +287,7 @@ public class EllipsoidToCentricTransform extends AbstractMathTransform implement
* - A "denormalization" transform for scaling (X,Y,Z) to the semi-major axis length.
*/
context.normalizeGeographicInputs(0);
- final DoubleDouble a = DoubleDouble.of(semiMajor);
+ final DoubleDouble a = DoubleDouble.of(semiMajor, true);
final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
for (int i=0; i<3; i++) {
denormalize.convertAfter(i, a, null);
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 179dc11510..dc1d283fe9 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
@@ -201,8 +201,8 @@ public final class MathTransforms extends Static {
case 1: {
final MatrixSIS m = MatrixSIS.castOrCopy(matrix);
return LinearTransform1D.create(
- DoubleDouble.of(m.getNumber(0,0)),
- DoubleDouble.of(m.getNumber(0,1)));
+ DoubleDouble.of(m.getNumber(0,0), true),
+ DoubleDouble.of(m.getNumber(0,1), true));
}
case 2: {
if (matrix instanceof ExtendedPrecisionMatrix) {
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java
index 6d313296a0..9b4d2b7176 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/InitializerTest.java
@@ -82,8 +82,8 @@ public final class InitializerTest extends TestCase{
* use the Initializer class.
*/
final double φ0 = toRadians(initializer.getAndStore(ObliqueStereographic.LATITUDE_OF_ORIGIN));
+ final double rc = initializer.radiusOfConformalSphere(sin(φ0)).doubleValue();
assertTrue(φ0 > 0);
- assertEquals("Conformal Sphere Radius", 6382644.571, 6377397.155 *
- initializer.radiusOfConformalSphere(sin(φ0)), Formulas.LINEAR_TOLERANCE);
+ assertEquals("Conformal Sphere Radius", 6382644.571, 6377397.155 * rc, Formulas.LINEAR_TOLERANCE);
}
}
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 2582187f1f..cbb6508ca0 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
@@ -215,32 +215,38 @@ public final class DoubleDouble extends Number {
* {@code DoubleDouble}, {@link BigDecimal}, {@link BigInteger} or {@link Fraction},
* then the error term will be taken in account.
*
- * @param n the value.
+ * @param value the value.
+ * @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the value as a double-double number.
*/
- public static DoubleDouble of(Number n) {
- if (n instanceof DoubleDouble) {
- return (DoubleDouble) n;
+ public static DoubleDouble of(Number value, final boolean decimal) {
+ if (value instanceof DoubleDouble) {
+ return (DoubleDouble) value;
}
- if (n instanceof Fraction) {
- final Fraction f = (Fraction) n;
+ if (value instanceof Fraction) {
+ final Fraction f = (Fraction) value;
return new DoubleDouble(f.numerator, 0)
.divide(f.denominator);
}
- final double value, error;
- if (n instanceof BigInteger) {
- n = new BigDecimal((BigInteger) n, MathContext.DECIMAL128);
+ final double v, error;
+ if (value instanceof Float) {
+ final float f = (Float) value;
+ value = decimal ? DecimalFunctions.floatToDouble(f) : f;
+ } else if (value instanceof BigInteger) {
+ value = new BigDecimal((BigInteger) value, MathContext.DECIMAL128);
}
- value = n.doubleValue();
- if (n instanceof Long) {
- error = n.longValue() - (long) value; // Need rounding toward zero.
- } else if (n instanceof BigDecimal) {
+ v = value.doubleValue();
+ if (value instanceof Long) {
+ error = value.longValue() - (long) v; // Need rounding toward zero.
+ } else if (value instanceof BigDecimal) {
// Really need new BigDecimal(value) below, not BigDecimal.valueOf(value).
- error = ((BigDecimal) n).subtract(new BigDecimal(value), MathContext.DECIMAL64).doubleValue();
+ error = ((BigDecimal) value).subtract(new BigDecimal(v), MathContext.DECIMAL64).doubleValue();
+ } else if (decimal) {
+ error = errorForWellKnownValue(v);
} else {
- error = errorForWellKnownValue(value);
+ error = 0;
}
- return new DoubleDouble(value, error);
+ return new DoubleDouble(v, error);
}
/**
@@ -265,30 +271,18 @@ public final class DoubleDouble extends Number {
}
/**
- * CReturns an instance for the given value and an error term inferred by
- * {@link #errorForWellKnownValue(double)}.
- *
- * <p><b>Tip:</b> if the other value is known to be an integer or a power of 2,
- * then invoking {@link #of0(double)} is more efficient.</p>
- *
- * @param value the initial value.
- * @return the value as a double-double number.
- */
- public static DoubleDouble of(final double value) {
- return new DoubleDouble(value, errorForWellKnownValue(value));
- }
-
- /**
- * Returns an instance for the given value verbatim, without inferring an error term for double-double arithmetic.
- * Use this method when the value has been computed using transcendental functions (cosine, logarithm, <i>etc.</i>)
- * in which case there is no way we can infer a meaningful error term.
- * Also used when the value is known to have an exact representation as a {@code double} primitive type.
+ * Returns an instance for the given value.
+ * If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
+ * {@code decimal} should be {@code false} when the value has been computed using transcendental functions
+ * (cosine, logarithm, <i>etc.</i>), in which case there is no way we can infer a meaningful error term.
+ * Should also be {@code false} (for performance reason) when the value is an exact representation in base 2.
*
- * @param value the value to wrap in a {@code DoubleDouble} instance.
+ * @param value the value.
+ * @param decimal whether the value was intended to be exact in base 10.
* @return the value as a double-double number.
*/
- public static DoubleDouble of0(final double value) {
- return new DoubleDouble(value, 0);
+ public static DoubleDouble of(final double value, final boolean decimal) {
+ return new DoubleDouble(value, decimal ? errorForWellKnownValue(value) : 0);
}
/**
@@ -340,8 +334,7 @@ public final class DoubleDouble extends Number {
}
/**
- * Return value + error.
- * The result should be identical to {@link #value},
+ * Return value + error. The result should be identical to {@link #value},
* but we nevertheless do the sum as a safety.
*
* @return {@link #value} + {@link #error}.
@@ -415,6 +408,8 @@ public final class DoubleDouble extends Number {
/**
* Returns the sum of the given numbers.
+ * The double-double accuracy is useful when the largest value is exact in base 2.
+ * A typical example is: 1 - (small value).
*
* <p>Source: [Hida & al.] page 4 algorithm 4, itself reproduced from [Shewchuk] page 314.</p>
*
@@ -430,6 +425,8 @@ public final class DoubleDouble extends Number {
/**
* Returns the product of the given numbers.
+ * Note that unless the given arguments are exact in base 2,
+ * the result is not more accurate than a {@code double}.
*
* <p>Source: [Hida & al.] page 5 algorithm 7, itself reproduced from [Shewchuk] page 326.
* This is the algorithm used when FMA instruction is available. For an algorithm without FMA,
@@ -472,11 +469,12 @@ public final class DoubleDouble extends Number {
* Adds a {@code Number} value to this {@code DoubleDouble}. If the given number is an instance
* of {@code DoubleDouble} or {@link Fraction}, then the error term will be taken in account.
*
- * @param other the other value to add to this {@code DoubleDouble}.
+ * @param other the other value to add to this {@code DoubleDouble}.
+ * @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the sum of {@code this} with the given number.
*/
- public DoubleDouble add(final Number other) {
- return add(of(other));
+ public DoubleDouble add(final Number other, final boolean decimal) {
+ return add(of(other, decimal));
}
/**
@@ -500,26 +498,15 @@ public final class DoubleDouble extends Number {
}
/**
- * Adds a {@code double} value to this {@code DoubleDouble} with a default error term.
- *
- * <p><b>Tip:</b> if the other value is known to be an integer or a power of 2,
- * then invoking {@link #add0(double)} is more efficient.</p>
+ * Adds a {@code double} value to this {@code DoubleDouble}.
+ * If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
*
- * @param other the other value to add to this {@code DoubleDouble}.
+ * @param other the other value to add to this {@code DoubleDouble}.
+ * @param decimal whether the value was intended to be exact in base 10.
* @return the sum of {@code this} with the given number.
*/
- public DoubleDouble add(final double other) {
- return add(other, errorForWellKnownValue(other));
- }
-
- /**
- * Adds a {@code double} value to this {@code DoubleDouble} without error term.
- *
- * @param other the other value to add to this {@code DoubleDouble}.
- * @return the sum of {@code this} with the given number.
- */
- public DoubleDouble add0(final double other) {
- return add(other, 0);
+ public DoubleDouble add(final double other, final boolean decimal) {
+ return add(other, decimal ? errorForWellKnownValue(other) : 0);
}
/**
@@ -584,11 +571,12 @@ public final class DoubleDouble extends Number {
* Subtracts a {@code Number} from this {@code DoubleDouble}. If the given number is an instance
* of {@code DoubleDouble} or {@link Fraction}, then the error term will be taken in account.
*
- * @param other the other value to subtract from this {@code DoubleDouble}.
+ * @param other the other value to subtract from this {@code DoubleDouble}.
+ * @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the difference between {@code this} and the given number.
*/
- public DoubleDouble subtract(final Number other) {
- return subtract(of(other));
+ public DoubleDouble subtract(final Number other, final boolean decimal) {
+ return subtract(of(other, decimal));
}
/**
@@ -613,26 +601,15 @@ public final class DoubleDouble extends Number {
/**
* Subtracts a {@code double} from this {@code DoubleDouble} with a default error term.
+ * If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
*
- * <p><b>Tip:</b> if the other value is known to be an integer or a power of 2,
- * then invoking {@link #subtract0(double)} is more efficient.</p>
- *
- * @param other the other value to subtract from this {@code DoubleDouble}.
+ * @param other the other value to subtract from this {@code DoubleDouble}.
+ * @param decimal whether the value was intended to be exact in base 10.
* @return the difference between {@code this} and the given number.
*/
- public DoubleDouble subtract(double other) {
+ public DoubleDouble subtract(double other, final boolean decimal) {
other = -other;
- return add(other, errorForWellKnownValue(other));
- }
-
- /**
- * Subtracts a {@code double} from this {@code DoubleDouble} without error term.
- *
- * @param other the other value to subtract from this {@code DoubleDouble}.
- * @return the difference between {@code this} and the given number.
- */
- public DoubleDouble subtract0(final double other) {
- return add(-other, 0);
+ return add(other, decimal ? errorForWellKnownValue(other) : 0);
}
/**
@@ -649,11 +626,12 @@ public final class DoubleDouble extends Number {
* Multiplies this {@code DoubleDouble} by a {@code Number}. If the given number is an instance
* of {@code DoubleDouble} or {@link Fraction}, then the error term will be taken in account.
*
- * @param other the other value to multiply by this {@code DoubleDouble}.
+ * @param other the other value to multiply by this {@code DoubleDouble}.
+ * @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the product of {@code this} with the given number.
*/
- public DoubleDouble multiply(final Number other) {
- return multiply(of(other));
+ public DoubleDouble multiply(final Number other, final boolean decimal) {
+ return multiply(of(other, decimal));
}
/**
@@ -678,25 +656,14 @@ public final class DoubleDouble extends Number {
/**
* Multiplies this {@code DoubleDouble} by a {@code double} with a default error term.
+ * If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
*
- * <p><b>Tip:</b> if the other value is known to be an integer or a power of 2,
- * then invoking {@link #multiply0(double)} is more efficient.</p>
- *
- * @param other the other value to multiply by this {@code DoubleDouble}.
+ * @param other the other value to multiply by this {@code DoubleDouble}.
+ * @param decimal whether the value was intended to be exact in base 10.
* @return the product of {@code this} with the given number.
*/
- public DoubleDouble multiply(final double other) {
- return multiply(other, errorForWellKnownValue(other));
- }
-
- /**
- * Multiplies this {@code DoubleDouble} by a {@code double} without error term.
- *
- * @param other the other value to multiply by this {@code DoubleDouble}.
- * @return the product of {@code this} with the given number.
- */
- public DoubleDouble multiply0(final double other) {
- return multiply(other, 0);
+ public DoubleDouble multiply(final double other, final boolean decimal) {
+ return multiply(other, decimal ? errorForWellKnownValue(other) : 0);
}
/**
@@ -759,11 +726,12 @@ public final class DoubleDouble extends Number {
* Divides this {@code DoubleDouble} by a {@code Number}. If the given number is an instance
* of {@code DoubleDouble} or {@link Fraction}, then the error term will be taken in account.
*
- * @param other the other value by which to divide this {@code DoubleDouble}.
+ * @param other the other value by which to divide this {@code DoubleDouble}.
+ * @param decimal whether {@code float} and {@code double} values were intended to be exact in base 10.
* @return the ratio between {@code this} and the given number.
*/
- public DoubleDouble divide(final Number other) {
- return divide(of(other));
+ public DoubleDouble divide(final Number other, final boolean decimal) {
+ return divide(of(other, decimal));
}
/**
@@ -788,25 +756,14 @@ public final class DoubleDouble extends Number {
/**
* Divides this {@code DoubleDouble} by a {@code double} with a default error term.
+ * If {@code decimal} is {@code true}, then an error term is inferred for well-known values.
*
- * <p><b>Tip:</b> if the other value is known to be an integer or a power of 2,
- * then invoking {@link #divide0(double)} is more efficient.</p>
- *
- * @param other the other value by which to divide this {@code DoubleDouble}.
+ * @param other the other value by which to divide this {@code DoubleDouble}.
+ * @param decimal whether the value was intended to be exact in base 10.
* @return the ratio between {@code this} and the given number.
*/
- public DoubleDouble divide(final double other) {
- return divide(other, errorForWellKnownValue(other));
- }
-
- /**
- * Divides this {@code DoubleDouble} by a {@code double} without error term.
- *
- * @param other the other value by which to divide this {@code DoubleDouble}.
- * @return the ratio between {@code this} and the given number.
- */
- public DoubleDouble divide0(final double other) {
- return divide(other, 0);
+ public DoubleDouble divide(final double other, final boolean decimal) {
+ return divide(other, decimal ? errorForWellKnownValue(other) : 0);
}
/**
@@ -864,6 +821,18 @@ public final class DoubleDouble extends Number {
return ONE.subtract(this).divide(add(1));
}
+ /**
+ * Returns {@code this} × 2ⁿ. Typical usages are
+ * {@code scalb(1)} for an efficient multiplication by 2 and
+ * {@code scalb(-1)} for an efficient division by 2.
+ *
+ * @param n power of 2 used to scale {@code this}.
+ * @return {@code this} × 2ⁿ.
+ */
+ public DoubleDouble scalb(final int n) {
+ return new DoubleDouble(Math.scalb(value, n), Math.scalb(error, n));
+ }
+
/**
* Computes the square of this value.
*
diff --git a/core/sis-utility/src/main/java/org/apache/sis/math/Line.java b/core/sis-utility/src/main/java/org/apache/sis/math/Line.java
index 369e5bb440..dd6db9020e 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/math/Line.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/math/Line.java
@@ -318,8 +318,8 @@ public class Line implements DoubleUnaryOperator, Cloneable, Serializable {
if (!isNaN(y = p.getOrdinate(1)) && // Test first the dimension which is most likely to contain NaN.
!isNaN(x = p.getOrdinate(0)))
{
- mean_x = mean_x.add0(x);
- mean_y = mean_y.add0(y);
+ mean_x = mean_x.add(x, false);
+ mean_y = mean_y.add(y, false);
n++;
}
}
@@ -346,9 +346,9 @@ public class Line implements DoubleUnaryOperator, Cloneable, Serializable {
if (!isNaN(y = p.getOrdinate(1)) && // Test first the dimension which is most likely to contain NaN.
!isNaN(x = p.getOrdinate(0)))
{
- var dx = DoubleDouble.of(x).subtract(mean_x); // Δx = x - mean_x
+ var dx = DoubleDouble.of(x, true).subtract(mean_x); // Δx = x - mean_x
mean_x2 = mean_x2.add(dx.square()); // mean_x² += (Δx)²
- mean_xy = mean_xy.add(dx.multiply(y)); // mean_xy += Δx * y
+ mean_xy = mean_xy.add(dx.multiply(y, true)); // mean_xy += Δx * y
mean_y2 = mean_y2.add(DoubleDouble.product(y, y)); // mean_y² += y²
}
}
diff --git a/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java b/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java
index 1e2407bf92..eb1831bf7c 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/math/Plane.java
@@ -410,9 +410,9 @@ public class Plane implements DoubleBinaryOperator, Cloneable, Serializable {
final double x = p.getOrdinate(0); if (Double.isNaN(x)) continue;
final double y = p.getOrdinate(1); if (Double.isNaN(y)) continue;
final double z = p.getOrdinate(2); if (Double.isNaN(z)) continue;
- sum_x = sum_x .add0(x);
- sum_y = sum_y .add0(y);
- sum_z = sum_z .add0(z);
+ sum_x = sum_x .add(x, false);
+ sum_y = sum_y .add(y, false);
+ sum_z = sum_z .add(z, false);
sum_xx = sum_xx.add(DoubleDouble.product(x, x));
sum_yy = sum_yy.add(DoubleDouble.product(y, y));
sum_xy = sum_xy.add(DoubleDouble.product(x, y));
@@ -445,17 +445,17 @@ public class Plane implements DoubleBinaryOperator, Cloneable, Serializable {
if (Double.isNaN(z)) {
throw new IllegalArgumentException(Errors.format(Errors.Keys.NotANumber_1, Strings.toIndexed("z", n)));
}
- sum_z = sum_z .add0(z);
+ sum_z = sum_z .add(z, false);
sum_zx = sum_zx.add(DoubleDouble.product(z, x));
sum_zy = sum_zy.add(DoubleDouble.product(z, y));
n++;
}
}
- sum_x = DoubleDouble.of0(n/2d).multiply (nx-1); // Division by 2 is exact.
- sum_y = DoubleDouble.of0(n/2d).multiply (ny-1);
- sum_xx = DoubleDouble.of (n) .multiply0(nx-0.5).multiply(nx-1).divide(3);
- sum_yy = DoubleDouble.of (n) .multiply0(ny-0.5).multiply(ny-1).divide(3);
- sum_xy = DoubleDouble.of0(n/4d).multiply (ny-1) .multiply(nx-1);
+ sum_x = DoubleDouble.of(n/2d, false).multiply(nx-1); // Division by 2 is exact.
+ sum_y = DoubleDouble.of(n/2d, false).multiply(ny-1);
+ sum_xx = DoubleDouble.of(n) .multiply(nx-0.5, false).multiply(nx-1).divide(3);
+ sum_yy = DoubleDouble.of(n) .multiply(ny-0.5, false).multiply(ny-1).divide(3);
+ sum_xy = DoubleDouble.of(n/4d, false).multiply(ny-1) .multiply(nx-1);
this.n = n;
}
diff --git a/core/sis-utility/src/main/java/org/apache/sis/math/Statistics.java b/core/sis-utility/src/main/java/org/apache/sis/math/Statistics.java
index 3c01123e48..9a792fb0aa 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/math/Statistics.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/math/Statistics.java
@@ -210,13 +210,13 @@ public class Statistics implements DoubleConsumer, LongConsumer, Cloneable, Seri
if (!Double.isNaN(mean) || !Double.isNaN(standardDeviation)) {
ArgumentChecks.ensureBetween("mean", minimum, maximum, mean);
}
- final DoubleDouble sd = DoubleDouble.of(mean).multiply(count);
+ final DoubleDouble sd = DoubleDouble.of(mean, true).multiply(count);
sum = sd.value;
lowBits = sd.error;
/*
* squareSum = standardDeviation² × (allPopulation ? count : count-1) + sum²/count
*/
- DoubleDouble sq = DoubleDouble.of(standardDeviation);
+ DoubleDouble sq = DoubleDouble.of(standardDeviation, true);
sq = sq.square().multiply(allPopulation ? count : count-1);
sq = sq.add(sd.square().divide(count));
squareSum = sq.value;
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 1770f77a64..607f714f5d 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
@@ -291,7 +291,7 @@ final class LinearConverter extends AbstractConverter implements LenientComparab
if (!isIdentity()) {
if (value instanceof DoubleDouble) {
var dd = (DoubleDouble) value;
- return dd.multiply(scale).add(offset).divide(divisor);
+ return dd.multiply(scale, true).add(offset, true).divide(divisor, true);
}
if (value instanceof BigInteger) {
value = new BigDecimal((BigInteger) value);
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 aa45ca2ef9..74f7d4ed2f 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
@@ -219,7 +219,7 @@ public final class DoubleDoubleTest extends TestCase {
@Test
@DependsOnMethod("testDivide")
public void testRatio_1m_1p() {
- final DoubleDouble t = DoubleDouble.of0(0.25).ratio_1m_1p();
+ final DoubleDouble t = DoubleDouble.of(0.25, false).ratio_1m_1p();
assertEquals((1 - 0.25) / (1 + 0.25), t.doubleValue(), STRICT);
}
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
index abc55503ab..362b05333a 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
@@ -230,12 +230,13 @@ final class GridGeometryBuilder extends GeoKeysLoader {
* scale = affine(i,i) — on the diagonal
*/
if (distance != Double.POSITIVE_INFINITY) {
+ final boolean decimal = true; // Whether values were intended to be exact in base 10.
final int numDim = affine.getNumRow() - 1;
final int trCol = affine.getNumCol() - 1;
for (int j=0; j<numDim; j++) {
final double src = -modelTiePoints.doubleValue(nearest + j);
final double tgt = modelTiePoints.doubleValue(nearest + j + Localization.RECORD_LENGTH / 2);
- final DoubleDouble t = DoubleDouble.of(src).multiply(affine.getNumber(j,j)).add(tgt);
+ var t = DoubleDouble.of(src, decimal).multiply(affine.getNumber(j,j), decimal).add(tgt, decimal);
affine.setNumber(j, trCol, t);
}
return true;