Author: desruisseaux
Date: Mon Aug 3 12:32:56 2015
New Revision: 1693890
URL: http://svn.apache.org/r1693890
Log:
Final adjustement (for now) about where to use double-double arithmetic and
where it is not worth.
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java?rev=1693890&r1=1693889&r2=1693890&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
[UTF-8] Mon Aug 3 12:32:56 2015
@@ -305,7 +305,7 @@ public abstract class MatrixSIS implemen
sum.clear();
for (int j=0; j<numRow; j++) {
get(j, i, dot);
- dot.multiply(dot);
+ dot.square();
sum.add(dot);
}
sum.sqrt();
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java?rev=1693890&r1=1693889&r2=1693890&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
[UTF-8] Mon Aug 3 12:32:56 2015
@@ -162,12 +162,12 @@ final class Initializer {
f.inverseDivide(1,0);
excentricitySquared.setFrom(f);
excentricitySquared.multiply(2,0);
- f.multiply(f);
+ f.square();
excentricitySquared.subtract(f);
} else {
final DoubleDouble rs = new DoubleDouble(b);
rs.divide(k); // rs = b/a
- rs.multiply(rs);
+ rs.square();
excentricitySquared.value = 1;
excentricitySquared.subtract(rs);
}
@@ -253,6 +253,17 @@ final class Initializer {
}
/**
+ * Returns {@code b/a} where {@code a} is the semi-major axis length and
{@code b} the semi-minor axis length.
+ * We retrieve this value from the excentricity with {@code b/a =
sqrt(1-ℯ²)}.
+ */
+ final DoubleDouble axisLengthRatio() {
+ final DoubleDouble b = new DoubleDouble(1,0);
+ b.subtract(excentricitySquared);
+ b.sqrt();
+ return b;
+ }
+
+ /**
* 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:
*
@@ -282,7 +293,7 @@ final class Initializer {
return verbatim(1 - excentricitySquared.value * (sinφ*sinφ));
}
final DoubleDouble t = verbatim(sinφ);
- t.multiply(t);
+ t.square();
t.multiply(excentricitySquared);
// Compute 1 - ℯ²⋅sin²φ. Since ℯ²⋅sin²φ may be small,
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java?rev=1693890&r1=1693889&r2=1693890&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
[UTF-8] Mon Aug 3 12:32:56 2015
@@ -21,7 +21,6 @@ import org.opengis.parameter.ParameterDe
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.OperationMethod;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
-import org.apache.sis.internal.referencing.provider.MapProjection;
import org.apache.sis.internal.referencing.provider.TransverseMercatorSouth;
import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.parameter.Parameters;
@@ -31,7 +30,6 @@ import org.apache.sis.util.Workaround;
import static java.lang.Math.*;
import static org.apache.sis.math.MathFunctions.asinh;
import static org.apache.sis.math.MathFunctions.atanh;
-import static org.apache.sis.internal.util.DoubleDouble.verbatim;
/**
@@ -119,21 +117,33 @@ public class TransverseMercator extends
super(initializer);
final double φ0 = toRadians(initializer.getAndStore(
org.apache.sis.internal.referencing.provider.TransverseMercator.LATITUDE_OF_ORIGIN));
- final double rs =
initializer.parameters.doubleValue(MapProjection.SEMI_MINOR)
- /
initializer.parameters.doubleValue(MapProjection.SEMI_MAJOR);
-
- final double n = (1 - rs) / (1 + rs); // Rewrite of n = f /
(2-f)
+ /*
+ * Opportunistically use double-double arithmetic for computation of B
since we will store
+ * it in the denormalization matrix, and there is no sine/cosine
functions involved here.
+ */
+ final double n;
+ final DoubleDouble B;
+ { // For keeping the 't' variable locale.
+ /*
+ * EPSG gives: n = f / (2-f)
+ * We rewrite as: n = (1 - b/a) / (1 + b/a)
+ */
+ final DoubleDouble t = initializer.axisLengthRatio(); // t =
b/a
+ t.ratio_1m_1p(); // t =
(1 - t) / (1 + t)
+ n = t.doubleValue();
+ /*
+ * Compute B = (1 + n²/4 + n⁴/64) / (1 + n)
+ */
+ B = new DoubleDouble(t); // B = n
+ B.square();
+ B.series(1, 0.25, 1./64); // B = (1 + n²/4 + n⁴/64)
+ t.add(1,0);
+ B.divide(t); // B = (1 + n²/4 + n⁴/64) / (1 +
n)
+ }
final double n2 = n * n;
final double n3 = n2 * n;
final double n4 = n2 * n2;
/*
- * Compute B = (n4/64 + n2/4 + 1) / (n + 1)
- * Opportunistically uses double-double arithmetic since we use it
anyway for denormalization matrix.
- */
- final DoubleDouble B = verbatim(n);
- B.add(1);
- B.inverseDivide(1, n4/64 + n2/4);
- /*
* Coefficients for direct projection.
* Add the smallest values first in order to reduce rounding errors.
*/
@@ -151,21 +161,22 @@ public class TransverseMercator extends
ih4 = (4397. / 161280)*n4;
/*
* Compute M₀ = B⋅(ξ₁ + ξ₂ + ξ₃ + ξ₄) and negate in anticipation for
what will be needed
- * in the denormalization matrix. We opportunistically use
double-double arithmetic, but
- * the precision is actually not better than double (in current SIS
version) because of
- * the precision of trigonometric functions. We may improve on that in
the future if it
- * seems useful.
+ * in the denormalization matrix. We opportunistically use
double-double arithmetic but
+ * only for the final multiplication by B, for consistency with the
translation term to
+ * be stored in the denormalization matrix. It is not worth to use
double-double in the
+ * sum of sine functions because the extra digits would be meaningless.
*
* NOTE: the EPSG documentation makes special cases for φ₀ = 0 or
±π/2. This is not
* needed here; we verified that the code below produces naturally the
expected values.
*/
final double Q = asinh(tan(φ0)) - excentricity * atanh(excentricity *
sin(φ0));
final double β = atan(sinh(Q));
- final DoubleDouble M0 = verbatim(β);
- M0.add(h1 * sin(2*β), 0);
- M0.add(h2 * sin(4*β), 0);
- M0.add(h3 * sin(6*β), 0);
- M0.add(h4 * sin(8*β), 0);
+ final DoubleDouble M0 = new DoubleDouble();
+ M0.value = h4 * sin(8*β)
+ + h3 * sin(6*β)
+ + h2 * sin(4*β)
+ + h1 * sin(2*β)
+ + β;
M0.multiply(B);
M0.negate();
/*
Modified:
sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java?rev=1693890&r1=1693889&r2=1693890&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] Mon Aug 3 12:32:56 2015
@@ -935,6 +935,24 @@ public final class DoubleDouble extends
}
/**
+ * Computes (1-x)/(1+x) where <var>x</var> is {@code this}.
+ * This pattern occurs in map projections.
+ */
+ public void ratio_1m_1p() {
+ final DoubleDouble numerator = new DoubleDouble(1, 0);
+ numerator.subtract(this);
+ add(1, 0);
+ inverseDivide(numerator);
+ }
+
+ /**
+ * Computes the square of this value.
+ */
+ public void square() {
+ multiply(value, error);
+ }
+
+ /**
* Sets this double-double value to its square root.
*
* <div class="section">Implementation</div>
@@ -968,6 +986,32 @@ public final class DoubleDouble extends
}
}
+ /**
+ * Computes c₀ + c₁x + c₂x² + c₃x³ + c₄x⁴ + … where <var>x</var> is {@code
this}.
+ * The given <var>c</var> coefficients are presumed accurate in base 2
+ * (i.e. this method does not try to apply a correction for base 10).
+ *
+ * @param coefficients The {@code c} coefficients. The array length must
be at least 1.
+ */
+ public void series(final double... coefficients) {
+ final DoubleDouble x = new DoubleDouble(this);
+ value = coefficients[0];
+ error = 0;
+ final int last = coefficients.length - 1;
+ if (last >= 1) {
+ final DoubleDouble xn = new DoubleDouble(x);
+ final DoubleDouble t = new DoubleDouble(xn);
+ for (int i=1; i<last; i++) {
+ t.multiply(coefficients[i], 0);
+ add(t);
+ xn.multiply(x);
+ t.setFrom(xn);
+ }
+ t.multiply(coefficients[last], 0);
+ add(t);
+ }
+ }
+
/**
* Returns a hash code value for this number.
*
Modified:
sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java?rev=1693890&r1=1693889&r2=1693890&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
[UTF-8] Mon Aug 3 12:32:56 2015
@@ -250,6 +250,17 @@ public final strictfp class DoubleDouble
}
/**
+ * Tests {@link DoubleDouble#ratio_1m_1p()}.
+ */
+ @Test
+ @DependsOnMethod("testDivide")
+ public void testRatio_1m_1p() {
+ final DoubleDouble t = new DoubleDouble(0.25, 0);
+ t.ratio_1m_1p();
+ assertEquals((1 - 0.25) / (1 + 0.25), t.doubleValue(), STRICT);
+ }
+
+ /**
* Tests {@link DoubleDouble#sqrt()} first with the square root of 2, then
with random values.
* In the {@code sqrt(2)} case:
*
@@ -278,7 +289,7 @@ public final strictfp class DoubleDouble
}
final double value = dd.value;
final double error = dd.error;
- dd.multiply(dd);
+ dd.square();
dd.sqrt();
dd.subtract(value, error);
assertEquals(0, dd.doubleValue(), 1E-29);
@@ -289,6 +300,17 @@ public final strictfp class DoubleDouble
}
/**
+ * Tests the {@link DoubleDouble#series(double...)} method.
+ */
+ @Test
+ @DependsOnMethod({"testMultiply", "testAdd"})
+ public void testSeries() {
+ final DoubleDouble t = new DoubleDouble(2);
+ t.series(1, 1./3, 1./9, 1./7, 1./13); // Random coefficient.
+ assertEquals(1 + 2./3 + 4./9 + 8./7 + 16./13, t.doubleValue(), STRICT);
+ }
+
+ /**
* List of all {@link DoubleDouble#VALUES} as string decimal
representations,
* together with some additional values to test.
*/
