From 1111070c0a68d0a15d1488415392658a75662be4 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Thu, 16 Apr 2020 15:31:11 -0700 Subject: [PATCH] Squashed commit of the following: commit ab1b004af900354f907f9a5d31ec514c2547ada4 Author: Paul Ramsey Date: Thu Apr 16 14:40:23 2020 -0700 remove ttmath in favour of DD commit 472c1f9a12df1a4b9628f61c93594ee162382db4 Merge: 8ccf3bf8 312c085b Author: Paul Ramsey Date: Thu Apr 16 13:24:39 2020 -0700 Merge branch 'master' of https://git.osgeo.org/gitea/geos/geos into master-dd commit 8ccf3bf874a2887b40afbaca57ee67b59b5eb40b Author: Paul Ramsey Date: Thu Apr 16 13:10:25 2020 -0700 add informational comment commit 8fd12e02f060f8e35d36a162d81d1ef94da1b784 Author: Paul Ramsey Date: Thu Apr 16 12:15:00 2020 -0700 add in all JTS unit tests for doubledouble calculations commit e24af3ba8ce7eca76840164c2edd5e066bbebe28 Author: Paul Ramsey Date: Wed Apr 15 13:47:56 2020 -0700 autotools build commit cb5942a13c5113ae42108e0422fed9a34a465abf Author: Paul Ramsey Date: Wed Apr 15 13:31:28 2020 -0700 fix doxygen complaint? commit 469037a4ca75b0f4df3638256c0fff2b61d37796 Author: Paul Ramsey Date: Wed Apr 15 12:47:14 2020 -0700 change name of ifdef guard commit 00559ec59b56a0ec6a0481b54f0ab597324d9439 Author: Paul Ramsey Date: Wed Apr 15 12:40:50 2020 -0700 allow DD to swap in for ttmath commit 75e70f7f28e751f26cc3db5b4964a86f3bae518f Author: Paul Ramsey Date: Wed Apr 15 11:28:15 2020 -0700 clean build of all DD functionality from JTS commit e42e9bc0ea3fd35bddc892fbe51b3de5677d8e57 Author: Paul Ramsey Date: Tue Apr 14 17:27:45 2020 -0700 DD wip --- .codecov.yml | 1 - Makefile.am | 1 - configure.ac | 3 +- doc/Doxyfile.in | 1 - include/geos/Makefile.am | 1 + include/geos/algorithm/CGAlgorithmsDD.h | 16 +- include/geos/algorithm/Makefile.am | 3 +- include/geos/algorithm/RayCrossingCounterDD.h | 1 - include/geos/algorithm/ttmath/COPYRIGHT | 28 - include/geos/algorithm/ttmath/Makefile.am | 24 - include/geos/algorithm/ttmath/README | 23 - include/geos/algorithm/ttmath/ttmath.h | 2880 -------- include/geos/algorithm/ttmath/ttmathbig.h | 6093 ----------------- include/geos/algorithm/ttmath/ttmathdec.h | 419 -- include/geos/algorithm/ttmath/ttmathint.h | 1923 ------ include/geos/algorithm/ttmath/ttmathmisc.h | 250 - include/geos/algorithm/ttmath/ttmathobjects.h | 812 --- include/geos/algorithm/ttmath/ttmathparser.h | 2777 -------- include/geos/algorithm/ttmath/ttmaththreads.h | 252 - include/geos/algorithm/ttmath/ttmathtypes.h | 718 -- include/geos/algorithm/ttmath/ttmathuint.h | 4189 ------------ .../geos/algorithm/ttmath/ttmathuint_noasm.h | 1038 --- .../geos/algorithm/ttmath/ttmathuint_x86.h | 1620 ----- .../geos/algorithm/ttmath/ttmathuint_x86_64.h | 1177 ---- .../ttmath/ttmathuint_x86_64_msvc.asm | 551 -- include/geos/math/DD.h | 204 + include/geos/math/Makefile.am | 11 + src/Makefile.am | 2 + src/algorithm/CGAlgorithmsDD.cpp | 7 +- src/algorithm/InteriorPointArea.cpp | 1 - src/geom/LineString.cpp | 2 +- src/math/DD.cpp | 403 ++ src/math/Makefile.am | 13 + tests/unit/Makefile.am | 1 + tests/unit/math/DDTest.cpp | 414 ++ 35 files changed, 1059 insertions(+), 24800 deletions(-) delete mode 100644 include/geos/algorithm/ttmath/COPYRIGHT delete mode 100644 include/geos/algorithm/ttmath/Makefile.am delete mode 100644 include/geos/algorithm/ttmath/README delete mode 100644 include/geos/algorithm/ttmath/ttmath.h delete mode 100644 include/geos/algorithm/ttmath/ttmathbig.h delete mode 100644 include/geos/algorithm/ttmath/ttmathdec.h delete mode 100644 include/geos/algorithm/ttmath/ttmathint.h delete mode 100644 include/geos/algorithm/ttmath/ttmathmisc.h delete mode 100644 include/geos/algorithm/ttmath/ttmathobjects.h delete mode 100644 include/geos/algorithm/ttmath/ttmathparser.h delete mode 100644 include/geos/algorithm/ttmath/ttmaththreads.h delete mode 100644 include/geos/algorithm/ttmath/ttmathtypes.h delete mode 100644 include/geos/algorithm/ttmath/ttmathuint.h delete mode 100644 include/geos/algorithm/ttmath/ttmathuint_noasm.h delete mode 100644 include/geos/algorithm/ttmath/ttmathuint_x86.h delete mode 100644 include/geos/algorithm/ttmath/ttmathuint_x86_64.h delete mode 100644 include/geos/algorithm/ttmath/ttmathuint_x86_64_msvc.asm create mode 100644 include/geos/math/DD.h create mode 100644 include/geos/math/Makefile.am create mode 100644 src/math/DD.cpp create mode 100644 src/math/Makefile.am create mode 100644 tests/unit/math/DDTest.cpp diff --git a/Makefile.am b/Makefile.am index dbc74514..f064d45d 100644 --- a/Makefile.am +++ b/Makefile.am @@ -73,5 +73,4 @@ valgrindcheck: check-local: ! find $(srcdir) -name '*.cpp' -o -name '*.h' | \ grep -v tests/xmltester/tinyxml | \ - grep -v include/geos/algorithm/ttmath | \ xargs grep -n '[[:space:]]$$' diff --git a/configure.ac b/configure.ac index 86778874..80d72e36 100644 --- a/configure.ac +++ b/configure.ac @@ -421,7 +421,6 @@ AC_OUTPUT([ include/geos/algorithm/Makefile include/geos/algorithm/locate/Makefile include/geos/algorithm/distance/Makefile - include/geos/algorithm/ttmath/Makefile include/geos/geom/Makefile include/geos/geom/prep/Makefile include/geos/geom/util/Makefile @@ -436,6 +435,7 @@ AC_OUTPUT([ include/geos/index/sweepline/Makefile include/geos/io/Makefile include/geos/linearref/Makefile + include/geos/math/Makefile include/geos/noding/Makefile include/geos/noding/snapround/Makefile include/geos/operation/Makefile @@ -468,6 +468,7 @@ AC_OUTPUT([ src/index/sweepline/Makefile src/io/Makefile src/linearref/Makefile + src/math/Makefile src/noding/Makefile src/noding/snapround/Makefile src/operation/Makefile diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index c8415081..03d9e5f7 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -853,7 +853,6 @@ EXCLUDE_PATTERNS = */examples/* \ */test/* \ */bigtest/* \ */io/markup/* \ - */ttmath/* \ config.h \ acconfig.h \ CoordinateList.cpp \ diff --git a/include/geos/Makefile.am b/include/geos/Makefile.am index cc55c4ac..04e85f8f 100644 --- a/include/geos/Makefile.am +++ b/include/geos/Makefile.am @@ -8,6 +8,7 @@ SUBDIRS = \ index \ io \ linearref \ + math \ noding \ operation \ planargraph \ diff --git a/include/geos/algorithm/CGAlgorithmsDD.h b/include/geos/algorithm/CGAlgorithmsDD.h index f596466c..583532c1 100644 --- a/include/geos/algorithm/CGAlgorithmsDD.h +++ b/include/geos/algorithm/CGAlgorithmsDD.h @@ -19,17 +19,7 @@ #ifndef GEOS_ALGORITHM_CGALGORITHMDD_H #define GEOS_ALGORITHM_CGALGORITHMDD_H #include -#include - -/// \file CGAlgorithmsDD.h - -/// \brief Close to DoubleDouble equivalent used by JTS -/// -/// Usage: `ttmath::Big` -typedef ttmath::Big DD; -//typedef ttmath::Big DD; -//typedef ttmath::Big DD; -//typedef ttmath::Big DD; +#include // Forward declarations namespace geos { @@ -39,6 +29,8 @@ class CoordinateSequence; } } +using namespace geos::math; + namespace geos { namespace algorithm { // geos::algorithm @@ -139,7 +131,7 @@ public: * The circumcentre does not necessarily lie within the triangle. For example, * the circumcentre of an obtuse isosceles triangle lies outside the triangle. * - * This method uses @ref DD extended-precision arithmetic to provide more accurate + * This method uses @ref geos::math::DD extended-precision arithmetic to provide more accurate * results than [circumcentre(Coordinate, Coordinate, Coordinate)] * (@ref geos::geom::Triangle::circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2)). * diff --git a/include/geos/algorithm/Makefile.am b/include/geos/algorithm/Makefile.am index 717e63d8..dbd5a619 100644 --- a/include/geos/algorithm/Makefile.am +++ b/include/geos/algorithm/Makefile.am @@ -3,8 +3,7 @@ # SUBDIRS = \ locate \ - distance \ - ttmath + distance EXTRA_DIST = diff --git a/include/geos/algorithm/RayCrossingCounterDD.h b/include/geos/algorithm/RayCrossingCounterDD.h index 2328ab3e..d7ba41e3 100644 --- a/include/geos/algorithm/RayCrossingCounterDD.h +++ b/include/geos/algorithm/RayCrossingCounterDD.h @@ -22,7 +22,6 @@ #include #include -#include #include diff --git a/include/geos/math/DD.h b/include/geos/math/DD.h new file mode 100644 index 00000000..7bcbb183 --- /dev/null +++ b/include/geos/math/DD.h @@ -0,0 +1,204 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * 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. + * + **********************************************************************/ + +/** + * Implements extended-precision floating-point numbers + * which maintain 106 bits (approximately 30 decimal digits) of precision. + *

+ * A DoubleDouble uses a representation containing two double-precision values. + * A number x is represented as a pair of doubles, x.hi and x.lo, + * such that the number represented by x is x.hi + x.lo, where + * 

+ *    |x.lo| <= 0.5*ulp(x.hi)
+ *
+ * and ulp(y) means "unit in the last place of y". + * The basic arithmetic operations are implemented using + * convenient properties of IEEE-754 floating-point arithmetic. + *

+ * The range of values which can be represented is the same as in IEEE-754. + * The precision of the representable numbers + * is twice as great as IEEE-754 double precision. + *

+ * The correctness of the arithmetic algorithms relies on operations + * being performed with standard IEEE-754 double precision and rounding. + * This is the Java standard arithmetic model, but for performance reasons + * Java implementations are not + * constrained to using this standard by default. + * Some processors (notably the Intel Pentium architecture) perform + * floating point operations in (non-IEEE-754-standard) extended-precision. + * A JVM implementation may choose to use the non-standard extended-precision + * as its default arithmetic mode. + * To prevent this from happening, this code uses the + * Java strictfp modifier, + * which forces all operations to take place in the standard IEEE-754 rounding model. + *

+ * The API provides both a set of value-oriented operations + * and a set of mutating operations. + * Value-oriented operations treat DoubleDouble values as + * immutable; operations on them return new objects carrying the result + * of the operation. This provides a simple and safe semantics for + * writing DoubleDouble expressions. However, there is a performance + * penalty for the object allocations required. + * The mutable interface updates object values in-place. + * It provides optimum memory performance, but requires + * care to ensure that aliasing errors are not created + * and constant values are not changed. + *

+ * For example, the following code example constructs three DD instances: + * two to hold the input values and one to hold the result of the addition. + * 

+ *     DD a = new DD(2.0);
+ *     DD b = new DD(3.0);
+ *     DD c = a.add(b);
+ *
+ * In contrast, the following approach uses only one object: + * 
+ *     DD a = new DD(2.0);
+ *     a.selfAdd(3.0);
+ *
+ *

+ * This implementation uses algorithms originally designed variously by + * Knuth, Kahan, Dekker, and Linnainmaa. + * Douglas Priest developed the first C implementation of these techniques. + * Other more recent C++ implementation are due to Keith M. Briggs and David Bailey et al. + * + *

References

+ *
    + *
  • Priest, D., Algorithms for Arbitrary Precision Floating Point Arithmetic, + * in P. Kornerup and D. Matula, Eds., Proc. 10th Symposium on Computer Arithmetic, + * IEEE Computer Society Press, Los Alamitos, Calif., 1991. + *
  • Yozo Hida, Xiaoye S. Li and David H. Bailey, + * Quad-Double Arithmetic: Algorithms, Implementation, and Application, + * manuscript, Oct 2000; Lawrence Berkeley National Laboratory Report BNL-46996. + *
  • David Bailey, High Precision Software Directory; + * http://crd.lbl.gov/~dhbailey/mpdist/index.html + *
+ * + * + * @author Martin Davis + * + */ + +#ifndef GEOS_MATH_DD_H +#define GEOS_MATH_DD_H + +#include + +namespace geos { +namespace math { // geos.math + +/** + * \class DD + * + * \brief + * Wrapper for DoubleDouble higher precision mathematics + * operations. + */ +class GEOS_DLL DD { + private: + static constexpr double SPLIT = 134217729.0; // 2^27+1, for IEEE double + double hi; + double lo; + + 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) {}; + DD(double x) : hi(x), lo(0.0) {}; + DD(const DD &dd) : hi(dd.hi), lo(dd.lo) {}; + DD() : hi(0.0), lo(0.0) {}; + + bool operator==(const DD &rhs) const + { + return hi == rhs.hi && lo == rhs.lo; + } + + bool operator!=(const DD &rhs) const + { + return hi != rhs.hi || lo != rhs.lo; + } + + bool operator<(const DD &rhs) const + { + return (hi < rhs.hi) || (hi == rhs.hi && lo < rhs.lo); + } + + bool operator<=(const DD &rhs) const + { + return (hi < rhs.hi) || (hi == rhs.hi && lo <= rhs.lo); + } + + bool operator>(const DD &rhs) const + { + return (hi > rhs.hi) || (hi == rhs.hi && lo > rhs.lo); + } + + bool operator>=(const DD &rhs) const + { + return (hi > rhs.hi) || (hi == rhs.hi && lo >= rhs.lo); + } + + friend DD operator+ (const DD &lhs, const DD &rhs); + friend DD operator+ (const DD &lhs, double rhs); + friend DD operator- (const DD &lhs, const DD &rhs); + friend DD operator- (const DD &lhs, double rhs); + friend DD operator* (const DD &lhs, const DD &rhs); + friend DD operator* (const DD &lhs, double rhs); + 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; + + void selfAdd(const DD &d); + void selfAdd(double p_hi, double p_lo); + void selfAdd(double y); + + void selfSubtract(const DD &d); + void selfSubtract(double p_hi, double p_lo); + void selfSubtract(double y); + + void selfMultiply(double p_hi, double p_lo); + void selfMultiply(const DD &d); + void selfMultiply(double y); + + void selfDivide(double p_hi, double p_lo); + void selfDivide(const DD &d); + void selfDivide(double y); +}; + + +} // namespace geos::math +} // namespace geos + + +#endif // GEOS_MATH_DD_H diff --git a/include/geos/math/Makefile.am b/include/geos/math/Makefile.am new file mode 100644 index 00000000..5396b951 --- /dev/null +++ b/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 diff --git a/src/Makefile.am b/src/Makefile.am index 1bba8e87..ea465248 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -7,6 +7,7 @@ SUBDIRS = \ index \ io \ linearref \ + math \ noding \ operation \ planargraph \ @@ -39,6 +40,7 @@ libgeos_la_LIBADD = \ index/libindex.la \ io/libio.la \ linearref/liblinearref.la \ + math/libmath.la \ noding/libnoding.la \ operation/liboperation.la \ planargraph/libplanargraph.la \ diff --git a/src/algorithm/CGAlgorithmsDD.cpp b/src/algorithm/CGAlgorithmsDD.cpp index c25c1f28..c484c1f6 100644 --- a/src/algorithm/CGAlgorithmsDD.cpp +++ b/src/algorithm/CGAlgorithmsDD.cpp @@ -35,7 +35,7 @@ namespace { double constexpr DP_SAFE_EPSILON = 1e-15; inline int -OrientationDD(DD const& dd) +OrientationDD(const DD &dd) { static DD const zero(0.0); if(dd < zero) { @@ -49,10 +49,7 @@ OrientationDD(DD const& dd) return CGAlgorithmsDD::STRAIGHT; } -// inline std::string ToStringDD(DD const& dd) -// { -// return dd.ToString(); -// } + } namespace geos { diff --git a/src/algorithm/InteriorPointArea.cpp b/src/algorithm/InteriorPointArea.cpp index 1ddddac9..61c36fd0 100644 --- a/src/algorithm/InteriorPointArea.cpp +++ b/src/algorithm/InteriorPointArea.cpp @@ -200,7 +200,6 @@ private: // edge intersects scan line, so add a crossing double xInt = intersection(p0, p1, scanY); crossings.push_back(xInt); - //checkIntersectionDD(p0, p1, scanY, xInt); } void findBestMidpoint(vector& crossings) diff --git a/src/geom/LineString.cpp b/src/geom/LineString.cpp index 65f86a0d..ac5d7fe6 100644 --- a/src/geom/LineString.cpp +++ b/src/geom/LineString.cpp @@ -45,7 +45,7 @@ using namespace geos::algorithm; namespace geos { namespace geom { // geos::geom -LineString::~LineString(){}; +LineString::~LineString(){} /*protected*/ LineString::LineString(const LineString& ls) diff --git a/src/math/DD.cpp b/src/math/DD.cpp new file mode 100644 index 00000000..673e9217 --- /dev/null +++ b/src/math/DD.cpp @@ -0,0 +1,403 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * 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. + * + **********************************************************************/ + +#include + +#include +#include + +namespace geos { +namespace math { // geos.util + + +/* private */ +int +DD::magnitude(double x) const +{ + 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() const +{ + return std::isnan(hi); +} +/* public */ +bool DD::isNegative() const +{ + return hi < 0.0 || (hi == 0.0 && lo < 0.0); +} +/* public */ +bool DD::isPositive() const +{ + return hi > 0.0 || (hi == 0.0 && lo > 0.0); +} +/* public */ +bool DD::isZero() const +{ + return hi == 0.0 && lo == 0.0; +} + +/* public */ +double DD::doubleValue() const +{ + return hi + lo; +} + +/* public */ +int DD::intValue() const +{ + return (int) hi; +} + +/* public */ +void DD::selfAdd(const DD &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 operator+(const DD &lhs, const DD &rhs) +{ + DD rv(lhs.hi, lhs.lo); + rv.selfAdd(rhs); + return rv; +} + +/* public */ +DD operator+(const DD &lhs, double rhs) +{ + DD rv(lhs.hi, lhs.lo); + rv.selfAdd(rhs); + return rv; +} + +/* public */ +void DD::selfSubtract(const DD &d) +{ + return selfAdd(-1*d.hi, -1*d.lo); +} + +/* public */ +void DD::selfSubtract(double p_hi, double p_lo) +{ + return selfAdd(-1*p_hi, -1*p_lo); +} + +/* public */ +void DD::selfSubtract(double y) +{ + return selfAdd(-1*y, 0.0); +} + +/* public */ +DD operator-(const DD &lhs, const DD &rhs) +{ + DD rv(lhs.hi, lhs.lo); + rv.selfSubtract(rhs); + return rv; +} + +/* public */ +DD operator-(const DD &lhs, double rhs) +{ + DD rv(lhs.hi, lhs.lo); + rv.selfSubtract(rhs); + return rv; +} + +/* public */ +void DD::selfMultiply(double yhi, double ylo) +{ + double hx, tx, hy, ty, C, c; + C = SPLIT * hi; hx = C-hi; c = SPLIT * yhi; + hx = C-hx; tx = hi-hx; hy = c-yhi; + C = hi*yhi; hy = c-hy; ty = yhi-hy; + c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+(hi*ylo+lo*yhi); + double zhi = C+c; hx = C-zhi; + double zlo = c+hx; + hi = zhi; + lo = zlo; + return; +} + +/* public */ +void DD::selfMultiply(DD const &d) +{ + return selfMultiply(d.hi, d.lo); +} + +/* public */ +void DD::selfMultiply(double y) +{ + return selfMultiply(y, 0.0); +} + +/* public */ +DD operator*(const DD &lhs, const DD &rhs) +{ + DD rv(lhs.hi, lhs.lo); + rv.selfMultiply(rhs); + return rv; +} + +/* public */ +DD operator*(const DD &lhs, double rhs) +{ + DD rv(lhs.hi, lhs.lo); + rv.selfMultiply(rhs); + return rv; +} + + +/* public */ +void DD::selfDivide(double yhi, double ylo) +{ + double hc, tc, hy, ty, C, c, U, u; + C = hi/yhi; c = SPLIT*C; hc =c-C; + u = SPLIT*yhi; hc = c-hc; + tc = C-hc; hy = u-yhi; U = C * yhi; + hy = u-hy; ty = yhi-hy; + u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty; + c = ((((hi-U)-u)+lo)-C*ylo)/yhi; + u = C+c; + hi = u; + lo = (C-u)+c; + return; +} + +/* public */ +void DD::selfDivide(const DD &d) +{ + return selfDivide(d.hi, d.lo); +} + +/* public */ +void DD::selfDivide(double y) +{ + return selfDivide(y, 0.0); +} + +/* public */ +DD operator/(const DD &lhs, const DD &rhs) +{ + DD rv(lhs.hi, lhs.lo); + rv.selfDivide(rhs); + return rv; +} + +/* public */ +DD operator/(const DD &lhs, double rhs) +{ + DD rv(lhs.hi, lhs.lo); + rv.selfDivide(rhs); + return rv; +} + +/* public */ +DD DD::negate() const +{ + DD rv(hi, lo); + if (rv.isNaN()) + { + return rv; + } + rv.hi = -hi; + rv.lo = -lo; + return rv; +} + +/* public static */ +DD DD::reciprocal() const +{ + double hc, tc, hy, ty, C, c, U, u; + C = 1.0/hi; + c = SPLIT*C; + hc = c-C; + u = SPLIT*hi; + hc = c-hc; tc = C-hc; hy = u-hi; U = C*hi; hy = u-hy; ty = hi-hy; + u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty; + c = ((((1.0-U)-u))-C*lo)/hi; + double zhi = C+c; + double zlo = (C-zhi)+c; + return DD(zhi, zlo); +} + +DD DD::floor() const +{ + DD rv(hi, lo); + if (isNaN()) return rv; + double fhi = std::floor(hi); + double flo = 0.0; + // Hi is already integral. Floor the low word + if (fhi == hi) { + flo = std::floor(lo); + } + // do we need to renormalize here? + rv.hi = fhi; + rv.lo = flo; + return rv; +} + +DD DD::ceil() const +{ + DD rv(hi, lo); + if (isNaN()) return rv; + double fhi = std::ceil(hi); + double flo = 0.0; + // Hi is already integral. Ceil the low word + if (fhi == hi) { + flo = std::ceil(lo); + // do we need to renormalize here? + } + rv.hi = fhi; + rv.lo = flo; + return rv; +} + +int DD::signum() const +{ + if (hi > 0) return 1; + if (hi < 0) return -1; + if (lo > 0) return 1; + if (lo < 0) return -1; + return 0; +} + +DD DD::rint() const +{ + DD rv(hi, lo); + if (isNaN()) return rv; + return (rv + 0.5).floor(); +} + +/* public static */ +DD DD::trunc(const DD &d) +{ + DD rv(d); + if (rv.isNaN()) return rv; + if (rv.isPositive()) + return rv.floor(); + return rv.ceil(); +} + +/* public static */ +DD DD::abs(const DD &d) +{ + DD rv(d); + if (rv.isNaN()) return rv; + if (rv.isNegative()) + return rv.negate(); + + return rv; +} + +/* public static */ +DD DD::determinant(const DD &x1, const DD &y1, const DD &x2, const DD &y2) +{ + return (x1 * y2) - (y1 * x2); +} + +/* public static */ +DD DD::determinant(double x1, double y1, double x2, double y2) +{ + 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/src/math/Makefile.am b/src/math/Makefile.am new file mode 100644 index 00000000..c30e97ea --- /dev/null +++ b/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 = diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am index 88c4816f..3c836fc8 100644 --- 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 \ diff --git a/tests/unit/math/DDTest.cpp b/tests/unit/math/DDTest.cpp new file mode 100644 index 00000000..00a29b8f --- /dev/null +++ b/tests/unit/math/DDTest.cpp @@ -0,0 +1,414 @@ +// +// Test Suite for geos::util::UniqueCoordinateArrayFilter class. + +// geos + +#include +#include + +// tut +#include +#include + +// std +#include +#include + +using namespace geos::math; + +namespace tut { +// +// Test Group +// + +// Common data used in test cases. +struct test_dd_data { + + double eps; + DD pi; + DD e; + + + void ensure_dd_equals(const char *str, const DD &d1, const DD &d2, double tolerance) + { + DD delta = DD::abs(d1 - d2); + double diff = delta.doubleValue(); + ensure(str, diff <= tolerance); + } + + void checkTrunc(const DD &x, const DD &expected) + { + DD trunc = DD::trunc(x); + ensure("checkTrunc", trunc == expected); + } + + void checkDeterminant(double x1, double y1, double x2, double y2, double expected, double errBound) + { + DD det = DD::determinant(x1, y1, x2, y2); + //ensure_equals("1", Angle::angle(Coordinate(10, 0)), 0.0, TOL); + ensure_dd_equals("checkDeterminant", det, DD(expected), errBound); + } + + void checkDeterminantDD(double x1, double y1, double x2, double y2, double expected, double errBound) + { + DD det = DD::determinant(DD(x1), DD(y1), DD(x2), DD(y2)); + ensure_dd_equals("checkDeterminantDD", det, DD(expected), errBound); + } + + void checkAddMult2(const DD &dd) + { + DD sum = dd + dd; + DD prod = dd * DD(2.0); + ensure_dd_equals("checkAddMult2", sum, prod, 0.0); + } + + void checkMultiplyDivide(const DD &a, const DD &b, double errBound) + { + DD a2 = (a * b) / b; + ensure_dd_equals("checkMultiplyDivide", a, a2, errBound); + } + + void checkDivideMultiply(const DD &a, const DD &b, double errBound) + { + DD a2 = (a / b) * b; + ensure_dd_equals("checkDivideMultiply", a, a2, errBound); + } + + /** + * Computes (a+b)^2 in two different ways and compares the result. + * For correct results, a and b should be integers. + */ + void checkBinomialSquare(double a, double b) + { + // binomial square + DD add(a); + DD bdd(b); + DD aPlusb = add + bdd; + DD abSq = aPlusb * aPlusb; + + // expansion + DD a2dd = add * add; + DD b2dd = bdd * bdd; + DD ab = add * bdd; + DD sum = b2dd + ab + ab; + DD diff = abSq - a2dd; + DD delta = diff - sum; + + ensure("isSame", diff == sum); + ensure("isDeltaZero", delta.isZero()); + } + + void checkBinomial2(double a, double b) + { + // binomial product + DD add(a); + DD bdd(b); + DD aPlusb = add + bdd; + DD aSubb = add - bdd; + DD abProd = aPlusb * aSubb; + + // expansion + DD a2dd = add * add; + DD b2dd = bdd * bdd; + + // this should equal b^2 + DD diff = (abProd - a2dd).negate(); + DD delta = diff - b2dd; + + ensure("isSame", diff == b2dd); + ensure("isDeltaZero", delta.isZero()); + } + + void checkReciprocal(double x, double errBound) + { + DD xdd(x); + DD rr = xdd.reciprocal().reciprocal(); + double err = (xdd - rr).doubleValue(); + ensure("checkReciprocal", err <= errBound); + } + + DD slowPow(const DD &x, int exp) + { + if (exp == 0) + return DD(1.0); + + int n = std::abs(exp); + // MD - could use binary exponentiation for better precision & speed + DD pow(x); + for (int i = 1; i < n; i++) { + pow = pow * x; + } + if (exp < 0) { + return pow.reciprocal(); + } + return pow; + } + + void checkPow(double x, int exp, double errBound) + { + DD xdd(x); + DD pow = DD::pow(xdd, exp); + DD pow2 = slowPow(xdd, exp); + double err = (pow - pow2).doubleValue(); + ensure("checkPow", err <= errBound); + } + + + DD arctan(DD x) + { + DD t = x; + DD t2 = t*t; + DD at(0.0); + DD two(2.0); + int k = 0; + DD d(1.0); + int sign = 1; + while (t.doubleValue() > eps) { + k++; + if (sign < 0) + at = at - (t / d); + else + at = at + (t / d); + + d = d + two; + t = t * t2; + sign = -sign; + } + return at; + } + + /** + * Uses Taylor series to compute e + * + * e = 1 + 1 + 1/2! + 1/3! + 1/4! + ... + */ + DD computeEByTaylorSeries() + { + DD s(2.0); + DD t(1.0); + double n = 1.0; + int i = 0; + while(t.doubleValue() > eps) + { + i++; + n += 1.0; + t = t / DD(n); + s = s + t; + } + return s; + } + + /** + * Uses Machin's arctangent formula to compute Pi: + * + * Pi / 4 = 4 arctan(1/5) - arctan(1/239) + */ + + DD computePiByMachin() + { + DD t1 = DD(1.0) / DD(5.0); + DD t2 = DD(1.0) / DD(239.0); + DD pi4 = (DD(4.0) * arctan(t1)) - arctan(t2); + DD pi = DD(4.0) * pi4; + return pi; + } + + test_dd_data(): + eps(1.23259516440783e-32), /* = 2^-106 */ + pi(DD(3.141592653589793116e+00, 1.224646799147353207e-16)), + e(DD(2.718281828459045091e+00, 1.445646891729250158e-16)) + {} + +}; + +typedef test_group group; +typedef group::object object; + +group test_dd_group("geos::math::DD"); + +// +// Test Cases +// + +// Test PI calculation +template<> +template<> +void object::test<1> +() +{ + DD testPi = computePiByMachin(); + double err = std::abs((testPi - pi).doubleValue()); + // std::cout << "Difference from PI = " << err << std::endl; + ensure("Test PI calculation", err < 8*eps); +} + +// Test E calculation +template<> +template<> +void object::test<2> +() +{ + DD testE = computeEByTaylorSeries(); + double err = std::abs((testE - e).doubleValue()); + // std::cout << "Difference from E = " << err << std::endl; + ensure("Test E calculation", err < eps); +} + + +// Test NaN +template<> +template<> +void object::test<3> +() +{ + DD nan = DD(1.0) / DD(0.0); + ensure("isNan", nan.isNaN()); + ensure("isNan", (DD(1.0) * nan).isNaN()); +} + +// testAddMult2 +template<> +template<> +void object::test<4> +() +{ + checkAddMult2(DD(3.0)); + checkAddMult2(DD(pi)); +} + + +// testMultiplyDivide +template<> +template<> +void object::test<5> +() +{ + checkMultiplyDivide(DD(pi), DD(e), 1e-30); + checkMultiplyDivide(DD(DD(2.0)*pi), DD(e), 1e-30); + checkMultiplyDivide(DD(DD(0.5)*pi), DD(e), 1e-30); + checkMultiplyDivide(DD(39.4), DD(10), 1e-30); +} + + +// testDivideMultiply +template<> +template<> +void object::test<6> +() +{ + checkDivideMultiply(DD(pi), DD(e), 1e-30); + checkDivideMultiply(DD(39.4), DD(10), 1e-30); +} + + +// testTrunc +template<> +template<> +void object::test<7> +() +{ + checkTrunc(DD(1e16) - DD(1), DD(1e16) - DD(1)); + // the appropriate error bound is determined empirically + checkTrunc(DD(pi), DD(3)); + checkTrunc(DD(999.999), DD(999)); + + checkTrunc(DD(e).negate(), DD(-2)); + checkTrunc(DD(-999.999), DD(-999)); +} + +// testPow +template<> +template<> +void object::test<8> +() +{ + checkPow(0, 3, 16 * eps); + checkPow(0, 3, 16 * eps); + checkPow(14, 3, 16 * eps); + checkPow(3, -5, 16 * eps); + checkPow(-3, 5, 16 * eps); + checkPow(-3, -5, 16 * eps); + checkPow(0.12345, -5, 1e5 * eps); +} + + +// testReciprocal +template<> +template<> +void object::test<9> +() +{ + checkReciprocal(3.0, 0); + checkReciprocal(99.0, 1e-29); + checkReciprocal(999.0, 0); + checkReciprocal(314159269.0, 0); +} + +// testDeterminant +template<> +template<> +void object::test<10> +() +{ + checkDeterminant(3, 8, 4, 6, -14, 0); + checkDeterminantDD(3, 8, 4, 6, -14, 0); +} + +// testDeterminantRobust +template<> +template<> +void object::test<11> +() +{ + checkDeterminant(1.0e9, 1.0e9 - 1, 1.0e9 - 1, 1.0e9 - 2, -1, 0); + checkDeterminantDD(1.0e9, 1.0e9 - 1, 1.0e9 - 1, 1.0e9 - 2, -1, 0); +} + + +// testBinom +template<> +template<> +void object::test<12> +() +{ + checkBinomialSquare(100.0, 1.0); + checkBinomialSquare(1000.0, 1.0); + checkBinomialSquare(10000.0, 1.0); + checkBinomialSquare(100000.0, 1.0); + checkBinomialSquare(1000000.0, 1.0); + checkBinomialSquare(1e8, 1.0); + checkBinomialSquare(1e10, 1.0); + checkBinomialSquare(1e14, 1.0); + // Following call will fail, because it requires 32 digits of precision + // checkBinomialSquare(1e16, 1.0); + + checkBinomialSquare(1e14, 291.0); + checkBinomialSquare(5e14, 291.0); + checkBinomialSquare(5e14, 345291.0); +} + +// testBinom2 +template<> +template<> +void object::test<13> +() +{ + checkBinomial2(100.0, 1.0); + checkBinomial2(1000.0, 1.0); + checkBinomial2(10000.0, 1.0); + checkBinomial2(100000.0, 1.0); + checkBinomial2(1000000.0, 1.0); + checkBinomial2(1e8, 1.0); + checkBinomial2(1e10, 1.0); + checkBinomial2(1e14, 1.0); + + checkBinomial2(1e14, 291.0); + + checkBinomial2(5e14, 291.0); + checkBinomial2(5e14, 345291.0); +} + + + +} // namespace tut + -- 2.26.2