add in all JTS unit tests for doubledouble calculations

libgeos · Apr 16, 2020 · 8fd12e0 · 8fd12e0
1 parent e24af3b
commit 8fd12e0
Show file tree

Hide file tree

Showing 4 changed files with 483 additions and 34 deletions.
diff --git a/include/geos/math/DD.h b/include/geos/math/DD.h
@@ -15,6 +15,8 @@
 #ifndef GEOS_MATH_DD_H
 #define GEOS_MATH_DD_H
 
+#include <cmath>
+
 namespace geos {
 namespace math { // geos.math
 
@@ -27,13 +29,14 @@ namespace math { // geos.math
  */
 class GEOS_DLL DD {
     private:
-
         double SPLIT = 134217729.0; // 2^27+1, for IEEE double
         double hi;
         double lo;
 
-        // DD parse(std::string &str);
-        int magnitude(double x);
+        int magnitude(double x) const;
+        int signum() const;
+        DD rint() const;
+
 
     public:
         DD(double p_hi, double p_lo) : hi(p_hi), lo(p_lo) {};
@@ -80,24 +83,23 @@ class GEOS_DLL DD {
         friend DD operator/ (const DD &lhs, const DD &rhs);
         friend DD operator/ (const DD &lhs, double rhs);
 
+        static DD determinant(const DD &x1, const DD &y1, const DD &x2, const DD &y2);
+        static DD determinant(double x1, double y1, double x2, double y2);
+        static DD abs(const DD &d);
+        static DD pow(const DD &d, int exp);
+        static DD trunc(const DD &d);
+
         bool isNaN() const;
         bool isNegative() const;
         bool isPositive() const;
         bool isZero() const;
         double doubleValue() const;
+        double ToDouble() const { return doubleValue(); }
         int intValue() const;
-
         DD negate() const;
         DD reciprocal() const;
         DD floor() const;
         DD ceil() const;
-        int signum() const;
-        DD rint() const;
-        DD trunc() const;
-        DD abs() const;
-        DD sqr() const;
-
-        void selfSqr();
 
         void selfAdd(const DD &d);
         void selfAdd(double p_hi, double p_lo);
@@ -114,11 +116,9 @@ class GEOS_DLL DD {
         void selfDivide(double p_hi, double p_lo);
         void selfDivide(const DD &d);
         void selfDivide(double y);
-
-        double ToDouble() const { return doubleValue(); }
-
 };
 
+
 } // namespace geos::math
 } // namespace geos
 

diff --git a/src/math/DD.cpp b/src/math/DD.cpp
@@ -31,7 +31,7 @@ namespace math { // geos.util
 
 /* private */
 int
-DD::magnitude(double x)
+DD::magnitude(double x) const
 {
     double xAbs = std::fabs(x);
     double xLog10 = std::log(xAbs) / std::log(10);
@@ -262,10 +262,11 @@ DD operator/(const DD &lhs, double rhs)
     return rv;
 }
 
+/* public */
 DD DD::negate() const
 {
     DD rv(hi, lo);
-    if (isNaN())
+    if (rv.isNaN())
     {
         return rv;
     }
@@ -274,6 +275,7 @@ DD DD::negate() const
     return rv;
 }
 
+/* public static */
 DD DD::reciprocal() const
 {
     double  hc, tc, hy, ty, C, c, U, u;
@@ -286,8 +288,7 @@ DD DD::reciprocal() const
     c = ((((1.0-U)-u))-C*lo)/hi;
     double zhi = C+c;
     double zlo = (C-zhi)+c;
-    DD rv(zhi, zlo);
-    return rv;
+    return DD(zhi, zlo);
 }
 
 DD DD::floor() const
@@ -335,42 +336,75 @@ DD DD::rint() const
 {
     DD rv(hi, lo);
     if (isNaN()) return rv;
-    return (rv + 0.5).floor();
+     return (rv + 0.5).floor();
 }
 
-DD DD::trunc() const
+/* public static */
+DD DD::trunc(const DD &d)
 {
-    DD rv(hi, lo);
-    if (isNaN()) return rv;
-    if (isPositive())
+    DD rv(d);
+    if (rv.isNaN()) return rv;
+    if (rv.isPositive())
         return rv.floor();
     return rv.ceil();
 }
 
-DD DD::abs() const
+/* public static */
+DD DD::abs(const DD &d)
 {
-    DD rv(hi, lo);
-    if (isNaN()) return rv;
-    if (isNegative())
+    DD rv(d);
+    if (rv.isNaN()) return rv;
+    if (rv.isNegative())
         return rv.negate();
 
     return rv;
 }
 
-DD DD::sqr() const
+/* public static */
+DD DD::determinant(const DD &x1, const DD &y1, const DD &x2, const DD &y2)
 {
-    DD rv(hi, lo);
-    return rv * rv;
+    return (x1 * y2) - (y1 * x2);
 }
 
-void DD::selfSqr()
+/* public static */
+DD DD::determinant(double x1, double y1, double x2, double y2)
 {
-    DD rhs(hi, lo);
-    selfMultiply(rhs);
-    return;
+    return determinant(DD(x1), DD(y1), DD(x2), DD(y2) );
 }
 
+/**
+* Computes the value of this number raised to an integral power.
+* Follows semantics of Java Math.pow as closely as possible.
+*/
+/* public static */
+DD DD::pow(const DD &d, int exp)
+{
+    if (exp == 0.0)
+        return DD(1.0);
+
+    DD r(d);
+    DD s(1.0);
+    int n = std::abs(exp);
 
+    if (n > 1) {
+        /* Use binary exponentiation */
+        while (n > 0) {
+        if (n % 2 == 1) {
+            s.selfMultiply(r);
+        }
+        n /= 2;
+        if (n > 0)
+            r = r*r;
+        }
+    } else {
+        s = r;
+    }
+
+    /* Compute the reciprocal if n is negative. */
+    if (exp < 0)
+        return s.reciprocal();
+    return s;
+}
 
 
 }

diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am
@@ -151,6 +151,7 @@ geos_unit_SOURCES = \
 	io/WKTWriterTest.cpp \
 	io/WriterTest.cpp \
 	linearref/LengthIndexedLineTest.cpp \
+	math/DDTest.cpp \
 	noding/BasicSegmentStringTest.cpp \
 	noding/NodedSegmentStringTest.cpp \
 	noding/OrientedCoordinateArrayTest.cpp \