sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1693890 - in /sis/branches/JDK8/core: sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/ sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ sis-utility/src/main/java/org/apache/sis/intern...
Date Mon, 03 Aug 2015 12:32:57 GMT
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.
      */



Mime
View raw message