40#ifndef GEOGRAM_NUMERICS_MULTI_PRECISION
41#define GEOGRAM_NUMERICS_MULTI_PRECISION
65 extern double expansion_splitter_;
66 extern double expansion_epsilon_;
77 inline void two_sum(
double a,
double b,
double& x,
double& y) {
80 double avirt = x - bvirt;
81 double bround = b - bvirt;
82 double around = a - avirt;
95 inline void two_diff(
double a,
double b,
double& x,
double& y) {
98 double avirt = x + bvirt;
99 double bround = bvirt - b;
100 double around = a - avirt;
115 double bvirt = x - a;
130 double bvirt = a - x;
143 inline void split(
double a,
double& ahi,
double& alo) {
144 double c = expansion_splitter_ * a;
159 inline void two_product(
double a,
double b,
double& x,
double& y) {
176 double err1 = x - (ahi * bhi);
177 double err2 = err1 - (alo * bhi);
178 double err3 = err2 - (ahi * blo);
179 y = (alo * blo) - err3;
191 inline void square(
double a,
double& x,
double& y) {
206 double err1 = x - (ahi * ahi);
207 double err3 = err1 - ((ahi + ahi) * alo);
208 y = (alo * alo) - err3;
253 length_ = new_length;
317 (capa + 1) *
sizeof(
double);
330#ifndef GEO_HAS_BIG_STACK
371#define new_expansion_on_stack(capa) \
372 (new (alloca(expansion::bytes_on_stack(capa)))expansion(capa))
459 two_sum(a, b, x_[1], x_[0]);
491 two_diff(a, b, x_[1], x_[0]);
523 two_product(a, b, x_[1], x_[0]);
552 square(a, x_[1], x_[0]);
850 product_capacity(a11, a22) +
851 product_capacity(a21, a12);
883 index_t c11_capa = det2x2_capacity(a22, a23, a32, a33);
884 index_t c12_capa = det2x2_capacity(a21, a23, a31, a33);
885 index_t c13_capa = det2x2_capacity(a21, a22, a31, a32);
925 det2x2_capacity(a22, a23, a32, a33) +
926 det2x2_capacity(a23, a21, a33, a31) +
927 det2x2_capacity(a21, a22, a31, a32);
1001 const double* p1,
const double* p2,
const double* p0,
1016 return square_capacity(x) + square_capacity(y) + square_capacity(z);
1046 for(
index_t i = 0; i < length_; ++i) {
1061 for(
index_t i = 0; i < length_; ++i) {
1073 double result = 0.0;
1074 for(
index_t i = 0; i < length(); ++i) {
1088 return geo_sgn(x_[length() - 1]);
1130 return (compare(rhs) ==
ZERO);
1140 return (compare(rhs) ==
ZERO);
1148 std::ostream&
show(std::ostream& out)
const {
1149 out <<
"expansion[" << length() <<
"] = [";
1150 for(
index_t i=0; i<length(); ++i) {
1151 out << (*this)[i] <<
" ";
1162 std::ostringstream out;
1192 return a_length * b_length * 2;
1226 static constexpr index_t MAX_CAPACITY_ON_STACK = 256;
1228 static constexpr index_t MAX_CAPACITY_ON_STACK = 1024;
1251#define expansion_create(a) \
1252 new_expansion_on_stack(1)->assign(a)
1268#define expansion_abs(e) \
1269 new_expansion_on_stack(e.length())->assign_abs(e)
1288#define expansion_sum(a, b) \
1289 new_expansion_on_stack( \
1290 expansion::sum_capacity(a, b) \
1310#define expansion_sum3(a, b, c) \
1311 new_expansion_on_stack( \
1312 expansion::sum_capacity(a, b, c) \
1313 )->assign_sum(a, b, c)
1335#define expansion_sum4(a, b, c, d) \
1336 new_expansion_on_stack( \
1337 expansion::sum_capacity(a, b, c, d) \
1338 )->assign_sum(a, b, c, d)
1357#define expansion_diff(a, b) \
1358 new_expansion_on_stack( \
1359 expansion::diff_capacity(a, b) \
1360 )->assign_diff(a, b)
1379#define expansion_product(a, b) \
1380 new_expansion_on_stack( \
1381 expansion::product_capacity(a, b) \
1382 )->assign_product(a, b)
1401#define expansion_product3(a, b, c) \
1402 new_expansion_on_stack( \
1403 expansion::product_capacity(a, b, c) \
1404 )->assign_product(a, b, c)
1421#define expansion_square(a) \
1422 new_expansion_on_stack( \
1423 expansion::square_capacity(a) \
1441#define expansion_det2x2(a11, a12, a21, a22) \
1442 new_expansion_on_stack( \
1443 expansion::det2x2_capacity(a11, a12, a21, a22) \
1444 )->assign_det2x2(a11, a12, a21, a22)
1462#define expansion_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33) \
1463 new_expansion_on_stack( \
1464 expansion::det3x3_capacity(a11,a12,a13,a21,a22,a23,a31,a32,a33) \
1465 )->assign_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33)
1484#define expansion_det_111_2x3(a21, a22, a23, a31, a32, a33) \
1485 new_expansion_on_stack( \
1486 expansion::det_111_2x3_capacity(a21, a22, a23, a31, a32, a33) \
1487 )->assign_det_111_2x3(a21, a22, a23, a31, a32, a33)
1507#define expansion_sq_dist(a, b, dim) \
1508 new_expansion_on_stack( \
1509 expansion::sq_dist_capacity(dim) \
1510 )->assign_sq_dist(a, b, dim)
1530#define expansion_dot_at(a, b, c, dim) \
1531 new_expansion_on_stack( \
1532 expansion::dot_at_capacity(dim) \
1533 )->assign_dot_at(a, b, c, dim)
1551#define expansion_length2(x,y,z) \
1552 new_expansion_on_stack( \
1553 expansion::length2_capacity(x,y,z) \
1554 )->assign_length2(x,y,z)
Assertion checking mechanism.
#define geo_debug_assert(x)
Verifies that a condition is met.
Expansion_nt (expansion Number Type) is used to compute the sign of polynoms exactly.
Represents numbers in arbitrary precision with a low-level API.
static index_t length2_capacity(const expansion &x, const expansion &y, const expansion &z)
Computes the required capacity to store the length of a 3d vector.
static index_t product_capacity(const expansion &a, const expansion &b)
Computes the required capacity of an expansion to store the exact product of two expansions.
void two_product(double a, double b, double &x, double &y)
Multiplies two doubles into a length 2 expansion.
expansion & assign_sum(const expansion &a, const expansion &b, const expansion &c, const expansion &d)
Assigns the sum of four expansions to this expansion (should not be used by client code).
static index_t product_capacity(double a, double b)
Computes the required capacity of an expansion to store the exact product of two doubles.
void two_diff(double a, double b, double &x, double &y)
Subtracts two doubles into a length 2 expansion.
void split(double a, double &ahi, double &alo)
Splits a number into two components, ready for computing a product.
expansion & assign(const expansion &rhs)
Copies an expansion to this expansion.
static index_t sum_capacity(double a, double b)
Computes the required capacity to store the sum of two doubles.
static index_t sum_capacity(const expansion &a, const expansion &b, const expansion &c, const expansion &d)
Computes the required capacity of an expansion to store the exact sum of four expansions.
index_t length() const
Gets the length of this expansion.
expansion & assign_diff(const expansion &a, const expansion &b)
Assigns the difference between two expansions to this expansion (should not be used by client code).
expansion & assign_product(const expansion &a, const expansion &b, const expansion &c)
Assigns the product of three expansions to this expansion (should not be used by client code).
static index_t sub_product_capacity(index_t a_length, index_t b_length)
Computes the required capacity of an expansion to store an exact sub-product.
expansion & scale_fast(double s)
Multiplies this expansion by a power of two.
static void delete_expansion_on_heap(expansion *e)
Deallocates an expansion on the heap.
static index_t sum_capacity(const expansion &a, const expansion &b)
Computes the required capacity of an expansion to store the exact sum of two expansions.
static void initialize()
Initializes the expansion class.
static index_t sum_capacity(const expansion &a, const expansion &b, const expansion &c)
Computes the required capacity of an expansion to store the exact sum of three expansions.
void square(double a, double &x, double &y)
Squares a number into a length 2 expansion.
static index_t diff_capacity(double a, double b)
Computes the required capacity of an expansion to store the exact difference of two doubles.
bool equals(const expansion &rhs) const
Compares two expansions.
Sign compare(const expansion &rhs) const
Compares two expansions.
expansion & assign_det2x2(const expansion &a11, const expansion &a12, const expansion &a21, const expansion &a22)
Assigns a 2x2 determinant to this expansion (should not be used by client code).
expansion & assign_det_111_2x3(const expansion &a21, const expansion &a22, const expansion &a23, const expansion &a31, const expansion &a32, const expansion &a33)
Assigns a 3x3 determinant to this expansion where the first row is 1 1 1(should not be used by client...
expansion & assign_product(const expansion &a, double b)
Assigns the product between an expansion and a double to this expansion (should not be used by client...
index_t capacity() const
Gets the capacity of this expansion.
std::ostream & show(std::ostream &out) const
Displays all the components of this expansion (for debugging purposes).
expansion & assign_abs(const expansion &rhs)
Copies the absolute value of an expansion to this expansion.
expansion & assign_dot_at(const double *p1, const double *p2, const double *p0, coord_index_t dim)
Assigns the dot product of two vectors to this expansion (should not be used by client code).
static index_t sum_capacity(const expansion &a, double b)
Computes the required capacity of an expansion to store the exact sum of an expansion and a double.
double * data()
Low level access to the array of components.
expansion & assign_product(const expansion &a, const expansion &b)
Assigns the product of two expansions to this expansion (should not be used by client code).
static index_t square_capacity(const expansion &a)
Computes the required capacity of an expansion to store the exact square of an expansion.
expansion & assign_sub_product(const double *a, index_t a_length, const expansion &b)
Assigns a sub-product to this expansion.
const double * data() const
Low level access to the array of components.
expansion & assign_diff(double a, double b)
Assigns the difference of two doubles to this expansion (should not be used by client code).
expansion & assign_square(const expansion &a)
Assigns the product of an expansions to this expansion (should not be used by client code).
void set_length(index_t new_length)
Changes the length of an expansion.
expansion & assign_det3x3(const expansion &a11, const expansion &a12, const expansion &a13, const expansion &a21, const expansion &a22, const expansion &a23, const expansion &a31, const expansion &a32, const expansion &a33)
Assigns a 3x3 determinant to this expansion (should not be used by client code).
bool is_same_as(double rhs) const
Compares an expansion and a double bit-by-bit.
expansion & assign(double a)
Assigns a number to this expansion.
expansion & assign_diff(const expansion &a, double b)
Assigns the difference between an expansion and a double to this expansion (should not be used by cli...
expansion & assign_sum(const expansion &a, double b)
Assigns the sum of an expansion and a double to this expansion (should not be used by client code).
static index_t det_111_2x3_capacity(const expansion &a21, const expansion &a22, const expansion &a23, const expansion &a31, const expansion &a32, const expansion &a33)
Computes the required capacity of an expansion to store an exact 3x3 determinant where the first row ...
expansion & assign_sum(double a, double b)
Assigns the sum of two doubles to this expansion (should not be used by client code).
bool equals(double rhs) const
Compares an expansion and a double.
static index_t det2x2_capacity(const expansion &a11, const expansion &a12, const expansion &a21, const expansion &a22)
Computes the required capacity of an expansion to store an exact 2x2 determinant.
expansion(index_t capa)
Client code should not use this constructor.
static index_t diff_capacity(const expansion &a, const expansion &b)
Computes the required capacity of an expansion to store the exact difference of two expansions.
expansion(const expansion &rhs)=delete
Expansions cannot be copied.
static size_t bytes_on_stack(index_t capa)
Computes the amount of memory required to store an expansion on the stack.
static index_t product_capacity(const expansion &a, double b)
Computes the required capacity of an expansion to store the exact product between an expansion and a ...
static index_t sq_dist_capacity(coord_index_t dim)
Computes the required capacity of an expansion to store the exact squared distance between two points...
expansion & negate()
Changes the sign of an expansion.
static index_t product_capacity(const expansion &a, const expansion &b, const expansion &c)
Computes the required capacity of an expansion to store the exact product of three expansions.
static size_t bytes(index_t capa)
Computes the amount of memory required to store an expansion.
static index_t det3x3_capacity(const expansion &a11, const expansion &a12, const expansion &a13, const expansion &a21, const expansion &a22, const expansion &a23, const expansion &a31, const expansion &a32, const expansion &a33)
Computes the required capacity of an expansion to store an exact 3x3 determinant.
double estimate() const
Computes an approximation of the stored value in this expansion.
static index_t square_capacity(double a)
Computes the required capacity of an expansion to store the exact square of a double.
expansion & assign_sum(const expansion &a, const expansion &b)
Assigns the sum of two expansions to this expansion (should not be used by client code).
static void show_all_stats()
Show global statistics.
static index_t dot_at_capacity(coord_index_t dim)
Computes the required capacity of an expansion to store the exact dot product between two vectors.
bool is_same_as(const expansion &rhs) const
Compares two expansions bit-by-bit.
expansion & assign_sum(const expansion &a, const expansion &b, const expansion &c)
Assigns the sum of three expansions to this expansion (should not be used by client code).
static index_t diff_capacity(const expansion &a, double b)
Computes the required capacity of an expansion to store the exact difference between an expansion and...
std::string to_string() const
Gets a string representation of this expansion.
static expansion * new_expansion_on_heap(index_t capa)
Allocates an expansion on the heap.
expansion & assign_length2(const expansion &x, const expansion &y, const expansion &z)
Assigns the length of a vector to this expansion (should not be used by client code)....
expansion & assign_square(double a)
Assigns the square of a double to this expansion (should not be used by client code).
expansion & assign_sq_dist(const double *p1, const double *p2, coord_index_t dim)
Assigns the squared distance between two points to this expansion (should not be used by client code)...
void two_sum(double a, double b, double &x, double &y)
Sums two doubles into a length 2 expansion.
void optimize()
Optimizes the internal representation without changing the represented value.
expansion & assign_product(double a, double b)
Assigns the product of two doubles to this expansion (should not be used by client code).
Sign compare(double rhs) const
Compares two expansions.
Sign sign() const
Gets the sign of the expansion.
Common include file, providing basic definitions. Should be included before anything else by all head...
Types and functions for memory manipulation.
Global Vorpaline namespace.
void scale_expansion_zeroelim(const expansion &e, double b, expansion &h)
Multiplies an expansion by a scalar, eliminating zero components from the output expansion.
void grow_expansion_zeroelim(const expansion &e, double b, expansion &h)
Adds a scalar to an expansion, eliminating zero components from the output expansion.
void geo_argused(const T &)
Suppresses compiler warnings about unused parameters.
void fast_two_diff(double a, double b, double &x, double &y)
Computes the difference of two doubles into a length 2 expansion.
void fast_expansion_sum_zeroelim(const expansion &e, const expansion &f, expansion &h)
Sums two expansions, eliminating zero components from the output expansion (sets h = e + f).
Sign geo_sgn(const T &x)
Gets the sign of a value.
geo_index_t index_t
The type for storing and manipulating indices.
void fast_expansion_diff_zeroelim(const expansion &e, const expansion &f, expansion &h)
Computes the difference of two expansions, eliminating zero components from the output expansion.
Sign
Integer constants that represent the sign of a value.
void fast_two_sum(double a, double b, double &x, double &y)
Computes the sum of two doubles into a length 2 expansion.
geo_coord_index_t coord_index_t
The type for storing coordinate indices, and iterating on the coordinates of a point.
Sign sign_of_expansion_determinant(const expansion &a00, const expansion &a01, const expansion &a10, const expansion &a11)
Computes the sign of a 2x2 determinant.
Types and functions for numbers manipulation.