sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1533507 - in /sis/branches/JDK7/core/sis-utility/src: main/java/org/apache/sis/internal/util/DoubleDouble.java test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
Date Fri, 18 Oct 2013 15:29:12 GMT
Author: desruisseaux
Date: Fri Oct 18 15:29:11 2013
New Revision: 1533507

URL: http://svn.apache.org/r1533507
Log:
Fix a hole in DoubleDouble.sqrt accuracy.

Modified:
    sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
    sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java

Modified: sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java?rev=1533507&r1=1533506&r2=1533507&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] (original)
+++ sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] Fri Oct 18 15:29:11 2013
@@ -698,11 +698,35 @@ public final class DoubleDouble extends 
     /**
      * Sets this double-double value to its square root.
      *
-     * @todo This method is not yet implemented with double-double precision.
+     * {@section Implementation}
+     * This method searches for a {@code (r + ε)} value where:
+     *
+     * <blockquote>(r + ε)²  =  {@linkplain #value} + {@linkplain #error}</blockquote>
+     *
+     * If we could compute {@code r = sqrt(value + error)} with enough precision, then ε
would be 0.
+     * But with the {@code double} type, we can only estimate {@code r ≈ sqrt(value)}.
However, since
+     * that <var>r</var> value should be close to the "true" value, then ε should
be small.
+     *
+     * <blockquote>value + error  =  (r + ε)²  =  r² + 2rε + ε²</blockquote>
+     *
+     * Neglecting ε² on the assumption that |ε| ≪ |r|:
+     *
+     * <blockquote>value + error  ≈  r² + 2rε</blockquote>
+     *
+     * Isolating ε:
+     *
+     * <blockquote>ε  ≈  (value + error - r²) / (2r)</blockquote>
      */
     public void sqrt() {
-        value = Math.sqrt(value);
-        error = 0;
+        final double thisValue = this.value;
+        final double thisError = this.error;
+        double r = Math.sqrt(thisValue);
+        setToProduct(r, r);
+        subtract(thisValue, thisError);
+        divide(-2*r, 0); // Multiplication by 2 does not cause any precision lost.
+        error = value;
+        value = r;
+        normalize();
     }
 
     /**

Modified: sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java?rev=1533507&r1=1533506&r2=1533507&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
[UTF-8] (original)
+++ sis/branches/JDK7/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
[UTF-8] Fri Oct 18 15:29:11 2013
@@ -81,7 +81,7 @@ public final strictfp class DoubleDouble
      * Fetches the next {@code DoubleDouble} random values and store them in the given object.
      */
     private void nextRandom(final DoubleDouble dd) {
-        dd.setToSum(nextRandom(), nextRandom());
+        dd.setToSum(nextRandom(), random.nextDouble());
     }
 
     /**
@@ -343,4 +343,40 @@ public final strictfp class DoubleDouble
         assertTrue("toRadians", DoubleDouble.errorForWellKnownValue(PI / 180)     != 0);
         assertTrue("toRadians", DoubleDouble.errorForWellKnownValue(toRadians(1)) != 0);
     }
+
+    /**
+     * Tests {@link DoubleDouble#sqrt()} first with the square root of 2, then with random
values.
+     * In the {@code sqrt(2)} case:
+     *
+     * <ul>
+     *   <li>The error using {@code double} arithmetic is approximatively 1E-16.</li>
+     *   <li>The error using double-double arithmetic is expected to be slightly less
that 1E-32.</li>
+     * </ul>
+     */
+    @Test
+    @DependsOnMethod({"testMultiply", "testDivide"})
+    public void testSqrt() {
+        final BigDecimal SQRT2 = new BigDecimal("1.414213562373095048801688724209698");
+        final DoubleDouble dd = new DoubleDouble();
+        dd.value = 2;
+        dd.sqrt();
+        assertEquals(0, SQRT2.subtract(toBigDecimal(dd)).doubleValue(), 1E-32);
+        /*
+         * If we have been able to compute √2, now test with random values.
+         * Since the range of values is approximatively [-1000 … 1000], use
+         * a tolerance value 1000 time the one that we used for √2.
+         */
+        for (int i=0; i<NUMBER_OF_REPETITIONS; i++) {
+            nextRandom(dd);
+            if (dd.value < 0) {
+                dd.negate();
+            }
+            final double value = dd.value;
+            final double error = dd.error;
+            dd.multiply(dd);
+            dd.sqrt();
+            dd.subtract(value, error);
+            assertEquals(0, dd.doubleValue(), 1E-29);
+        }
+    }
 }



Mime
View raw message