/**********************************************************************
*
* 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.
* <p>
* 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
* <pre>
* |x.lo| <= 0.5*ulp(x.hi)
* </pre>
* 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.
* <p>
* 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.
* <p>
* 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 <tt>strictfp</tt> modifier,
* which forces all operations to take place in the standard IEEE-754 rounding model.
* <p>
* 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.
* <p>
* 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.
* <pre>
* DD a = new DD(2.0);
* DD b = new DD(3.0);
* DD c = a.add(b);
* </pre>
* In contrast, the following approach uses only one object:
* <pre>
* DD a = new DD(2.0);
* a.selfAdd(3.0);
* </pre>
* <p>
* 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.
*
* <h3>References</h3>
* <ul>
* <li>Priest, D., <i>Algorithms for Arbitrary Precision Floating Point Arithmetic</i>,
* in P. Kornerup and D. Matula, Eds., Proc. 10th Symposium on Computer Arithmetic,
* IEEE Computer Society Press, Los Alamitos, Calif., 1991.
* <li>Yozo Hida, Xiaoye S. Li and David H. Bailey,
* <i>Quad-Double Arithmetic: Algorithms, Implementation, and Application</i>,
* manuscript, Oct 2000; Lawrence Berkeley National Laboratory Report BNL-46996.
* <li>David Bailey, <i>High Precision Software Directory</i>;
* <tt>http://crd.lbl.gov/~dhbailey/mpdist/index.html</tt>
* </ul>
*
*
* @author Martin Davis
*
*/
#pragma once
#include <cmath>
#include <geos/export.h>
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