DD wip

include/geos/Makefile.am

@@ -8,6 +8,7 @@ SUBDIRS = \
     index \
     io \
     linearref \
+    math \
     noding \
     operation \
     planargraph \

include/geos/math/DD.h

@@ -0,0 +1,94 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2001-2002 Vivid Solutions Inc.
+ * Copyright (C) 2020 Crunchy Data
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#ifndef GEOS_MATH_TOPOLOGYEXCEPTION_H
+#define GEOS_MATH_TOPOLOGYEXCEPTION_H
+
+#include <cassert>
+
+namespace geos {
+namespace math { // geos.math
+
+/**
+ * \class DD
+ *
+ * \brief
+ * Wrapper for DoubleDouble higher precision mathematics
+ * operations.
+ */
+class GEOS_DLL DD {
+    private:
+        double hi;
+        double lo;
+
+        DD parse(std::string &str);
+        int magnitude(double x);
+
+    public:
+        DD(double p_hi, double p_lo) : hi(p_hi), lo(p_lo) {};
+        DD(double x) : hi(x), lo(0.0) {};
+        DD(DD &dd) : hi(dd.hi), lo(dd.lo) {};
+
+        bool operator==(DD const &rhs) const
+        {
+            return hi == rhs.hi && lo == rhs.lo;
+        }
+
+        bool operator!=(DD const &rhs) const
+        {
+            return hi != rhs.hi || lo != rhs.lo;
+        }
+
+        bool operator<(DD const &rhs) const
+        {
+            return (hi < rhs.hi) || (hi == rhs.hi && lo < rhs.lo);
+        }
+
+        bool operator<=(DD const &rhs) const
+        {
+            return (hi < rhs.hi) || (hi == rhs.hi && lo <= rhs.lo);
+        }
+
+        bool operator>(DD const &rhs) const
+        {
+            return (hi > rhs.hi) || (hi == rhs.hi && lo > rhs.lo);
+        }
+
+        bool operator>=(DD const &rhs) const
+        {
+            return (hi > rhs.hi) || (hi == rhs.hi && lo >= rhs.lo);
+        }
+
+        DD operator+(DD const &rhs) const;
+        DD operator+(double rhs) const;
+
+
+        bool isNaN();
+        bool isNegative();
+        bool isPositive();
+        bool isZero();
+        double doubleValue();
+        int intValue();
+
+        void selfAdd(DD const &d);
+        void selfAdd(double p_hi, double p_lo);
+        void selfAdd(double y);
+};
+
+} // namespace geos::math
+} // namespace geos
+
+
+#endif // GEOS_MATH_TOPOLOGYEXCEPTION_H

include/geos/math/Makefile.am

@@ -0,0 +1,11 @@
+#
+# This file is part of project GEOS (http://trac.osgeo.org/geos/) 
+#
+#SUBDIRS = 
+
+#EXTRA_DIST = 
+
+geosdir = $(includedir)/geos/math
+
+geos_HEADERS = \
+    DD.h 

src/Makefile.am

@@ -7,6 +7,7 @@ SUBDIRS = \
     index \
     io \
     linearref \
+    math \
     noding \
     operation \
     planargraph \

src/math/DD.cpp

@@ -0,0 +1,136 @@
+#include <geos/profiler.h>
+#include <string>
+#include <sstream>
+#include <utility>
+
+using namespace std;
+
+#include <cmath>
+#include <geos/math/DD.h>
+
+namespace geos {
+namespace math { // geos.util
+
+
+/* private */
+DD
+DD::parse(std::string &str)
+{
+    DD d(0.0, 0.0);
+    return d;
+}
+
+/* private */
+int
+DD::magnitude(double x)
+{
+    double xAbs = std::fabs(x);
+    double xLog10 = std::log(xAbs) / std::log(10);
+    int xMag = (int) std::floor(xLog10);
+    /**
+     * Since log computation is inexact, there may be an off-by-one error
+     * in the computed magnitude.
+     * Following tests that magnitude is correct, and adjusts it if not
+     */
+    double xApprox = std::pow(10, xMag);
+    if (xApprox * 10 <= xAbs)
+      xMag += 1;
+
+    return xMag;
+}
+
+/* public */
+bool DD::isNaN()
+{
+    return std::isnan(hi);
+}
+/* public */
+bool DD::isNegative()
+{
+    return hi < 0.0 || (hi == 0.0 && lo < 0.0);
+}
+/* public */
+bool DD::isPositive()
+{
+    return hi > 0.0 || (hi == 0.0 && lo > 0.0);
+}
+/* public */
+bool DD::isZero()
+{
+    return hi == 0.0 && lo == 0.0;
+}
+
+/* public */
+double DD::doubleValue()
+{
+    return hi + lo;
+}
+
+/* public */
+int DD::intValue()
+{
+    return (int) hi;
+}
+
+/* public */
+void DD::selfAdd(DD const &y)
+{
+    return selfAdd(y.hi, y.lo);
+}
+
+/* public */
+void DD::selfAdd(double yhi, double ylo)
+{
+    double H, h, T, t, S, s, e, f;
+    S = hi + yhi;
+    T = lo + ylo;
+    e = S - hi;
+    f = T - lo;
+    s = S-e;
+    t = T-f;
+    s = (yhi-e)+(hi-s);
+    t = (ylo-f)+(lo-t);
+    e = s+T; H = S+e; h = e+(S-H); e = t+h;
+
+    double zhi = H + e;
+    double zlo = e + (H - zhi);
+    hi = zhi;
+    lo = zlo;
+    return;
+}
+
+/* public */
+void DD::selfAdd(double y)
+{
+    double H, h, S, s, e, f;
+    S = hi + y;
+    e = S - hi;
+    s = S - e;
+    s = (y - e) + (hi - s);
+    f = s + lo;
+    H = S + f;
+    h = f + (S - H);
+    hi = H + h;
+    lo = h + (H - hi);
+    return;
+}
+
+/* public */
+DD DD::operator+(DD const &rhs) const
+{
+    DD lhs(hi, lo);
+    lhs.selfAdd(rhs);
+    return lhs;
+}
+
+/* public */
+DD DD::operator+(double rhs) const
+{
+    DD lhs(hi, lo);
+    lhs.selfAdd(rhs);
+    return lhs;
+}
+
+
+}
+}

src/math/Makefile.am

@@ -0,0 +1,13 @@
+#
+# This file is part of project GEOS (http://trac.osgeo.org/geos/) 
+#
+SUBDIRS = 
+
+noinst_LTLIBRARIES = libmath.la
+
+AM_CPPFLAGS = -I$(top_srcdir)/include 
+
+libmath_la_SOURCES = \
+	DD.cpp 
+
+libmath_la_LIBADD =