Geogram Version 1.9.7
A programming library of geometric algorithms
Loading...
Searching...
No Matches
multi_precision.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40#ifndef GEOGRAM_NUMERICS_MULTI_PRECISION
41#define GEOGRAM_NUMERICS_MULTI_PRECISION
42
47#include <iostream>
48#include <sstream>
49#include <new>
50#include <math.h>
51
63namespace GEO {
64
65 extern double expansion_splitter_;
66 extern double expansion_epsilon_;
67
77 inline void two_sum(double a, double b, double& x, double& y) {
78 x = a + b;
79 double bvirt = x - a;
80 double avirt = x - bvirt;
81 double bround = b - bvirt;
82 double around = a - avirt;
83 y = around + bround;
84 }
85
95 inline void two_diff(double a, double b, double& x, double& y) {
96 x = a - b;
97 double bvirt = a - x;
98 double avirt = x + bvirt;
99 double bround = bvirt - b;
100 double around = a - avirt;
101 y = around + bround;
102 }
103
113 inline void fast_two_sum(double a, double b, double& x, double& y) {
114 x = a + b;
115 double bvirt = x - a;
116 y = b - bvirt;
117 }
118
128 inline void fast_two_diff(double a, double b, double& x, double& y) {
129 x = a - b;
130 double bvirt = a - x;
131 y = bvirt - b;
132 }
133
143 inline void split(double a, double& ahi, double& alo) {
144 double c = expansion_splitter_ * a;
145 double abig = c - a;
146 ahi = c - abig;
147 alo = a - ahi;
148 }
149
159 inline void two_product(double a, double b, double& x, double& y) {
160#ifdef FP_FAST_FMA
161 // If the target processor supports the FMA (Fused Multiply Add)
162 // instruction, then the product of two doubles into a length-2
163 // expansion can be implemented as follows. Thanks to Marc Glisse
164 // for the information.
165 // Note: under gcc, automatic generations of fma() for a*b+c needs
166 // to be deactivated, using -ffp-contract=off, else it may break
167 // other functions such as fast_expansion_sum_zeroelim().
168 x = a*b;
169 y = fma(a,b,-x);
170#else
171 x = a * b;
172 double ahi, alo;
173 split(a, ahi, alo);
174 double bhi, blo;
175 split(b, bhi, blo);
176 double err1 = x - (ahi * bhi);
177 double err2 = err1 - (alo * bhi);
178 double err3 = err2 - (ahi * blo);
179 y = (alo * blo) - err3;
180#endif
181 }
182
191 inline void square(double a, double& x, double& y) {
192#ifdef FP_FAST_FMA
193 // If the target processor supports the FMA (Fused Multiply Add)
194 // instruction, then the product of two doubles into a length-2
195 // expansion can be implemented as follows. Thanks to Marc Glisse
196 // for the information.
197 // Note: under gcc, automatic generations of fma() for a*b+c needs
198 // to be deactivated, using -ffp-contract=off, else it may break
199 // other functions such as fast_expansion_sum_zeroelim().
200 x = a*a;
201 y = fma(a,a,-x);
202#else
203 x = a * a;
204 double ahi, alo;
205 split(a, ahi, alo);
206 double err1 = x - (ahi * ahi);
207 double err3 = err1 - ((ahi + ahi) * alo);
208 y = (alo * alo) - err3;
209#endif
210 }
211
212 /************************************************************************/
213
227 class GEOGRAM_API expansion {
228 public:
233 index_t length() const {
234 return length_;
235 }
236
243 return capacity_;
244 }
245
251 void set_length(index_t new_length) {
252 geo_debug_assert(new_length <= capacity());
253 length_ = new_length;
254 }
255
261 const double& operator[] (index_t i) const {
262 // Note: we allocate capacity+1 storage
263 // systematically, since basic functions
264 // may access one additional value (without
265 // using it)
266 geo_debug_assert(i <= capacity_);
267 return x_[i];
268 }
269
275 double& operator[] (index_t i) {
276 // Note: we allocate capacity+1 storage
277 // systematically, since basic functions
278 // may access one additional value (without
279 // using it)
280 geo_debug_assert(i <= capacity_);
281 return x_[i];
282 }
283
289 double* data() {
290 return x_;
291 }
292
298 const double* data() const {
299 return x_;
300 }
301
309 static size_t bytes(index_t capa) {
310 // --> 2*sizeof(double) because x_ is declared of size [2]
311 // to avoid compiler's warning.
312 // --> capa+1 to have an additional 'sentry' at the end
313 // because fast_expansion_sum_zeroelim() may access
314 // an entry past the end (without using it).
315 return
316 sizeof(expansion) - 2 * sizeof(double) +
317 (capa + 1) * sizeof(double);
318 }
319
329 static size_t bytes_on_stack(index_t capa) {
330#ifndef GEO_HAS_BIG_STACK
331 // Note: standard predicates need at least 512, hence the min.
332 // index_t(MAX_CAPACITY_ON_STACK) is necessary, else with
333 // MAX_CAPACITY_ON_STACK alone the compiler tries to generate a
334 // reference to NOT_IN_LIST resulting in a link error.
335 // (weird, even with constexpr, I do not understand...)
336 // Probably when the function excepts a *reference*
338 capa <= std::max(index_t(MAX_CAPACITY_ON_STACK),index_t(512))
339 );
340#endif
341 return bytes(capa);
342 }
343
354 length_(0),
355 capacity_(capa) {
356 }
357
365#ifdef CPPCHECK
366 // cppcheck does not understand that the result
367 // of alloca() is passed to the placement syntax
368 // of operator new.
369 expansion& new_expansion_on_stack(index_t capa);
370#else
371#define new_expansion_on_stack(capa) \
372 (new (alloca(expansion::bytes_on_stack(capa)))expansion(capa))
373#endif
374
382
390
391 // ========================== Initialization from doubles
392
398 expansion& assign(double a) {
399 set_length(1);
400 x_[0] = a;
401 return *this;
402 }
403
410 geo_debug_assert(capacity() >= rhs.length());
411 set_length(rhs.length());
412 for(index_t i=0; i<rhs.length(); ++i) {
413 x_[i] = rhs.x_[i];
414 }
415 return *this;
416 }
417
424 assign(rhs);
425 if(sign() == NEGATIVE) {
426 negate();
427 }
428 return *this;
429 }
430
441 static index_t sum_capacity(double a, double b) {
442 geo_argused(a);
443 geo_argused(b);
444 return 2;
445 }
446
457 expansion& assign_sum(double a, double b) {
458 set_length(2);
459 two_sum(a, b, x_[1], x_[0]);
460 return *this;
461 }
462
473 static index_t diff_capacity(double a, double b) {
474 geo_argused(a);
475 geo_argused(b);
476 return 2;
477 }
478
489 expansion& assign_diff(double a, double b) {
490 set_length(2);
491 two_diff(a, b, x_[1], x_[0]);
492 return *this;
493 }
494
505 static index_t product_capacity(double a, double b) {
506 geo_argused(a);
507 geo_argused(b);
508 return 2;
509 }
510
521 expansion& assign_product(double a, double b) {
522 set_length(2);
523 two_product(a, b, x_[1], x_[0]);
524 return *this;
525 }
526
536 static index_t square_capacity(double a) {
537 geo_argused(a);
538 return 2;
539 }
540
551 set_length(2);
552 square(a, x_[1], x_[0]);
553 return *this;
554 }
555
556 // ====== Initialization from expansion and double
557
568 static index_t sum_capacity(const expansion& a, double b) {
569 geo_argused(b);
570 return a.length() + 1;
571 }
572
583 expansion& assign_sum(const expansion& a, double b);
584
595 static index_t diff_capacity(const expansion& a, double b) {
596 geo_argused(b);
597 return a.length() + 1;
598 }
599
610 expansion& assign_diff(const expansion& a, double b);
611
622 static index_t product_capacity(const expansion& a, double b) {
623 geo_argused(b);
624 // TODO: implement special case where the double argument
625 // is a power of two.
626 return a.length() * 2;
627 }
628
639 expansion& assign_product(const expansion& a, double b);
640
641 // ========================== Initialization from expansions
642
651 static index_t sum_capacity(const expansion& a, const expansion& b) {
652 return a.length() + b.length();
653 }
654
666
677 const expansion& a, const expansion& b, const expansion& c
678 ) {
679 return a.length() + b.length() + c.length();
680 }
681
694 const expansion& a, const expansion& b, const expansion& c
695 );
696
708 const expansion& a, const expansion& b,
709 const expansion& c, const expansion& d
710 ) {
711 return a.length() + b.length() + c.length() + d.length();
712 }
713
727 const expansion& a, const expansion& b,
728 const expansion& c, const expansion& d
729 );
730
739 static index_t diff_capacity(const expansion& a, const expansion& b) {
740 return a.length() + b.length();
741 }
742
754
764 const expansion& a, const expansion& b
765 ) {
766 return a.length() * b.length() * 2;
767 }
768
780
791 const expansion& a, const expansion& b, const expansion& c
792 ) {
793 return a.length() * b.length() * c.length() * 4;
794 }
795
808 const expansion& a, const expansion& b, const expansion& c
809 );
810
819 if(a.length() == 2) {
820 return 6;
821 } // see two_square()
822 return a.length() * a.length() * 2;
823 }
824
835
836 // ====== Determinants =============================
837
846 const expansion& a11, const expansion& a12,
847 const expansion& a21, const expansion& a22
848 ) {
849 return
850 product_capacity(a11, a22) +
851 product_capacity(a21, a12);
852 }
853
865 const expansion& a11, const expansion& a12,
866 const expansion& a21, const expansion& a22
867 );
868
878 const expansion& a11, const expansion& a12, const expansion& a13,
879 const expansion& a21, const expansion& a22, const expansion& a23,
880 const expansion& a31, const expansion& a32, const expansion& a33
881 ) {
882 // Development w.r.t. first row
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);
886 return 2 * (
887 a11.length() * c11_capa +
888 a12.length() * c12_capa +
889 a13.length() * c13_capa
890 );
891 }
892
906 const expansion& a11, const expansion& a12, const expansion& a13,
907 const expansion& a21, const expansion& a22, const expansion& a23,
908 const expansion& a31, const expansion& a32, const expansion& a33
909 );
910
921 const expansion& a21, const expansion& a22, const expansion& a23,
922 const expansion& a31, const expansion& a32, const expansion& a33
923 ) {
924 return
925 det2x2_capacity(a22, a23, a32, a33) +
926 det2x2_capacity(a23, a21, a33, a31) +
927 det2x2_capacity(a21, a22, a31, a32);
928 }
929
942 const expansion& a21, const expansion& a22, const expansion& a23,
943 const expansion& a31, const expansion& a32, const expansion& a33
944 );
945
946 // ======= Geometry-specific initializations =======
947
957 return index_t(dim) * 6;
958 }
959
973 const double* p1, const double* p2, coord_index_t dim
974 );
975
985 return index_t(dim) * 8;
986 }
987
1001 const double* p1, const double* p2, const double* p0,
1002 coord_index_t dim
1003 );
1004
1005
1014 const expansion& x, const expansion& y, const expansion& z
1015 ) {
1016 return square_capacity(x) + square_capacity(y) + square_capacity(z);
1017 }
1018
1028 const expansion& x, const expansion& y, const expansion& z
1029 );
1030
1031 // =============== some general purpose functions =========
1032
1039 static void initialize();
1040
1046 for(index_t i = 0; i < length_; ++i) {
1047 x_[i] = -x_[i];
1048 }
1049 return *this;
1050 }
1051
1060 // TODO: debug assert is_power_of_two(s)
1061 for(index_t i = 0; i < length_; ++i) {
1062 x_[i] *= s;
1063 }
1064 return *this;
1065 }
1066
1072 double estimate() const {
1073 double result = 0.0;
1074 for(index_t i = 0; i < length(); ++i) {
1075 result += x_[i];
1076 }
1077 return result;
1078 }
1079
1084 Sign sign() const {
1085 if(length() == 0) {
1086 return ZERO;
1087 }
1088 return geo_sgn(x_[length() - 1]);
1089 }
1090
1098 bool is_same_as(const expansion& rhs) const;
1099
1108 bool is_same_as(double rhs) const;
1109
1110
1115 Sign compare(const expansion& rhs) const;
1116
1121 Sign compare(double rhs) const;
1122
1129 bool equals(const expansion& rhs) const {
1130 return (compare(rhs) == ZERO);
1131 }
1132
1139 bool equals(double rhs) const {
1140 return (compare(rhs) == ZERO);
1141 }
1142
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] << " ";
1152 }
1153 out << "]";
1154 return out;
1155 }
1156
1161 std::string to_string() const {
1162 std::ostringstream out;
1163 show(out);
1164 return out.str();
1165 }
1166
1172 void optimize();
1173
1177 static void show_all_stats();
1178
1179 protected:
1190 index_t a_length, index_t b_length
1191 ) {
1192 return a_length * b_length * 2;
1193 }
1194
1205 const double* a, index_t a_length, const expansion& b
1206 );
1207
1211 expansion(const expansion& rhs) = delete;
1212
1216 expansion& operator= (const expansion& rhs) = delete;
1217
1218 private:
1219
1225#ifdef GEO_OS_APPLE
1226 static constexpr index_t MAX_CAPACITY_ON_STACK = 256;
1227#else
1228 static constexpr index_t MAX_CAPACITY_ON_STACK = 1024;
1229#endif
1230 index_t length_;
1231 index_t capacity_;
1232 double x_[2]; // x_ is in fact of size [capacity_]
1233
1234 friend class expansion_nt;
1235 };
1236
1237 // =============== arithmetic operations ===========================
1238
1251#define expansion_create(a) \
1252 new_expansion_on_stack(1)->assign(a)
1253
1254
1268#define expansion_abs(e) \
1269 new_expansion_on_stack(e.length())->assign_abs(e)
1270
1288#define expansion_sum(a, b) \
1289 new_expansion_on_stack( \
1290 expansion::sum_capacity(a, b) \
1291 )->assign_sum(a, b)
1292
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)
1314
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)
1339
1357#define expansion_diff(a, b) \
1358 new_expansion_on_stack( \
1359 expansion::diff_capacity(a, b) \
1360 )->assign_diff(a, b)
1361
1379#define expansion_product(a, b) \
1380 new_expansion_on_stack( \
1381 expansion::product_capacity(a, b) \
1382 )->assign_product(a, b)
1383
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)
1405
1421#define expansion_square(a) \
1422 new_expansion_on_stack( \
1423 expansion::square_capacity(a) \
1424 )->assign_square(a)
1425
1426 // =============== determinants =====================================
1427
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)
1445
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)
1466
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)
1488
1489 // =============== geometric functions ==============================
1490
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)
1511
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)
1534
1535
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)
1555
1556 /************************************************************************/
1557
1565 const expansion& a00,const expansion& a01,
1566 const expansion& a10,const expansion& a11
1567 );
1568
1576 const expansion& a00,const expansion& a01,const expansion& a02,
1577 const expansion& a10,const expansion& a11,const expansion& a12,
1578 const expansion& a20,const expansion& a21,const expansion& a22
1579 );
1580
1588 const expansion& a00,const expansion& a01,
1589 const expansion& a02,const expansion& a03,
1590 const expansion& a10,const expansion& a11,
1591 const expansion& a12,const expansion& a13,
1592 const expansion& a20,const expansion& a21,
1593 const expansion& a22,const expansion& a23,
1594 const expansion& a30,const expansion& a31,
1595 const expansion& a32,const expansion& a33
1596 );
1597
1598 /************************************************************************/
1599
1614 void GEOGRAM_API grow_expansion_zeroelim(
1615 const expansion& e, double b, expansion& h
1616 );
1617
1634 const expansion& e, double b, expansion& h
1635 );
1636
1653 const expansion& e, const expansion& f, expansion& h
1654 );
1655
1656
1672 const expansion& e, const expansion& f, expansion& h
1673 );
1674
1675 /************************************************************************/
1676}
1677
1678#endif
Assertion checking mechanism.
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition assert.h:196
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.
Definition basic.h:55
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.
Definition argused.h:60
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.
Definition numeric.h:106
geo_index_t index_t
The type for storing and manipulating indices.
Definition numeric.h:329
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.
Definition numeric.h:68
@ ZERO
Definition numeric.h:72
@ NEGATIVE
Definition numeric.h:70
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.
Definition numeric.h:363
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.