Geogram Version 1.8.5
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 split(double a, double& ahi, double& alo) {
114 double c = expansion_splitter_ * a;
115 double abig = c - a;
116 ahi = c - abig;
117 alo = a - ahi;
118 }
119
129 inline void two_product(double a, double b, double& x, double& y) {
130#ifdef FP_FAST_FMA
131 // If the target processor supports the FMA (Fused Multiply Add)
132 // instruction, then the product of two doubles into a length-2
133 // expansion can be implemented as follows. Thanks to Marc Glisse
134 // for the information.
135 // Note: under gcc, automatic generations of fma() for a*b+c needs
136 // to be deactivated, using -ffp-contract=off, else it may break
137 // other functions such as fast_expansion_sum_zeroelim().
138 x = a*b;
139 y = fma(a,b,-x);
140#else
141 x = a * b;
142 double ahi, alo;
143 split(a, ahi, alo);
144 double bhi, blo;
145 split(b, bhi, blo);
146 double err1 = x - (ahi * bhi);
147 double err2 = err1 - (alo * bhi);
148 double err3 = err2 - (ahi * blo);
149 y = (alo * blo) - err3;
150#endif
151 }
152
161 inline void square(double a, double& x, double& y) {
162#ifdef FP_FAST_FMA
163 // If the target processor supports the FMA (Fused Multiply Add)
164 // instruction, then the product of two doubles into a length-2
165 // expansion can be implemented as follows. Thanks to Marc Glisse
166 // for the information.
167 // Note: under gcc, automatic generations of fma() for a*b+c needs
168 // to be deactivated, using -ffp-contract=off, else it may break
169 // other functions such as fast_expansion_sum_zeroelim().
170 x = a*a;
171 y = fma(a,a,-x);
172#else
173 x = a * a;
174 double ahi, alo;
175 split(a, ahi, alo);
176 double err1 = x - (ahi * ahi);
177 double err3 = err1 - ((ahi + ahi) * alo);
178 y = (alo * alo) - err3;
179#endif
180 }
181
182 /************************************************************************/
183
197 class GEOGRAM_API expansion {
198 public:
203 index_t length() const {
204 return length_;
205 }
206
213 return capacity_;
214 }
215
221 void set_length(index_t new_length) {
222 geo_debug_assert(new_length <= capacity());
223 length_ = new_length;
224 }
225
231 const double& operator[] (index_t i) const {
232 // Note: we allocate capacity+1 storage
233 // systematically, since basic functions
234 // may access one additional value (without
235 // using it)
236 geo_debug_assert(i <= capacity_);
237 return x_[i];
238 }
239
245 double& operator[] (index_t i) {
246 // Note: we allocate capacity+1 storage
247 // systematically, since basic functions
248 // may access one additional value (without
249 // using it)
250 geo_debug_assert(i <= capacity_);
251 return x_[i];
252 }
253
259 double* data() {
260 return x_;
261 }
262
268 const double* data() const {
269 return x_;
270 }
271
279 static size_t bytes(index_t capa) {
280 // --> 2*sizeof(double) because x_ is declared of size [2]
281 // to avoid compiler's warning.
282 // --> capa+1 to have an additional 'sentry' at the end
283 // because fast_expansion_sum_zeroelim() may access
284 // an entry past the end (without using it).
285 return
286 sizeof(expansion) - 2 * sizeof(double) +
287 (capa + 1) * sizeof(double);
288 }
289
300 length_(0),
301 capacity_(capa) {
302 }
303
311#ifdef CPPCHECK
312 // cppcheck does not understand that the result
313 // of alloca() is passed to the placement syntax
314 // of operator new.
315 expansion& new_expansion_on_stack(index_t capa);
316#else
317#define new_expansion_on_stack(capa) \
318 (new (alloca(expansion::bytes(capa)))expansion(capa))
319#endif
320
328
336
337 // ========================== Initialization from doubles
338
344 expansion& assign(double a) {
345 set_length(1);
346 x_[0] = a;
347 return *this;
348 }
349
356 geo_debug_assert(capacity() >= rhs.length());
357 set_length(rhs.length());
358 for(index_t i=0; i<rhs.length(); ++i) {
359 x_[i] = rhs.x_[i];
360 }
361 return *this;
362 }
363
370 assign(rhs);
371 if(sign() == NEGATIVE) {
372 negate();
373 }
374 return *this;
375 }
376
387 static index_t sum_capacity(double a, double b) {
388 geo_argused(a);
389 geo_argused(b);
390 return 2;
391 }
392
403 expansion& assign_sum(double a, double b) {
404 set_length(2);
405 two_sum(a, b, x_[1], x_[0]);
406 return *this;
407 }
408
419 static index_t diff_capacity(double a, double b) {
420 geo_argused(a);
421 geo_argused(b);
422 return 2;
423 }
424
435 expansion& assign_diff(double a, double b) {
436 set_length(2);
437 two_diff(a, b, x_[1], x_[0]);
438 return *this;
439 }
440
451 static index_t product_capacity(double a, double b) {
452 geo_argused(a);
453 geo_argused(b);
454 return 2;
455 }
456
467 expansion& assign_product(double a, double b) {
468 set_length(2);
469 two_product(a, b, x_[1], x_[0]);
470 return *this;
471 }
472
482 static index_t square_capacity(double a) {
483 geo_argused(a);
484 return 2;
485 }
486
497 set_length(2);
498 square(a, x_[1], x_[0]);
499 return *this;
500 }
501
502 // ====== Initialization from expansion and double
503
514 static index_t sum_capacity(const expansion& a, double b) {
515 geo_argused(b);
516 return a.length() + 1;
517 }
518
529 expansion& assign_sum(const expansion& a, double b);
530
541 static index_t diff_capacity(const expansion& a, double b) {
542 geo_argused(b);
543 return a.length() + 1;
544 }
545
556 expansion& assign_diff(const expansion& a, double b);
557
568 static index_t product_capacity(const expansion& a, double b) {
569 geo_argused(b);
570 // TODO: implement special case where the double argument
571 // is a power of two.
572 return a.length() * 2;
573 }
574
585 expansion& assign_product(const expansion& a, double b);
586
587 // ========================== Initialization from expansions
588
597 static index_t sum_capacity(const expansion& a, const expansion& b) {
598 return a.length() + b.length();
599 }
600
612
623 const expansion& a, const expansion& b, const expansion& c
624 ) {
625 return a.length() + b.length() + c.length();
626 }
627
640 const expansion& a, const expansion& b, const expansion& c
641 );
642
654 const expansion& a, const expansion& b,
655 const expansion& c, const expansion& d
656 ) {
657 return a.length() + b.length() + c.length() + d.length();
658 }
659
673 const expansion& a, const expansion& b,
674 const expansion& c, const expansion& d
675 );
676
685 static index_t diff_capacity(const expansion& a, const expansion& b) {
686 return a.length() + b.length();
687 }
688
700
710 const expansion& a, const expansion& b
711 ) {
712 return a.length() * b.length() * 2;
713 }
714
726
737 const expansion& a, const expansion& b, const expansion& c
738 ) {
739 return a.length() * b.length() * c.length() * 4;
740 }
741
754 const expansion& a, const expansion& b, const expansion& c
755 );
756
765 if(a.length() == 2) {
766 return 6;
767 } // see two_square()
768 return a.length() * a.length() * 2;
769 }
770
781
782 // ====== Determinants =============================
783
792 const expansion& a11, const expansion& a12,
793 const expansion& a21, const expansion& a22
794 ) {
795 return
796 product_capacity(a11, a22) +
797 product_capacity(a21, a12);
798 }
799
811 const expansion& a11, const expansion& a12,
812 const expansion& a21, const expansion& a22
813 );
814
824 const expansion& a11, const expansion& a12, const expansion& a13,
825 const expansion& a21, const expansion& a22, const expansion& a23,
826 const expansion& a31, const expansion& a32, const expansion& a33
827 ) {
828 // Development w.r.t. first row
829 index_t c11_capa = det2x2_capacity(a22, a23, a32, a33);
830 index_t c12_capa = det2x2_capacity(a21, a23, a31, a33);
831 index_t c13_capa = det2x2_capacity(a21, a22, a31, a32);
832 return 2 * (
833 a11.length() * c11_capa +
834 a12.length() * c12_capa +
835 a13.length() * c13_capa
836 );
837 }
838
852 const expansion& a11, const expansion& a12, const expansion& a13,
853 const expansion& a21, const expansion& a22, const expansion& a23,
854 const expansion& a31, const expansion& a32, const expansion& a33
855 );
856
867 const expansion& a21, const expansion& a22, const expansion& a23,
868 const expansion& a31, const expansion& a32, const expansion& a33
869 ) {
870 return
871 det2x2_capacity(a22, a23, a32, a33) +
872 det2x2_capacity(a23, a21, a33, a31) +
873 det2x2_capacity(a21, a22, a31, a32);
874 }
875
888 const expansion& a21, const expansion& a22, const expansion& a23,
889 const expansion& a31, const expansion& a32, const expansion& a33
890 );
891
892 // ======= Geometry-specific initializations =======
893
903 return index_t(dim) * 6;
904 }
905
919 const double* p1, const double* p2, coord_index_t dim
920 );
921
931 return index_t(dim) * 8;
932 }
933
947 const double* p1, const double* p2, const double* p0,
948 coord_index_t dim
949 );
950
951
960 const expansion& x, const expansion& y, const expansion& z
961 ) {
962 return square_capacity(x) + square_capacity(y) + square_capacity(z);
963 }
964
974 const expansion& x, const expansion& y, const expansion& z
975 );
976
977 // =============== some general purpose functions =========
978
985 static void initialize();
986
992 for(index_t i = 0; i < length_; ++i) {
993 x_[i] = -x_[i];
994 }
995 return *this;
996 }
997
1006 // TODO: debug assert is_power_of_two(s)
1007 for(index_t i = 0; i < length_; ++i) {
1008 x_[i] *= s;
1009 }
1010 return *this;
1011 }
1012
1018 double estimate() const {
1019 double result = 0.0;
1020 for(index_t i = 0; i < length(); ++i) {
1021 result += x_[i];
1022 }
1023 return result;
1024 }
1025
1030 Sign sign() const {
1031 if(length() == 0) {
1032 return ZERO;
1033 }
1034 return geo_sgn(x_[length() - 1]);
1035 }
1036
1044 bool is_same_as(const expansion& rhs) const;
1045
1054 bool is_same_as(double rhs) const;
1055
1056
1061 Sign compare(const expansion& rhs) const;
1062
1067 Sign compare(double rhs) const;
1068
1075 bool equals(const expansion& rhs) const {
1076 return (compare(rhs) == ZERO);
1077 }
1078
1085 bool equals(double rhs) const {
1086 return (compare(rhs) == ZERO);
1087 }
1088
1094 std::ostream& show(std::ostream& out) const {
1095 out << "expansion[" << length() << "] = [";
1096 for(index_t i=0; i<length(); ++i) {
1097 out << (*this)[i] << " ";
1098 }
1099 out << "]";
1100 return out;
1101 }
1102
1107 std::string to_string() const {
1108 std::ostringstream out;
1109 show(out);
1110 return out.str();
1111 }
1112
1118 void optimize();
1119
1120 protected:
1131 index_t a_length, index_t b_length
1132 ) {
1133 return a_length * b_length * 2;
1134 }
1135
1144 const double* a, index_t a_length, const expansion& b
1145 );
1146
1163#define expansion_sub_product(a, a_length, b) \
1164 new_expansion_on_stack( \
1165 sub_product_capacity(a_length, b.length()) \
1166 )->assign_sub_product(a, a_length, b)
1167
1168 private:
1172 expansion(const expansion& rhs);
1173
1177 expansion& operator= (const expansion& rhs);
1178
1179 private:
1180 index_t length_;
1181 index_t capacity_;
1182 double x_[2]; // x_ is in fact of size [capacity_]
1183
1184 friend class expansion_nt;
1185 };
1186
1187 // =============== arithmetic operations ===========================
1188
1201#define expansion_create(a) \
1202 new_expansion_on_stack(1)->assign(a)
1203
1204
1218#define expansion_abs(e) \
1219 new_expansion_on_stack(e.length())->assign_abs(e)
1220
1238#define expansion_sum(a, b) \
1239 new_expansion_on_stack( \
1240 expansion::sum_capacity(a, b) \
1241 )->assign_sum(a, b)
1242
1260#define expansion_sum3(a, b, c) \
1261 new_expansion_on_stack( \
1262 expansion::sum_capacity(a, b, c) \
1263 )->assign_sum(a, b, c)
1264
1285#define expansion_sum4(a, b, c, d) \
1286 new_expansion_on_stack( \
1287 expansion::sum_capacity(a, b, c, d) \
1288 )->assign_sum(a, b, c, d)
1289
1307#define expansion_diff(a, b) \
1308 new_expansion_on_stack( \
1309 expansion::diff_capacity(a, b) \
1310 )->assign_diff(a, b)
1311
1329#define expansion_product(a, b) \
1330 new_expansion_on_stack( \
1331 expansion::product_capacity(a, b) \
1332 )->assign_product(a, b)
1333
1351#define expansion_product3(a, b, c) \
1352 new_expansion_on_stack( \
1353 expansion::product_capacity(a, b, c) \
1354 )->assign_product(a, b, c)
1355
1371#define expansion_square(a) \
1372 new_expansion_on_stack( \
1373 expansion::square_capacity(a) \
1374 )->assign_square(a)
1375
1376 // =============== determinants =====================================
1377
1391#define expansion_det2x2(a11, a12, a21, a22) \
1392 new_expansion_on_stack( \
1393 expansion::det2x2_capacity(a11, a12, a21, a22) \
1394 )->assign_det2x2(a11, a12, a21, a22)
1395
1412#define expansion_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33) \
1413 new_expansion_on_stack( \
1414 expansion::det3x3_capacity(a11,a12,a13,a21,a22,a23,a31,a32,a33) \
1415 )->assign_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33)
1416
1434#define expansion_det_111_2x3(a21, a22, a23, a31, a32, a33) \
1435 new_expansion_on_stack( \
1436 expansion::det_111_2x3_capacity(a21, a22, a23, a31, a32, a33) \
1437 )->assign_det_111_2x3(a21, a22, a23, a31, a32, a33)
1438
1439 // =============== geometric functions ==============================
1440
1457#define expansion_sq_dist(a, b, dim) \
1458 new_expansion_on_stack( \
1459 expansion::sq_dist_capacity(dim) \
1460 )->assign_sq_dist(a, b, dim)
1461
1480#define expansion_dot_at(a, b, c, dim) \
1481 new_expansion_on_stack( \
1482 expansion::dot_at_capacity(dim) \
1483 )->assign_dot_at(a, b, c, dim)
1484
1485
1501#define expansion_length2(x,y,z) \
1502 new_expansion_on_stack( \
1503 expansion::length2_capacity(x,y,z) \
1504 )->assign_length2(x,y,z)
1505
1506 /************************************************************************/
1507
1515 const expansion& a00,const expansion& a01,
1516 const expansion& a10,const expansion& a11
1517 );
1518
1526 const expansion& a00,const expansion& a01,const expansion& a02,
1527 const expansion& a10,const expansion& a11,const expansion& a12,
1528 const expansion& a20,const expansion& a21,const expansion& a22
1529 );
1530
1538 const expansion& a00,const expansion& a01,
1539 const expansion& a02,const expansion& a03,
1540 const expansion& a10,const expansion& a11,
1541 const expansion& a12,const expansion& a13,
1542 const expansion& a20,const expansion& a21,
1543 const expansion& a22,const expansion& a23,
1544 const expansion& a30,const expansion& a31,
1545 const expansion& a32,const expansion& a33
1546 );
1547
1548 /************************************************************************/
1549
1564 void GEOGRAM_API grow_expansion_zeroelim(
1565 const expansion& e, double b, expansion& h
1566 );
1567
1584 const expansion& e, double b, expansion& h
1585 );
1586
1603 const expansion& e, const expansion& f, expansion& h
1604 );
1605
1606
1622 const expansion& e, const expansion& f, expansion& h
1623 );
1624
1625 /************************************************************************/
1626}
1627
1628#endif
1629
Assertion checking mechanism.
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition assert.h:195
Common include file, providing basic definitions. Should be included before anything else by all head...
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.
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 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.
Types and functions for memory manipulation.
Global Vorpaline namespace.
Definition algorithm.h:64
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_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:246
geo_index_t index_t
The type for storing and manipulating indices.
Definition numeric.h:287
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:224
@ ZERO
Definition numeric.h:228
@ NEGATIVE
Definition numeric.h:226
geo_coord_index_t coord_index_t
The type for storing coordinate indices, and iterating on the coordinates of a point.
Definition numeric.h:321
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.