/********************************************************************** * * 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

* * * * @author Martin Davis * */ #pragma once #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 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() : 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 GEOS_DLL DD operator+ (const DD &lhs, const DD &rhs); friend GEOS_DLL DD operator+ (const DD &lhs, double rhs); friend GEOS_DLL DD operator- (const DD &lhs, const DD &rhs); friend GEOS_DLL DD operator- (const DD &lhs, double rhs); friend GEOS_DLL DD operator* (const DD &lhs, const DD &rhs); friend GEOS_DLL DD operator* (const DD &lhs, double rhs); friend GEOS_DLL DD operator/ (const DD &lhs, const DD &rhs); friend GEOS_DLL 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