| Line | Branch | Exec | Source |
|---|---|---|---|
| 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 | |||
| 43 | #include <geogram/basic/common.h> | ||
| 44 | #include <geogram/basic/numeric.h> | ||
| 45 | #include <geogram/basic/memory.h> | ||
| 46 | #include <geogram/basic/assert.h> | ||
| 47 | #include <iostream> | ||
| 48 | #include <sstream> | ||
| 49 | #include <new> | ||
| 50 | #include <math.h> | ||
| 51 | |||
| 52 | /** | ||
| 53 | * \file geogram/numerics/multi_precision.h | ||
| 54 | * \brief Implementation of multi-precision arithmetics | ||
| 55 | * \details | ||
| 56 | * Multi-precision arithmetics based on expansions, as described by | ||
| 57 | * Jonathan Shewchuk in: | ||
| 58 | * Adaptive Precision Floating-Point Arithmetic and Fast Robust | ||
| 59 | * Geometric Predicates, | ||
| 60 | * Discrete & Computational Geometry 18(3):305-363, October 1997 | ||
| 61 | */ | ||
| 62 | |||
| 63 | namespace GEO { | ||
| 64 | |||
| 65 | extern bool expansion_initialized_; | ||
| 66 | extern double expansion_splitter_; | ||
| 67 | extern double expansion_epsilon_; | ||
| 68 | |||
| 69 | /** | ||
| 70 | * \brief Sums two doubles into a length 2 expansion. | ||
| 71 | * \details By Jonathan Shewchuk. | ||
| 72 | * \param[in] a one of the numbers to sum. | ||
| 73 | * \param[in] b the other numbers to sum. | ||
| 74 | * \param[out] x high-magnitude component of the result. | ||
| 75 | * \param[out] y low-magnitude component of the result. | ||
| 76 | * \relates expansion | ||
| 77 | */ | ||
| 78 | 2017427205 | inline void two_sum(double a, double b, double& x, double& y) { | |
| 79 | 2017427205 | x = a + b; | |
| 80 | 2017427205 | double bvirt = x - a; | |
| 81 | 2017427205 | double avirt = x - bvirt; | |
| 82 | 2017427205 | double bround = b - bvirt; | |
| 83 | 2017427205 | double around = a - avirt; | |
| 84 | 2017427205 | y = around + bround; | |
| 85 | 2017427205 | } | |
| 86 | |||
| 87 | /** | ||
| 88 | * \brief Subtracts two doubles into a length 2 expansion. | ||
| 89 | * \details By Jonathan Shewchuk. | ||
| 90 | * \param[in] a first number. | ||
| 91 | * \param[in] b the number to subtract from \p a. | ||
| 92 | * \param[out] x high-magnitude component of the result. | ||
| 93 | * \param[out] y low-magnitude component of the result. | ||
| 94 | * \relates expansion | ||
| 95 | */ | ||
| 96 | 97789217 | inline void two_diff(double a, double b, double& x, double& y) { | |
| 97 | 97789217 | x = a - b; | |
| 98 | 97789217 | double bvirt = a - x; | |
| 99 | 97789217 | double avirt = x + bvirt; | |
| 100 | 97789217 | double bround = bvirt - b; | |
| 101 | 97789217 | double around = a - avirt; | |
| 102 | 97789217 | y = around + bround; | |
| 103 | 97789217 | } | |
| 104 | |||
| 105 | /** | ||
| 106 | * \brief Computes the sum of two doubles into a length 2 expansion. | ||
| 107 | * \details By Jonathan Shewchuk. | ||
| 108 | * \param[in] a first argument | ||
| 109 | * \param[in] b second argument | ||
| 110 | * \param[out] x high-magnitude component of the result | ||
| 111 | * \param[out] y low-magnitude component of the result | ||
| 112 | * \pre |\p a| > |\p b| | ||
| 113 | */ | ||
| 114 | 256858842 | inline void fast_two_sum(double a, double b, double& x, double& y) { | |
| 115 | 256858842 | x = a + b; | |
| 116 | 256858842 | double bvirt = x - a; | |
| 117 | 256858842 | y = b - bvirt; | |
| 118 | 256858842 | } | |
| 119 | |||
| 120 | /** | ||
| 121 | * \brief Computes the difference of two doubles into a length 2 expansion. | ||
| 122 | * \details By Jonathan Shewchuk. | ||
| 123 | * \param[in] a first argument | ||
| 124 | * \param[in] b second argument | ||
| 125 | * \param[out] x high-magnitude component of the result | ||
| 126 | * \param[out] y low-magnitude component of the result | ||
| 127 | * \pre | \p a| > | \p b | | ||
| 128 | */ | ||
| 129 | inline void fast_two_diff(double a, double b, double& x, double& y) { | ||
| 130 | x = a - b; | ||
| 131 | double bvirt = a - x; | ||
| 132 | y = bvirt - b; | ||
| 133 | } | ||
| 134 | |||
| 135 | /** | ||
| 136 | * \brief Splits a number into two components, ready for | ||
| 137 | * computing a product. | ||
| 138 | * \details By Jonathan Shewchuk. | ||
| 139 | * \param[in] a input number. | ||
| 140 | * \param[out] ahi split number, high-magnitude part. | ||
| 141 | * \param[out] alo split number, low-magnitude part. | ||
| 142 | * \relates expansion | ||
| 143 | */ | ||
| 144 | inline void split(double a, double& ahi, double& alo) { | ||
| 145 | geo_debug_assert(expansion_initialized_); | ||
| 146 | double c = expansion_splitter_ * a; | ||
| 147 | double abig = c - a; | ||
| 148 | ahi = c - abig; | ||
| 149 | alo = a - ahi; | ||
| 150 | } | ||
| 151 | |||
| 152 | /** | ||
| 153 | * \brief Multiplies two doubles into a length 2 expansion. | ||
| 154 | * \details By Jonathan Shewchuk. | ||
| 155 | * \param[in] a first number to multiply. | ||
| 156 | * \param[in] b second number to multiply. | ||
| 157 | * \param[out] x high-magnitude component of the result. | ||
| 158 | * \param[out] y low-magnitude component of the result. | ||
| 159 | * \relates expansion | ||
| 160 | */ | ||
| 161 | 452119389 | inline void two_product(double a, double b, 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 | 452119389 | x = a*b; | |
| 171 | 452119389 | y = fma(a,b,-x); | |
| 172 | #else | ||
| 173 | x = a * b; | ||
| 174 | double ahi, alo; | ||
| 175 | split(a, ahi, alo); | ||
| 176 | double bhi, blo; | ||
| 177 | split(b, bhi, blo); | ||
| 178 | double err1 = x - (ahi * bhi); | ||
| 179 | double err2 = err1 - (alo * bhi); | ||
| 180 | double err3 = err2 - (ahi * blo); | ||
| 181 | y = (alo * blo) - err3; | ||
| 182 | #endif | ||
| 183 | 452119389 | } | |
| 184 | |||
| 185 | /** | ||
| 186 | * \brief Squares a number into a length 2 expansion. | ||
| 187 | * \details By Jonathan Shewchuk. | ||
| 188 | * \param[in] a number to square. | ||
| 189 | * \param[out] x high-magnitude component of the result. | ||
| 190 | * \param[out] y low-magnitude component of the result. | ||
| 191 | * \relates expansion | ||
| 192 | */ | ||
| 193 | 12106188 | inline void square(double a, double& x, double& y) { | |
| 194 | #ifdef FP_FAST_FMA | ||
| 195 | // If the target processor supports the FMA (Fused Multiply Add) | ||
| 196 | // instruction, then the product of two doubles into a length-2 | ||
| 197 | // expansion can be implemented as follows. Thanks to Marc Glisse | ||
| 198 | // for the information. | ||
| 199 | // Note: under gcc, automatic generations of fma() for a*b+c needs | ||
| 200 | // to be deactivated, using -ffp-contract=off, else it may break | ||
| 201 | // other functions such as fast_expansion_sum_zeroelim(). | ||
| 202 | 12106188 | x = a*a; | |
| 203 | 12106188 | y = fma(a,a,-x); | |
| 204 | #else | ||
| 205 | x = a * a; | ||
| 206 | double ahi, alo; | ||
| 207 | split(a, ahi, alo); | ||
| 208 | double err1 = x - (ahi * ahi); | ||
| 209 | double err3 = err1 - ((ahi + ahi) * alo); | ||
| 210 | y = (alo * alo) - err3; | ||
| 211 | #endif | ||
| 212 | 12106188 | } | |
| 213 | |||
| 214 | /************************************************************************/ | ||
| 215 | |||
| 216 | /** | ||
| 217 | * \brief Represents numbers in arbitrary precision with a low-level API. | ||
| 218 | * \details The three basic operations sum, difference and product | ||
| 219 | * are implemented, as well as some geometric functions | ||
| 220 | * (squared distance and dot product). The sign of | ||
| 221 | * an expansion can be exactly computed. expansion | ||
| 222 | * is useful to implement exact geometric predicates. | ||
| 223 | * Some of Jonathan Shewchuk's expansion manipulation functions | ||
| 224 | * are used. | ||
| 225 | * A higher-level (but less efficient) interface is available | ||
| 226 | * through the \ref expansion_nt class (expansion number type, that | ||
| 227 | * overloads operators). | ||
| 228 | */ | ||
| 229 | class GEOGRAM_API expansion { | ||
| 230 | public: | ||
| 231 | /** | ||
| 232 | * \brief Gets the length of this expansion. | ||
| 233 | * \return the number of components of this expansion | ||
| 234 | */ | ||
| 235 | 2617112242 | index_t length() const { | |
| 236 | 2617112242 | return length_; | |
| 237 | } | ||
| 238 | |||
| 239 | /** | ||
| 240 | * \brief Gets the capacity of this expansion. | ||
| 241 | * \return the maximum number of components | ||
| 242 | * that can be stored in this expansion | ||
| 243 | */ | ||
| 244 | 574653117 | index_t capacity() const { | |
| 245 | 574653117 | return capacity_; | |
| 246 | } | ||
| 247 | |||
| 248 | /** | ||
| 249 | * \brief Changes the length of an expansion. | ||
| 250 | * \param[in] new_length new length of the expansion | ||
| 251 | * \pre new_length <= capacity() | ||
| 252 | */ | ||
| 253 | 310744258 | void set_length(index_t new_length) { | |
| 254 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 310744258 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
310744258 | geo_debug_assert(new_length <= capacity()); |
| 255 | 310744258 | length_ = new_length; | |
| 256 | 310744258 | } | |
| 257 | |||
| 258 | /** | ||
| 259 | * \brief Low level access to a component. | ||
| 260 | * \return a const reference to the \p i%th component | ||
| 261 | * of this expansion | ||
| 262 | */ | ||
| 263 | 1507629840 | const double& operator[] (index_t i) const { | |
| 264 | // Note: we allocate capacity+1 storage | ||
| 265 | // systematically, since basic functions | ||
| 266 | // may access one additional value (without | ||
| 267 | // using it) | ||
| 268 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1507629840 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1507629840 | geo_debug_assert(i <= capacity_); |
| 269 | 1507629840 | return x_[i]; | |
| 270 | } | ||
| 271 | |||
| 272 | /** | ||
| 273 | * \brief Low level access to a component. | ||
| 274 | * \return a reference to the \p i%th component | ||
| 275 | * of this expansion | ||
| 276 | */ | ||
| 277 | 508389554 | double& operator[] (index_t i) { | |
| 278 | // Note: we allocate capacity+1 storage | ||
| 279 | // systematically, since basic functions | ||
| 280 | // may access one additional value (without | ||
| 281 | // using it) | ||
| 282 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 508389554 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
508389554 | geo_debug_assert(i <= capacity_); |
| 283 | 508389554 | return x_[i]; | |
| 284 | } | ||
| 285 | |||
| 286 | /** | ||
| 287 | * \brief Low level access to the array of components. | ||
| 288 | * \return a pointer to the array that stores | ||
| 289 | * the components | ||
| 290 | */ | ||
| 291 | double* data() { | ||
| 292 | return x_; | ||
| 293 | } | ||
| 294 | |||
| 295 | /** | ||
| 296 | * \brief Low level access to the array of components. | ||
| 297 | * \return a const pointer to the array that stores | ||
| 298 | * the components | ||
| 299 | */ | ||
| 300 | 128856893 | const double* data() const { | |
| 301 | 128856893 | return x_; | |
| 302 | } | ||
| 303 | |||
| 304 | /** | ||
| 305 | * \brief Computes the amount of memory required to store | ||
| 306 | * an expansion. | ||
| 307 | * \param[in] capa the required capacity | ||
| 308 | * \return the total number of bytes required to store | ||
| 309 | * an expansion of capacity \p capa. | ||
| 310 | */ | ||
| 311 | 296637731 | static size_t bytes(index_t capa) { | |
| 312 | // --> 2*sizeof(double) because x_ is declared of size [2] | ||
| 313 | // to avoid compiler's warning. | ||
| 314 | // --> capa+1 to have an additional 'sentry' at the end | ||
| 315 | // because fast_expansion_sum_zeroelim() may access | ||
| 316 | // an entry past the end (without using it). | ||
| 317 | return | ||
| 318 | 296637731 | sizeof(expansion) - 2 * sizeof(double) + | |
| 319 | 296637731 | (capa + 1) * sizeof(double); | |
| 320 | } | ||
| 321 | |||
| 322 | /** | ||
| 323 | * \brief Computes the amount of memory required to store | ||
| 324 | * an expansion on the stack | ||
| 325 | * \details Behaves like bytes() but in debug mode checks | ||
| 326 | * that this will fit on the stack. | ||
| 327 | * \param[in] capa the required capacity | ||
| 328 | * \return the total number of bytes required to store | ||
| 329 | * an expansion of capacity \p capa. | ||
| 330 | */ | ||
| 331 | 265021983 | static size_t bytes_on_stack(index_t capa) { | |
| 332 | #ifndef GEO_HAS_BIG_STACK | ||
| 333 | // Note: standard predicates need at least 512, hence the min. | ||
| 334 | // index_t(MAX_CAPACITY_ON_STACK) is necessary, else with | ||
| 335 | // MAX_CAPACITY_ON_STACK alone the compiler tries to generate a | ||
| 336 | // reference to NOT_IN_LIST resulting in a link error. | ||
| 337 | // (weird, even with constexpr, I do not understand...) | ||
| 338 | // Probably when the function excepts a *reference* | ||
| 339 | geo_debug_assert( | ||
| 340 | capa <= std::max(index_t(MAX_CAPACITY_ON_STACK),index_t(512)) | ||
| 341 | ); | ||
| 342 | #endif | ||
| 343 | 265021983 | return bytes(capa); | |
| 344 | } | ||
| 345 | |||
| 346 | /** | ||
| 347 | * \brief Client code should not use this constructor. | ||
| 348 | * \details This constructor should not be used by client code, | ||
| 349 | * use new_expansion_on_stack() and new_expansion_on_heap() | ||
| 350 | * instead. This constructor initializes this expansion's length | ||
| 351 | * and capacity. Note that it should be created with enough space, | ||
| 352 | * using placement syntax of operator new. | ||
| 353 | * \param[in] capa the capacity | ||
| 354 | */ | ||
| 355 | 296637731 | expansion(index_t capa) : | |
| 356 | 296637731 | length_(0), | |
| 357 | 296637731 | capacity_(capa) { | |
| 358 | 296637731 | } | |
| 359 | |||
| 360 | /** | ||
| 361 | * \brief Allocates an expansion on the stack. | ||
| 362 | * \details It can only be a macro (and not an inline function) | ||
| 363 | * since alloca() cannot be called from inline functions. | ||
| 364 | * \arg capa required capacity of the expansion. | ||
| 365 | * \relates GEO::expansion | ||
| 366 | */ | ||
| 367 | #ifdef CPPCHECK | ||
| 368 | // cppcheck does not understand that the result | ||
| 369 | // of alloca() is passed to the placement syntax | ||
| 370 | // of operator new. | ||
| 371 | expansion& new_expansion_on_stack(index_t capa); | ||
| 372 | #else | ||
| 373 | #define new_expansion_on_stack(capa) \ | ||
| 374 | (new (alloca(expansion::bytes_on_stack(capa)))expansion(capa)) | ||
| 375 | #endif | ||
| 376 | |||
| 377 | /** | ||
| 378 | * \brief Allocates an expansion on the heap. | ||
| 379 | * \details Allocates also some space for the reference counter. | ||
| 380 | * \param[in] capa capacity (i.e. maximum length) of the expansion. | ||
| 381 | * \return a pointer to the newly allocated expansion | ||
| 382 | */ | ||
| 383 | 31615748 | static inline expansion* new_expansion_on_heap(index_t capa) { | |
| 384 | 31615748 | void* addr = malloc(bytes(capa)); | |
| 385 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 31615748 times.
|
31615748 | return new(addr)expansion(capa); |
| 386 | } | ||
| 387 | |||
| 388 | /** | ||
| 389 | * \brief Deallocates an expansion on the heap. | ||
| 390 | * \param[in] e the expansion | ||
| 391 | * \pre \p e should have been previously allocated | ||
| 392 | * by new_expansion_on_heap() | ||
| 393 | */ | ||
| 394 | 31615748 | static inline void delete_expansion_on_heap(expansion* e) { | |
| 395 | 31615748 | free(e); | |
| 396 | 31615748 | } | |
| 397 | |||
| 398 | // ========================== Initialization from doubles | ||
| 399 | |||
| 400 | /** | ||
| 401 | * \brief Assigns a number to this expansion. | ||
| 402 | * \param[in] a the number | ||
| 403 | * \return the new value of this expansion (\p a) | ||
| 404 | */ | ||
| 405 | 12132 | expansion& assign(double a) { | |
| 406 | 12132 | set_length(1); | |
| 407 | 12132 | x_[0] = a; | |
| 408 | 12132 | return *this; | |
| 409 | } | ||
| 410 | |||
| 411 | /** | ||
| 412 | * \brief Copies an expansion to this expansion | ||
| 413 | * \param[in] rhs the expansion to be copied | ||
| 414 | * \return the new value of this expansion (\p rhs) | ||
| 415 | */ | ||
| 416 | 222771 | expansion& assign(const expansion& rhs) { | |
| 417 |
1/6✗ Branch 2 not taken.
✓ Branch 3 taken 222771 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
222771 | geo_debug_assert(capacity() >= rhs.length()); |
| 418 | 222771 | set_length(rhs.length()); | |
| 419 |
2/2✓ Branch 1 taken 475660 times.
✓ Branch 2 taken 222771 times.
|
698431 | for(index_t i=0; i<rhs.length(); ++i) { |
| 420 | 475660 | x_[i] = rhs.x_[i]; | |
| 421 | } | ||
| 422 | 222771 | return *this; | |
| 423 | } | ||
| 424 | |||
| 425 | /** | ||
| 426 | * \brief Copies the absolute value of an expansion to this expansion | ||
| 427 | * \param[in] rhs the expansion to be copied | ||
| 428 | * \return the new value of this expansion that is, abs(\p rhs) | ||
| 429 | */ | ||
| 430 | expansion& assign_abs(const expansion& rhs) { | ||
| 431 | assign(rhs); | ||
| 432 | if(sign() == NEGATIVE) { | ||
| 433 | negate(); | ||
| 434 | } | ||
| 435 | return *this; | ||
| 436 | } | ||
| 437 | |||
| 438 | /** | ||
| 439 | * \brief Computes the required capacity to store the | ||
| 440 | * sum of two doubles. | ||
| 441 | * \param[in] a first number | ||
| 442 | * \param[in] b second number | ||
| 443 | * \return the required capacity of an expansion | ||
| 444 | * to store the exact sum of two doubles | ||
| 445 | * \note The result does not depend on the values of the | ||
| 446 | * two numbers. | ||
| 447 | */ | ||
| 448 | ✗ | static index_t sum_capacity(double a, double b) { | |
| 449 | ✗ | geo_argused(a); | |
| 450 | ✗ | geo_argused(b); | |
| 451 | ✗ | return 2; | |
| 452 | } | ||
| 453 | |||
| 454 | /** | ||
| 455 | * \brief Assigns the sum of two doubles to this expansion | ||
| 456 | * (should not be used by client code). | ||
| 457 | * \details Do not use directly, | ||
| 458 | * use expansion_sum() macro instead. | ||
| 459 | * \param[in] a first number to sum | ||
| 460 | * \param[in] b second number to sum | ||
| 461 | * \return the new value of this expansion (\p a + \p b) | ||
| 462 | * \pre capacity() >= sum_capacity(a,b) | ||
| 463 | */ | ||
| 464 | ✗ | expansion& assign_sum(double a, double b) { | |
| 465 | ✗ | set_length(2); | |
| 466 | ✗ | two_sum(a, b, x_[1], x_[0]); | |
| 467 | ✗ | return *this; | |
| 468 | } | ||
| 469 | |||
| 470 | /** | ||
| 471 | * \brief Computes the required capacity of an expansion | ||
| 472 | * to store the exact difference of two doubles. | ||
| 473 | * \param[in] a first number | ||
| 474 | * \param[in] b second number | ||
| 475 | * \return the required capacity of an expansion | ||
| 476 | * to store the exact difference of two doubles | ||
| 477 | * \note The result does not depend on the values of the | ||
| 478 | * two numbers. | ||
| 479 | */ | ||
| 480 | 154957819 | static index_t diff_capacity(double a, double b) { | |
| 481 | 154957819 | geo_argused(a); | |
| 482 | 154957819 | geo_argused(b); | |
| 483 | 154957819 | return 2; | |
| 484 | } | ||
| 485 | |||
| 486 | /** | ||
| 487 | * \brief Assigns the difference of two doubles to this expansion | ||
| 488 | * (should not be used by client code). | ||
| 489 | * \details Do not use directly, | ||
| 490 | * use expansion_diff() macro instead. | ||
| 491 | * \param[in] a first number | ||
| 492 | * \param[in] b second number | ||
| 493 | * \return the new value of this expansion (\p a - \p b) | ||
| 494 | * \pre capacity() >= diff_capacity(a,b) | ||
| 495 | */ | ||
| 496 | 79422929 | expansion& assign_diff(double a, double b) { | |
| 497 | 79422929 | set_length(2); | |
| 498 | 79422929 | two_diff(a, b, x_[1], x_[0]); | |
| 499 | 79422929 | return *this; | |
| 500 | } | ||
| 501 | |||
| 502 | /** | ||
| 503 | * \brief Computes the required capacity of an expansion | ||
| 504 | * to store the exact product of two doubles. | ||
| 505 | * \param[in] a first number | ||
| 506 | * \param[in] b second number | ||
| 507 | * \return the required capacity of an expansion | ||
| 508 | * to store the exact product of two doubles | ||
| 509 | * \note The result does not depend on the values of the | ||
| 510 | * two numbers. | ||
| 511 | */ | ||
| 512 | ✗ | static index_t product_capacity(double a, double b) { | |
| 513 | ✗ | geo_argused(a); | |
| 514 | ✗ | geo_argused(b); | |
| 515 | ✗ | return 2; | |
| 516 | } | ||
| 517 | |||
| 518 | /** | ||
| 519 | * \brief Assigns the product of two doubles to this expansion | ||
| 520 | * (should not be used by client code). | ||
| 521 | * \details Do not use directly, | ||
| 522 | * use expansion_product() macro instead. | ||
| 523 | * \param[in] a first number | ||
| 524 | * \param[in] b second number | ||
| 525 | * \return the new value of this expansion (\p a * \p b) | ||
| 526 | * \pre capacity() >= product_capacity(a,b) | ||
| 527 | */ | ||
| 528 | ✗ | expansion& assign_product(double a, double b) { | |
| 529 | ✗ | set_length(2); | |
| 530 | ✗ | two_product(a, b, x_[1], x_[0]); | |
| 531 | ✗ | return *this; | |
| 532 | } | ||
| 533 | |||
| 534 | /** | ||
| 535 | * \brief Computes the required capacity of an expansion | ||
| 536 | * to store the exact square of a double. | ||
| 537 | * \param[in] a the number to be squared | ||
| 538 | * \return the required capacity of an expansion | ||
| 539 | * to store the exact square of a double. | ||
| 540 | * \note The result does not depend on the value of the | ||
| 541 | * number \p a. | ||
| 542 | */ | ||
| 543 | static index_t square_capacity(double a) { | ||
| 544 | geo_argused(a); | ||
| 545 | return 2; | ||
| 546 | } | ||
| 547 | |||
| 548 | /** | ||
| 549 | * \brief Assigns the square of a double to this expansion | ||
| 550 | * (should not be used by client code). | ||
| 551 | * \details Do not use directly, | ||
| 552 | * use expansion_square() macro instead. | ||
| 553 | * \param[in] a the number to be squared | ||
| 554 | * \return the new value of this expansion (\p a * \p a) | ||
| 555 | * \pre capacity() >= square_capacity(a) | ||
| 556 | */ | ||
| 557 | expansion& assign_square(double a) { | ||
| 558 | set_length(2); | ||
| 559 | square(a, x_[1], x_[0]); | ||
| 560 | return *this; | ||
| 561 | } | ||
| 562 | |||
| 563 | // ====== Initialization from expansion and double | ||
| 564 | |||
| 565 | /** | ||
| 566 | * \brief Computes the required capacity of an expansion | ||
| 567 | * to store the exact sum of an expansion and a double. | ||
| 568 | * \param[in] a first number as an expansion | ||
| 569 | * \param[in] b second number as a double | ||
| 570 | * \return the required capacity of an expansion | ||
| 571 | * to store the exact sum of \p a and \p b | ||
| 572 | * \note The result does not depend on the value of the | ||
| 573 | * double argument \p b. | ||
| 574 | */ | ||
| 575 | ✗ | static index_t sum_capacity(const expansion& a, double b) { | |
| 576 | ✗ | geo_argused(b); | |
| 577 | ✗ | return a.length() + 1; | |
| 578 | } | ||
| 579 | |||
| 580 | /** | ||
| 581 | * \brief Assigns the sum of an expansion and a double | ||
| 582 | * to this expansion (should not be used by client code). | ||
| 583 | * \details Do not use directly, | ||
| 584 | * use expansion_sum() macro instead. | ||
| 585 | * \param[in] a the expansion | ||
| 586 | * \param[in] b the double | ||
| 587 | * \return the new value of this expansion (\p a + \p b) | ||
| 588 | * \pre capacity() >= sum_capacity(a,b) | ||
| 589 | */ | ||
| 590 | expansion& assign_sum(const expansion& a, double b); | ||
| 591 | |||
| 592 | /** | ||
| 593 | * \brief Computes the required capacity of an expansion | ||
| 594 | * to store the exact difference between an expansion and a double. | ||
| 595 | * \param[in] a first number as an expansion | ||
| 596 | * \param[in] b second number as a double | ||
| 597 | * \return the required capacity of an expansion | ||
| 598 | * to store the exact difference of \p a and \p b | ||
| 599 | * \note The result does not depend on the value of the | ||
| 600 | * double argument \p b. | ||
| 601 | */ | ||
| 602 | ✗ | static index_t diff_capacity(const expansion& a, double b) { | |
| 603 | ✗ | geo_argused(b); | |
| 604 | ✗ | return a.length() + 1; | |
| 605 | } | ||
| 606 | |||
| 607 | /** | ||
| 608 | * \brief Assigns the difference between an expansion and a double | ||
| 609 | * to this expansion (should not be used by client code). | ||
| 610 | * \details Do not use directly, | ||
| 611 | * use expansion_diff() macro instead. | ||
| 612 | * \param[in] a the expansion | ||
| 613 | * \param[in] b the double | ||
| 614 | * \return the new value of this expansion (\p a - \p b) | ||
| 615 | * \pre capacity() >= diff_capacity(a,b) | ||
| 616 | */ | ||
| 617 | expansion& assign_diff(const expansion& a, double b); | ||
| 618 | |||
| 619 | /** | ||
| 620 | * \brief Computes the required capacity of an expansion | ||
| 621 | * to store the exact product between an expansion and a double. | ||
| 622 | * \param[in] a the expansion | ||
| 623 | * \param[in] b the double | ||
| 624 | * \return the required capacity of an expansion to store | ||
| 625 | * the exact product between \p a and \p b | ||
| 626 | * \note The result does not depend on the value of the | ||
| 627 | * double argument \p b. | ||
| 628 | */ | ||
| 629 | 42777912 | static index_t product_capacity(const expansion& a, double b) { | |
| 630 | 42777912 | geo_argused(b); | |
| 631 | // TODO: implement special case where the double argument | ||
| 632 | // is a power of two. | ||
| 633 | 42777912 | return a.length() * 2; | |
| 634 | } | ||
| 635 | |||
| 636 | /** | ||
| 637 | * \brief Assigns the product between an expansion and a double | ||
| 638 | * to this expansion (should not be used by client code). | ||
| 639 | * \details Do not use directly, | ||
| 640 | * use expansion_product() macro instead. | ||
| 641 | * \param[in] a the expansion | ||
| 642 | * \param[in] b the double | ||
| 643 | * \return the new value of this expansion (\p a * \p b) | ||
| 644 | * \pre capacity() >= product_capacity(a,b) | ||
| 645 | */ | ||
| 646 | expansion& assign_product(const expansion& a, double b); | ||
| 647 | |||
| 648 | // ========================== Initialization from expansions | ||
| 649 | |||
| 650 | /** | ||
| 651 | * \brief Computes the required capacity of an expansion | ||
| 652 | * to store the exact sum of two expansions. | ||
| 653 | * \param[in] a first expansion | ||
| 654 | * \param[in] b second expansion | ||
| 655 | * \return the required capacity of an expansion | ||
| 656 | * to store the exact sum of \p a and \p b | ||
| 657 | */ | ||
| 658 | 62913619 | static index_t sum_capacity(const expansion& a, const expansion& b) { | |
| 659 | 62913619 | return a.length() + b.length(); | |
| 660 | } | ||
| 661 | |||
| 662 | /** | ||
| 663 | * \brief Assigns the sum of two expansions | ||
| 664 | * to this expansion (should not be used by client code). | ||
| 665 | * \details Do not use directly, | ||
| 666 | * use expansion_sum() macro instead. | ||
| 667 | * \param[in] a first expansion | ||
| 668 | * \param[in] b second expansion | ||
| 669 | * \return the new value of this expansion (\p a + \p b) | ||
| 670 | * \pre capacity() >= sum_capacity(a,b) | ||
| 671 | */ | ||
| 672 | expansion& assign_sum(const expansion& a, const expansion& b); | ||
| 673 | |||
| 674 | /** | ||
| 675 | * \brief Computes the required capacity of an expansion | ||
| 676 | * to store the exact sum of three expansions. | ||
| 677 | * \param[in] a first expansion | ||
| 678 | * \param[in] b second expansion | ||
| 679 | * \param[in] c third expansion | ||
| 680 | * \return the required capacity of an expansion | ||
| 681 | * to store the exact sum of \p a, \p b and \p c | ||
| 682 | */ | ||
| 683 | 12215372 | static index_t sum_capacity( | |
| 684 | const expansion& a, const expansion& b, const expansion& c | ||
| 685 | ) { | ||
| 686 | 12215372 | return a.length() + b.length() + c.length(); | |
| 687 | } | ||
| 688 | |||
| 689 | /** | ||
| 690 | * \brief Assigns the sum of three expansions | ||
| 691 | * to this expansion (should not be used by client code). | ||
| 692 | * \details Do not use directly, | ||
| 693 | * use expansion_sum3() macro instead. | ||
| 694 | * \param[in] a first expansion | ||
| 695 | * \param[in] b second expansion | ||
| 696 | * \param[in] c third expansion | ||
| 697 | * \return the new value of this expansion (\p a + \p b + \p c) | ||
| 698 | * \pre capacity() >= sum_capacity(a,b,c) | ||
| 699 | */ | ||
| 700 | expansion& assign_sum( | ||
| 701 | const expansion& a, const expansion& b, const expansion& c | ||
| 702 | ); | ||
| 703 | |||
| 704 | /** | ||
| 705 | * \brief Computes the required capacity of an expansion | ||
| 706 | * to store the exact sum of four expansions. | ||
| 707 | * \param[in] a first expansion | ||
| 708 | * \param[in] b second expansion | ||
| 709 | * \param[in] c third expansion | ||
| 710 | * \param[in] d fourth expansion | ||
| 711 | * \return the required capacity of an expansion | ||
| 712 | * to store the exact sum of \p a, \p b, \p c and \p d | ||
| 713 | */ | ||
| 714 | 748836 | static index_t sum_capacity( | |
| 715 | const expansion& a, const expansion& b, | ||
| 716 | const expansion& c, const expansion& d | ||
| 717 | ) { | ||
| 718 | 748836 | return a.length() + b.length() + c.length() + d.length(); | |
| 719 | } | ||
| 720 | |||
| 721 | /** | ||
| 722 | * \brief Assigns the sum of four expansions | ||
| 723 | * to this expansion (should not be used by client code). | ||
| 724 | * \details Do not use directly, | ||
| 725 | * use expansion_sum4() macro instead. | ||
| 726 | * \param[in] a first expansion | ||
| 727 | * \param[in] b second expansion | ||
| 728 | * \param[in] c third expansion | ||
| 729 | * \param[in] d fourth expansion | ||
| 730 | * \return the new value of this expansion (\p a + \p b + \p c + \p d) | ||
| 731 | * \pre capacity() >= sum_capacity(a,b,c,d) | ||
| 732 | */ | ||
| 733 | expansion& assign_sum( | ||
| 734 | const expansion& a, const expansion& b, | ||
| 735 | const expansion& c, const expansion& d | ||
| 736 | ); | ||
| 737 | |||
| 738 | /** | ||
| 739 | * \brief Computes the required capacity of an expansion | ||
| 740 | * to store the exact difference of two expansions. | ||
| 741 | * \param[in] a first expansion | ||
| 742 | * \param[in] b second expansion | ||
| 743 | * \return the required capacity of an expansion | ||
| 744 | * to store the exact difference between \p a and \p b | ||
| 745 | */ | ||
| 746 | 51302721 | static index_t diff_capacity(const expansion& a, const expansion& b) { | |
| 747 | 51302721 | return a.length() + b.length(); | |
| 748 | } | ||
| 749 | |||
| 750 | /** | ||
| 751 | * \brief Assigns the difference between two expansions | ||
| 752 | * to this expansion (should not be used by client code). | ||
| 753 | * \details Do not use directly, | ||
| 754 | * use expansion_diff() macro instead. | ||
| 755 | * \param[in] a first expansion | ||
| 756 | * \param[in] b second expansion | ||
| 757 | * \return the new value of this expansion (\p a - \p b) | ||
| 758 | * \pre capacity() >= diff_capacity(a,b) | ||
| 759 | */ | ||
| 760 | expansion& assign_diff(const expansion& a, const expansion& b); | ||
| 761 | |||
| 762 | /** | ||
| 763 | * \brief Computes the required capacity of an expansion | ||
| 764 | * to store the exact product of two expansions. | ||
| 765 | * \param[in] a first expansion | ||
| 766 | * \param[in] b second expansion | ||
| 767 | * \return the required capacity of an expansion | ||
| 768 | * to store the exact product of \p a and \p b | ||
| 769 | */ | ||
| 770 | 485414003 | static index_t product_capacity( | |
| 771 | const expansion& a, const expansion& b | ||
| 772 | ) { | ||
| 773 | 485414003 | return a.length() * b.length() * 2; | |
| 774 | } | ||
| 775 | |||
| 776 | /** | ||
| 777 | * \brief Assigns the product of two expansions | ||
| 778 | * to this expansion (should not be used by client code). | ||
| 779 | * \details Do not use directly, | ||
| 780 | * use expansion_product() macro instead. | ||
| 781 | * \param[in] a first expansion | ||
| 782 | * \param[in] b second expansion | ||
| 783 | * \return the new value of this expansion (\p a * \p b) | ||
| 784 | * \pre capacity() >= product_capacity(a,b) | ||
| 785 | */ | ||
| 786 | expansion& assign_product(const expansion& a, const expansion& b); | ||
| 787 | |||
| 788 | /** | ||
| 789 | * \brief Computes the required capacity of an expansion | ||
| 790 | * to store the exact product of three expansions. | ||
| 791 | * \param[in] a first expansion | ||
| 792 | * \param[in] b second expansion | ||
| 793 | * \param[in] c third expansion | ||
| 794 | * \return the required capacity of an expansion | ||
| 795 | * to store the exact product of \p a, \p b and \p c | ||
| 796 | */ | ||
| 797 | ✗ | static index_t product_capacity( | |
| 798 | const expansion& a, const expansion& b, const expansion& c | ||
| 799 | ) { | ||
| 800 | ✗ | return a.length() * b.length() * c.length() * 4; | |
| 801 | } | ||
| 802 | |||
| 803 | /** | ||
| 804 | * \brief Assigns the product of three expansions | ||
| 805 | * to this expansion (should not be used by client code). | ||
| 806 | * \details Do not use directly, | ||
| 807 | * use expansion_product3() macro instead. | ||
| 808 | * \param[in] a first expansion | ||
| 809 | * \param[in] b second expansion | ||
| 810 | * \param[in] c third expansion | ||
| 811 | * \return the new value of this expansion (\p a * \p b * \p c) | ||
| 812 | * \pre capacity() >= product_capacity(a,b,c) | ||
| 813 | */ | ||
| 814 | expansion& assign_product( | ||
| 815 | const expansion& a, const expansion& b, const expansion& c | ||
| 816 | ); | ||
| 817 | |||
| 818 | /** | ||
| 819 | * \brief Computes the required capacity of an expansion | ||
| 820 | * to store the exact square of an expansion. | ||
| 821 | * \param[in] a the expansion to be squared | ||
| 822 | * \return the required capacity of an expansion | ||
| 823 | * to store the exact square \p a * \p a | ||
| 824 | */ | ||
| 825 | ✗ | static index_t square_capacity(const expansion& a) { | |
| 826 | ✗ | if(a.length() == 2) { | |
| 827 | ✗ | return 6; | |
| 828 | } // see two_square() | ||
| 829 | ✗ | return a.length() * a.length() * 2; | |
| 830 | } | ||
| 831 | |||
| 832 | /** | ||
| 833 | * \brief Assigns the product of an expansions | ||
| 834 | * to this expansion (should not be used by client code). | ||
| 835 | * \details Do not use directly, | ||
| 836 | * use expansion_square() macro instead. | ||
| 837 | * \param[in] a the expansion to be squared | ||
| 838 | * \return the new value of this expansion (\p a * \p a ) | ||
| 839 | * \pre capacity() >= square_capacity(a) | ||
| 840 | */ | ||
| 841 | expansion& assign_square(const expansion& a); | ||
| 842 | |||
| 843 | // ====== Determinants ============================= | ||
| 844 | |||
| 845 | /** | ||
| 846 | * \brief Computes the required capacity of an expansion | ||
| 847 | * to store an exact 2x2 determinant. | ||
| 848 | * \param[in] a11 , a12 , a21 , a22 coefficients of the determinant | ||
| 849 | * \return the required capacity of an expansion to store | ||
| 850 | * the exact determinant \p a11 * \p a22 - \p a21 * \p a12 | ||
| 851 | */ | ||
| 852 | 99545656 | static index_t det2x2_capacity( | |
| 853 | const expansion& a11, const expansion& a12, | ||
| 854 | const expansion& a21, const expansion& a22 | ||
| 855 | ) { | ||
| 856 | return | ||
| 857 | 99545656 | product_capacity(a11, a22) + | |
| 858 | 99545656 | product_capacity(a21, a12); | |
| 859 | } | ||
| 860 | |||
| 861 | /** | ||
| 862 | * \brief Assigns a 2x2 determinant to this expansion | ||
| 863 | * (should not be used by client code). | ||
| 864 | * \details Do not use directly, use expansion_det2x2() | ||
| 865 | * macro instead. | ||
| 866 | * \param[in] a11 , a12 , a21 , a22 coefficients of the determinant | ||
| 867 | * \return the new value of this expansion, with | ||
| 868 | * the exact determinant \p a11 * \p a22 - \p a21 * \p a12 | ||
| 869 | * \pre capacity() >= det_2x2_capacity(a11,a12,,a21,a22) | ||
| 870 | */ | ||
| 871 | expansion& assign_det2x2( | ||
| 872 | const expansion& a11, const expansion& a12, | ||
| 873 | const expansion& a21, const expansion& a22 | ||
| 874 | ); | ||
| 875 | |||
| 876 | /** | ||
| 877 | * \brief Computes the required capacity of an expansion | ||
| 878 | * to store an exact 3x3 determinant. | ||
| 879 | * \param[in] a11 , a12 , a13 , a21 , a22 , a23 , a31 , a32 , a33 | ||
| 880 | * coefficients of the determinant | ||
| 881 | * \return the required capacity of an expansion to store | ||
| 882 | * the exact value of the determinant | ||
| 883 | */ | ||
| 884 | 12007104 | static index_t det3x3_capacity( | |
| 885 | const expansion& a11, const expansion& a12, const expansion& a13, | ||
| 886 | const expansion& a21, const expansion& a22, const expansion& a23, | ||
| 887 | const expansion& a31, const expansion& a32, const expansion& a33 | ||
| 888 | ) { | ||
| 889 | // Development w.r.t. first row | ||
| 890 | 12007104 | index_t c11_capa = det2x2_capacity(a22, a23, a32, a33); | |
| 891 | 12007104 | index_t c12_capa = det2x2_capacity(a21, a23, a31, a33); | |
| 892 | 12007104 | index_t c13_capa = det2x2_capacity(a21, a22, a31, a32); | |
| 893 | return 2 * ( | ||
| 894 | 12007104 | a11.length() * c11_capa + | |
| 895 | 12007104 | a12.length() * c12_capa + | |
| 896 | 12007104 | a13.length() * c13_capa | |
| 897 | 12007104 | ); | |
| 898 | } | ||
| 899 | |||
| 900 | /** | ||
| 901 | * \brief Assigns a 3x3 determinant to this expansion | ||
| 902 | * (should not be used by client code). | ||
| 903 | * \details Do not use directly, use expansion_det3x3() | ||
| 904 | * macro instead. | ||
| 905 | * \param[in] a11 , a12 , a13 , a21 , a22 , a23 , a31 , a32 , a33 | ||
| 906 | * coefficients of the determinant | ||
| 907 | * \return the new value of this expansion, with | ||
| 908 | * the exact 3x3 determinant | ||
| 909 | * \pre capacity() >= | ||
| 910 | * det_3x3_capacity(a11,a12,a13,a21,a22,a23,a31,a32,a33) | ||
| 911 | */ | ||
| 912 | expansion& assign_det3x3( | ||
| 913 | const expansion& a11, const expansion& a12, const expansion& a13, | ||
| 914 | const expansion& a21, const expansion& a22, const expansion& a23, | ||
| 915 | const expansion& a31, const expansion& a32, const expansion& a33 | ||
| 916 | ); | ||
| 917 | |||
| 918 | /** | ||
| 919 | * \brief Computes the required capacity of an expansion | ||
| 920 | * to store an exact 3x3 determinant where the | ||
| 921 | * first row is 1 1 1. | ||
| 922 | * \param[in] a21 , a22 , a23 , a31 , a32 , a33 coefficients | ||
| 923 | * of the determinant | ||
| 924 | * \return the required capacity of an expansion to store | ||
| 925 | * the exact value of the determinant | ||
| 926 | */ | ||
| 927 | ✗ | static index_t det_111_2x3_capacity( | |
| 928 | const expansion& a21, const expansion& a22, const expansion& a23, | ||
| 929 | const expansion& a31, const expansion& a32, const expansion& a33 | ||
| 930 | ) { | ||
| 931 | return | ||
| 932 | ✗ | det2x2_capacity(a22, a23, a32, a33) + | |
| 933 | ✗ | det2x2_capacity(a23, a21, a33, a31) + | |
| 934 | ✗ | det2x2_capacity(a21, a22, a31, a32); | |
| 935 | } | ||
| 936 | |||
| 937 | /** | ||
| 938 | * \brief Assigns a 3x3 determinant to this expansion | ||
| 939 | * where the first row is 1 1 1(should not be used by client code). | ||
| 940 | * \details Do not use directly, use expansion_det_111_3x3() | ||
| 941 | * macro instead. | ||
| 942 | * \param[in] a21 , a22 , a23 , a31 , a32 , a33 coefficients | ||
| 943 | * of the determinant | ||
| 944 | * \return the new value of this expansion, with | ||
| 945 | * the exact 3x3 determinant | ||
| 946 | * \pre capacity() >= det__111_2x3capacity(a21,a22,a23,a31,a32,a33) | ||
| 947 | */ | ||
| 948 | expansion& assign_det_111_2x3( | ||
| 949 | const expansion& a21, const expansion& a22, const expansion& a23, | ||
| 950 | const expansion& a31, const expansion& a32, const expansion& a33 | ||
| 951 | ); | ||
| 952 | |||
| 953 | // ======= Geometry-specific initializations ======= | ||
| 954 | |||
| 955 | /** | ||
| 956 | * \brief Computes the required capacity of an expansion | ||
| 957 | * to store the exact squared distance between two points | ||
| 958 | * of specified dimension. | ||
| 959 | * \param[in] dim dimension of the points | ||
| 960 | * \return the required capacity of an expansion to store | ||
| 961 | * the exact squared distance between points of dimension \p dim | ||
| 962 | */ | ||
| 963 | 30265386 | static index_t sq_dist_capacity(coord_index_t dim) { | |
| 964 | 30265386 | return index_t(dim) * 6; | |
| 965 | } | ||
| 966 | |||
| 967 | /** | ||
| 968 | * \brief Assigns the squared distance between two points to | ||
| 969 | * this expansion (should not be used by client code). | ||
| 970 | * \details Do not use directly, | ||
| 971 | * use expansion_sq_dist() macro instead. | ||
| 972 | * \param[in] p1 pointer to the coordinates of the first point | ||
| 973 | * \param[in] p2 pointer to the coordinates of the second point | ||
| 974 | * \param[in] dim dimension of the points | ||
| 975 | * \return the new value of this expansion, with the squared distance | ||
| 976 | * between \p p1 and \p p2 | ||
| 977 | * \pre capacity() >= sq_dist_capacity(dim) | ||
| 978 | */ | ||
| 979 | expansion& assign_sq_dist( | ||
| 980 | const double* p1, const double* p2, coord_index_t dim | ||
| 981 | ); | ||
| 982 | |||
| 983 | /** | ||
| 984 | * \brief Computes the required capacity of an expansion | ||
| 985 | * to store the exact dot product between two vectors. | ||
| 986 | * \param[in] dim dimension of the vectors | ||
| 987 | * \return the required capacity of an expansion to store | ||
| 988 | * the dot product between two \p dim%-dimensional vectors | ||
| 989 | * specified by differences of doubles. | ||
| 990 | */ | ||
| 991 | 30782733 | static index_t dot_at_capacity(coord_index_t dim) { | |
| 992 | 30782733 | return index_t(dim) * 8; | |
| 993 | } | ||
| 994 | |||
| 995 | /** | ||
| 996 | * \brief Assigns the dot product of two vectors to | ||
| 997 | * this expansion (should not be used by client code). | ||
| 998 | * \details Do not use directly, | ||
| 999 | * use expansion_dot_at() macro instead. | ||
| 1000 | * \param[in] p1 pointer to the coordinates of a point | ||
| 1001 | * \param[in] p2 pointer to the coordinates of a point | ||
| 1002 | * \param[in] p0 pointer to the coordinates of a point | ||
| 1003 | * \param[in] dim dimension of the points | ||
| 1004 | * \return the new value of this expansion, with the dot | ||
| 1005 | * product (p1-p0).(p2-p0) | ||
| 1006 | */ | ||
| 1007 | expansion& assign_dot_at( | ||
| 1008 | const double* p1, const double* p2, const double* p0, | ||
| 1009 | coord_index_t dim | ||
| 1010 | ); | ||
| 1011 | |||
| 1012 | |||
| 1013 | /** | ||
| 1014 | * \brief Computes the required capacity to store the | ||
| 1015 | * length of a 3d vector. | ||
| 1016 | * \param[in] x , y , z coordinates of the vector | ||
| 1017 | * \return the capacity required to store the squared norm | ||
| 1018 | * of [x,y,z] | ||
| 1019 | */ | ||
| 1020 | static index_t length2_capacity( | ||
| 1021 | const expansion& x, const expansion& y, const expansion& z | ||
| 1022 | ) { | ||
| 1023 | return square_capacity(x) + square_capacity(y) + square_capacity(z); | ||
| 1024 | } | ||
| 1025 | |||
| 1026 | /** | ||
| 1027 | * \brief Assigns the length of a vector to this expansion | ||
| 1028 | * (should not be used by client code). Do not call this | ||
| 1029 | * function directly, use expansion_length2() macro instead. | ||
| 1030 | * \param[in] x , y , z coordinates of the vector | ||
| 1031 | * \return the new value of this expansion, with the squared | ||
| 1032 | * length of [x,y,z] | ||
| 1033 | */ | ||
| 1034 | expansion& assign_length2( | ||
| 1035 | const expansion& x, const expansion& y, const expansion& z | ||
| 1036 | ); | ||
| 1037 | |||
| 1038 | // =============== some general purpose functions ========= | ||
| 1039 | |||
| 1040 | /** | ||
| 1041 | * \brief Initializes the expansion class. | ||
| 1042 | * \details This function needs to be called once in the program, | ||
| 1043 | * before using any expansion object and operation (it computes | ||
| 1044 | * some internally-used constants). | ||
| 1045 | */ | ||
| 1046 | static void initialize(); | ||
| 1047 | |||
| 1048 | /** | ||
| 1049 | * \brief Changes the sign of an expansion. | ||
| 1050 | * \return the new value of this expansion | ||
| 1051 | */ | ||
| 1052 | 5054761 | expansion& negate() { | |
| 1053 |
2/2✓ Branch 0 taken 16709445 times.
✓ Branch 1 taken 5054761 times.
|
21764206 | for(index_t i = 0; i < length_; ++i) { |
| 1054 | 16709445 | x_[i] = -x_[i]; | |
| 1055 | } | ||
| 1056 | 5054761 | return *this; | |
| 1057 | } | ||
| 1058 | |||
| 1059 | /** | ||
| 1060 | * \brief Multiplies this expansion by a power of two. | ||
| 1061 | * \param[in] s the factor to be used to scale this expansion. | ||
| 1062 | * \return the new value of this expansion | ||
| 1063 | * \note Does not check for overflows/underflows | ||
| 1064 | * \pre \p s should be a (possibly negative) power of two. | ||
| 1065 | */ | ||
| 1066 | 2052283 | expansion& scale_fast(double s) { | |
| 1067 | // TODO: debug assert is_power_of_two(s) | ||
| 1068 |
2/2✓ Branch 0 taken 2377845 times.
✓ Branch 1 taken 2052283 times.
|
4430128 | for(index_t i = 0; i < length_; ++i) { |
| 1069 | 2377845 | x_[i] *= s; | |
| 1070 | } | ||
| 1071 | 2052283 | return *this; | |
| 1072 | } | ||
| 1073 | |||
| 1074 | /** | ||
| 1075 | * \brief Computes an approximation of the stored | ||
| 1076 | * value in this expansion. | ||
| 1077 | * \return an approximation of the stored value. | ||
| 1078 | */ | ||
| 1079 | 99996 | double estimate() const { | |
| 1080 | 99996 | double result = 0.0; | |
| 1081 |
2/2✓ Branch 1 taken 270800 times.
✓ Branch 2 taken 99996 times.
|
370796 | for(index_t i = 0; i < length(); ++i) { |
| 1082 | 270800 | result += x_[i]; | |
| 1083 | } | ||
| 1084 | 99996 | return result; | |
| 1085 | } | ||
| 1086 | |||
| 1087 | /** | ||
| 1088 | * \brief Gets the sign of the expansion. | ||
| 1089 | * \return the sign of this expansion, computed exactly. | ||
| 1090 | */ | ||
| 1091 | 55986096 | Sign sign() const { | |
| 1092 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 55986096 times.
|
55986096 | if(length() == 0) { |
| 1093 | ✗ | return ZERO; | |
| 1094 | } | ||
| 1095 | 55986096 | return geo_sgn(x_[length() - 1]); | |
| 1096 | } | ||
| 1097 | |||
| 1098 | /** | ||
| 1099 | * \brief Compares two expansions bit-by-bit | ||
| 1100 | * \details This function may return false even if the | ||
| 1101 | * expansion and \p rhs represent the same number | ||
| 1102 | * \retval true if the two expansions are the same | ||
| 1103 | * \retval false otherwise | ||
| 1104 | */ | ||
| 1105 | bool is_same_as(const expansion& rhs) const; | ||
| 1106 | |||
| 1107 | /** | ||
| 1108 | * \brief Compares an expansion and a double bit-by-bit | ||
| 1109 | * \details This function may return false even if the | ||
| 1110 | * expansion and \p rhs represent the same number | ||
| 1111 | * \retval true if the expansion has a single component | ||
| 1112 | * with value \p rhs | ||
| 1113 | * \retval false otherwise | ||
| 1114 | */ | ||
| 1115 | bool is_same_as(double rhs) const; | ||
| 1116 | |||
| 1117 | |||
| 1118 | /** | ||
| 1119 | * \brief Compares two expansions | ||
| 1120 | * \return the sign of this expansion minus rhs. | ||
| 1121 | */ | ||
| 1122 | Sign compare(const expansion& rhs) const; | ||
| 1123 | |||
| 1124 | /** | ||
| 1125 | * \brief Compares two expansions | ||
| 1126 | * \return the sign of this expansion minus rhs. | ||
| 1127 | */ | ||
| 1128 | Sign compare(double rhs) const; | ||
| 1129 | |||
| 1130 | /** | ||
| 1131 | * \brief Compares two expansions | ||
| 1132 | * \retval true if the two expansions represent the same | ||
| 1133 | * number | ||
| 1134 | * \retval false otherwise | ||
| 1135 | */ | ||
| 1136 | 2576435 | bool equals(const expansion& rhs) const { | |
| 1137 | 2576435 | return (compare(rhs) == ZERO); | |
| 1138 | } | ||
| 1139 | |||
| 1140 | /** | ||
| 1141 | * \brief Compares an expansion and a double | ||
| 1142 | * \retval true if the expansion and \p rhs represent | ||
| 1143 | * the same number | ||
| 1144 | * \retval false otherwise | ||
| 1145 | */ | ||
| 1146 | ✗ | bool equals(double rhs) const { | |
| 1147 | ✗ | return (compare(rhs) == ZERO); | |
| 1148 | } | ||
| 1149 | |||
| 1150 | /** | ||
| 1151 | * \brief Displays all the components of this expansion | ||
| 1152 | * (for debugging purposes). | ||
| 1153 | * \param[out] out an output stream used to print the components | ||
| 1154 | */ | ||
| 1155 | std::ostream& show(std::ostream& out) const { | ||
| 1156 | out << "expansion[" << length() << "] = ["; | ||
| 1157 | for(index_t i=0; i<length(); ++i) { | ||
| 1158 | out << (*this)[i] << " "; | ||
| 1159 | } | ||
| 1160 | out << "]"; | ||
| 1161 | return out; | ||
| 1162 | } | ||
| 1163 | |||
| 1164 | /** | ||
| 1165 | * \brief Gets a string representation of this expansion | ||
| 1166 | * \return a string with the length and components | ||
| 1167 | */ | ||
| 1168 | std::string to_string() const { | ||
| 1169 | std::ostringstream out; | ||
| 1170 | show(out); | ||
| 1171 | return out.str(); | ||
| 1172 | } | ||
| 1173 | |||
| 1174 | /** | ||
| 1175 | * \brief Optimizes the internal representation without changing the | ||
| 1176 | * represented value | ||
| 1177 | * \details this function can reduce the length of an expansion | ||
| 1178 | */ | ||
| 1179 | void optimize(); | ||
| 1180 | |||
| 1181 | /** | ||
| 1182 | * \brief Show global statistics | ||
| 1183 | */ | ||
| 1184 | static void show_all_stats(); | ||
| 1185 | |||
| 1186 | protected: | ||
| 1187 | /** | ||
| 1188 | * \brief Computes the required capacity of an expansion | ||
| 1189 | * to store an exact sub-product. | ||
| 1190 | * \param[in] a_length number of components in first sub-expansion | ||
| 1191 | * \param[in] b_length number of components in second sub-expansion | ||
| 1192 | * \return the required capacity of an expansion to store the | ||
| 1193 | * exact product of two expansions of lengths \p a_length | ||
| 1194 | * and \p b_length | ||
| 1195 | */ | ||
| 1196 | 132727 | static index_t sub_product_capacity( | |
| 1197 | index_t a_length, index_t b_length | ||
| 1198 | ) { | ||
| 1199 | 132727 | return a_length * b_length * 2; | |
| 1200 | } | ||
| 1201 | |||
| 1202 | /** | ||
| 1203 | * \brief Assigns a sub-product to this expansion. | ||
| 1204 | * \details Used by assign_product() when operating in balanced | ||
| 1205 | * distillation mode. Recursively assembles sub-sums. | ||
| 1206 | * \param[in] a a pointer to the first component of the first term | ||
| 1207 | * \param[in] a_length number of components in first term | ||
| 1208 | * \param[in] b second term | ||
| 1209 | * \return the new value of this expansion, with [a1...a_length]*b. | ||
| 1210 | */ | ||
| 1211 | expansion& assign_sub_product( | ||
| 1212 | const double* a, index_t a_length, const expansion& b | ||
| 1213 | ); | ||
| 1214 | |||
| 1215 | /** | ||
| 1216 | * \brief Expansion%s cannot be copied. | ||
| 1217 | */ | ||
| 1218 | expansion(const expansion& rhs) = delete; | ||
| 1219 | |||
| 1220 | /** | ||
| 1221 | * \brief Expansion%s cannot be copied. | ||
| 1222 | */ | ||
| 1223 | expansion& operator= (const expansion& rhs) = delete; | ||
| 1224 | |||
| 1225 | private: | ||
| 1226 | |||
| 1227 | /** | ||
| 1228 | * \brief Threshold in terms of expansion length for | ||
| 1229 | * allocating an expansion on the stack (if smaller) | ||
| 1230 | * or on the heap (if larger). | ||
| 1231 | */ | ||
| 1232 | #ifdef GEO_OS_APPLE | ||
| 1233 | static constexpr index_t MAX_CAPACITY_ON_STACK = 256; | ||
| 1234 | #else | ||
| 1235 | static constexpr index_t MAX_CAPACITY_ON_STACK = 1024; | ||
| 1236 | #endif | ||
| 1237 | index_t length_; | ||
| 1238 | index_t capacity_; | ||
| 1239 | double x_[2]; // x_ is in fact of size [capacity_] | ||
| 1240 | |||
| 1241 | friend class expansion_nt; | ||
| 1242 | }; | ||
| 1243 | |||
| 1244 | // =============== arithmetic operations =========================== | ||
| 1245 | |||
| 1246 | /** | ||
| 1247 | * \brief Creates an expansion from a double. | ||
| 1248 | * \param[in] a the double | ||
| 1249 | * \return a reference to an expansion, allocated on the | ||
| 1250 | * stack. | ||
| 1251 | * \code | ||
| 1252 | * const expansion& e1 = expansion_create(a); | ||
| 1253 | * \endcode | ||
| 1254 | * \warning Do not return or use the returned reference outside the | ||
| 1255 | * calling function. | ||
| 1256 | * \relates GEO::expansion | ||
| 1257 | */ | ||
| 1258 | #define expansion_create(a) \ | ||
| 1259 | new_expansion_on_stack(1)->assign(a) | ||
| 1260 | |||
| 1261 | |||
| 1262 | /** | ||
| 1263 | * \brief Creates an expansion from the absolute value of another | ||
| 1264 | * expansion | ||
| 1265 | * \param[in] e the expansion | ||
| 1266 | * \return a reference to an expansion, allocated on the | ||
| 1267 | * stack. | ||
| 1268 | * \code | ||
| 1269 | * const expansion& e1 = expansion_abs(e); | ||
| 1270 | * \endcode | ||
| 1271 | * \warning Do not return or use the returned reference outside the | ||
| 1272 | * calling function. | ||
| 1273 | * \relates GEO::expansion | ||
| 1274 | */ | ||
| 1275 | #define expansion_abs(e) \ | ||
| 1276 | new_expansion_on_stack(e.length())->assign_abs(e) | ||
| 1277 | |||
| 1278 | /** | ||
| 1279 | * \brief Computes an expansion that represents the exact | ||
| 1280 | * sum of its arguments. | ||
| 1281 | * \param[in] a a double or an expansion | ||
| 1282 | * \param[in] b a double or an expansion | ||
| 1283 | * \return a reference to an expansion, allocated on the stack. | ||
| 1284 | * \code | ||
| 1285 | * expansion& e1 = ...; | ||
| 1286 | * expansion& e2 = ...; | ||
| 1287 | * expansion& e3 = expansion_sum(e1,e2); | ||
| 1288 | * double x = ...; | ||
| 1289 | * expansion& e4 = expansion_sum(e1,x); | ||
| 1290 | * \endcode | ||
| 1291 | * \warning Do not return or use the returned reference outside the | ||
| 1292 | * calling function. | ||
| 1293 | * \relates GEO::expansion | ||
| 1294 | */ | ||
| 1295 | #define expansion_sum(a, b) \ | ||
| 1296 | new_expansion_on_stack( \ | ||
| 1297 | expansion::sum_capacity(a, b) \ | ||
| 1298 | )->assign_sum(a, b) | ||
| 1299 | |||
| 1300 | /** | ||
| 1301 | * \brief Computes an expansion that represents the exact | ||
| 1302 | * sum of its arguments. | ||
| 1303 | * \param[in] a an expansion | ||
| 1304 | * \param[in] b an expansion | ||
| 1305 | * \param[in] c an expansion | ||
| 1306 | * \return a reference to an expansion, allocated on the stack. | ||
| 1307 | * \code | ||
| 1308 | * expansion& e1 = ...; | ||
| 1309 | * expansion& e2 = ...; | ||
| 1310 | * expansion& e3 = ...; | ||
| 1311 | * expansion& e4 = expansion_sum3(e1,e2,e3); | ||
| 1312 | * \endcode | ||
| 1313 | * \warning Do not return or use the returned reference outside the | ||
| 1314 | * calling function. | ||
| 1315 | * \relates GEO::expansion | ||
| 1316 | */ | ||
| 1317 | #define expansion_sum3(a, b, c) \ | ||
| 1318 | new_expansion_on_stack( \ | ||
| 1319 | expansion::sum_capacity(a, b, c) \ | ||
| 1320 | )->assign_sum(a, b, c) | ||
| 1321 | |||
| 1322 | /** | ||
| 1323 | * \brief Computes an expansion that represents the exact | ||
| 1324 | * sum of its arguments. | ||
| 1325 | * \param[in] a an expansion | ||
| 1326 | * \param[in] b an expansion | ||
| 1327 | * \param[in] c an expansion | ||
| 1328 | * \param[in] d an expansion | ||
| 1329 | * \return a reference to an expansion, allocated on the stack. | ||
| 1330 | * \code | ||
| 1331 | * expansion& e1 = ...; | ||
| 1332 | * expansion& e2 = ...; | ||
| 1333 | * expansion& e3 = ...; | ||
| 1334 | * expansion& e4 = ...; | ||
| 1335 | * expansion& e5 = expansion_sum4(e1,e2,e3,e4); | ||
| 1336 | * \endcode | ||
| 1337 | * \warning Do not return or use the returned reference outside the | ||
| 1338 | * calling function. | ||
| 1339 | * \relates GEO::expansion | ||
| 1340 | */ | ||
| 1341 | |||
| 1342 | #define expansion_sum4(a, b, c, d) \ | ||
| 1343 | new_expansion_on_stack( \ | ||
| 1344 | expansion::sum_capacity(a, b, c, d) \ | ||
| 1345 | )->assign_sum(a, b, c, d) | ||
| 1346 | |||
| 1347 | /** | ||
| 1348 | * \brief Computes an expansion that represents the exact | ||
| 1349 | * difference of its arguments. | ||
| 1350 | * \param[in] a a double or an expansion | ||
| 1351 | * \param[in] b a double or an expansion | ||
| 1352 | * \return a reference to an expansion, allocated on the stack. | ||
| 1353 | * \code | ||
| 1354 | * expansion& e1 = ...; | ||
| 1355 | * expansion& e2 = ...; | ||
| 1356 | * expansion& e3 = expansion_diff(e1,e2); | ||
| 1357 | * double x = ...; | ||
| 1358 | * expansion& e4 = expansion_diff(e1,x); | ||
| 1359 | * \endcode | ||
| 1360 | * \warning Do not return or use the returned reference outside the | ||
| 1361 | * calling function. | ||
| 1362 | * \relates GEO::expansion | ||
| 1363 | */ | ||
| 1364 | #define expansion_diff(a, b) \ | ||
| 1365 | new_expansion_on_stack( \ | ||
| 1366 | expansion::diff_capacity(a, b) \ | ||
| 1367 | )->assign_diff(a, b) | ||
| 1368 | |||
| 1369 | /** | ||
| 1370 | * \brief Computes an expansion that represents the exact | ||
| 1371 | * product of its arguments. | ||
| 1372 | * \param[in] a a double or an expansion | ||
| 1373 | * \param[in] b a double or an expansion | ||
| 1374 | * \return a reference to an expansion, allocated on the stack. | ||
| 1375 | * \code | ||
| 1376 | * expansion& e1 = ...; | ||
| 1377 | * expansion& e2 = ...; | ||
| 1378 | * expansion& e3 = expansion_product(e1,e2); | ||
| 1379 | * double x = ...; | ||
| 1380 | * expansion& e4 = expansion_product(e1,x); | ||
| 1381 | * \endcode | ||
| 1382 | * \warning Do not return or use the returned reference outside the | ||
| 1383 | * calling function. | ||
| 1384 | * \relates GEO::expansion | ||
| 1385 | */ | ||
| 1386 | #define expansion_product(a, b) \ | ||
| 1387 | new_expansion_on_stack( \ | ||
| 1388 | expansion::product_capacity(a, b) \ | ||
| 1389 | )->assign_product(a, b) | ||
| 1390 | |||
| 1391 | /** | ||
| 1392 | * \brief Computes an expansion that represents the exact | ||
| 1393 | * product of its arguments. | ||
| 1394 | * \param[in] a an expansion | ||
| 1395 | * \param[in] b an expansion | ||
| 1396 | * \param[in] c an expansion | ||
| 1397 | * \return a reference to an expansion, allocated on the stack. | ||
| 1398 | * \code | ||
| 1399 | * expansion& e1 = ...; | ||
| 1400 | * expansion& e2 = ...; | ||
| 1401 | * expansion& e3 = ...; | ||
| 1402 | * expansion& e4 = expansion_product3(e1,e2,e3); | ||
| 1403 | * \endcode | ||
| 1404 | * \warning Do not return or use the returned reference outside the | ||
| 1405 | * calling function. | ||
| 1406 | * \relates GEO::expansion | ||
| 1407 | */ | ||
| 1408 | #define expansion_product3(a, b, c) \ | ||
| 1409 | new_expansion_on_stack( \ | ||
| 1410 | expansion::product_capacity(a, b, c) \ | ||
| 1411 | )->assign_product(a, b, c) | ||
| 1412 | |||
| 1413 | /** | ||
| 1414 | * \brief Computes an expansion that represents the exact | ||
| 1415 | * square of its argument. | ||
| 1416 | * \param[in] a a double or an expansion | ||
| 1417 | * \return a reference to an expansion, allocated on the stack. | ||
| 1418 | * \code | ||
| 1419 | * expansion& e1 = ...; | ||
| 1420 | * expansion& e2 = expansion_square(e1); | ||
| 1421 | * double x = ...; | ||
| 1422 | * expansion& e3 = expansion_square(x); | ||
| 1423 | * \endcode | ||
| 1424 | * \warning Do not return or use the returned reference outside the | ||
| 1425 | * calling function. | ||
| 1426 | * \relates GEO::expansion | ||
| 1427 | */ | ||
| 1428 | #define expansion_square(a) \ | ||
| 1429 | new_expansion_on_stack( \ | ||
| 1430 | expansion::square_capacity(a) \ | ||
| 1431 | )->assign_square(a) | ||
| 1432 | |||
| 1433 | // =============== determinants ===================================== | ||
| 1434 | |||
| 1435 | /** | ||
| 1436 | * \brief Computes an expansion that represents the exact | ||
| 1437 | * 2x2 determinant of its arguments. | ||
| 1438 | * \return a reference to an expansion, allocated on the stack. | ||
| 1439 | * \code | ||
| 1440 | * const expansion& a11 = ...; | ||
| 1441 | * const expansion& a12 = ...; | ||
| 1442 | * const expansion& a21 = ...; | ||
| 1443 | * const expansion& a22 = ...; | ||
| 1444 | * expansion& d12 = expansion_set2x2(a11,a12,a21,a22); | ||
| 1445 | * \endcode | ||
| 1446 | * \relates GEO::expansion | ||
| 1447 | */ | ||
| 1448 | #define expansion_det2x2(a11, a12, a21, a22) \ | ||
| 1449 | new_expansion_on_stack( \ | ||
| 1450 | expansion::det2x2_capacity(a11, a12, a21, a22) \ | ||
| 1451 | )->assign_det2x2(a11, a12, a21, a22) | ||
| 1452 | |||
| 1453 | /** | ||
| 1454 | * \brief Computes an expansion that represents the exact | ||
| 1455 | * 3x3 determinant of its arguments. | ||
| 1456 | * \return a reference to an expansion, allocated on the stack. | ||
| 1457 | * \code | ||
| 1458 | * const expansion& a11 = ...; | ||
| 1459 | * ... | ||
| 1460 | * const expansion& a33 = ...; | ||
| 1461 | * expansion& d = expansion_det3x3( | ||
| 1462 | * a11,a12,a13,a21,a22,a23,a31,a32,a33 | ||
| 1463 | * ); | ||
| 1464 | * \endcode | ||
| 1465 | * \warning Do not return or use the returned reference outside the | ||
| 1466 | * calling function. | ||
| 1467 | * \relates GEO::expansion | ||
| 1468 | */ | ||
| 1469 | #define expansion_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33) \ | ||
| 1470 | new_expansion_on_stack( \ | ||
| 1471 | expansion::det3x3_capacity(a11,a12,a13,a21,a22,a23,a31,a32,a33) \ | ||
| 1472 | )->assign_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33) | ||
| 1473 | |||
| 1474 | /** | ||
| 1475 | * \brief Computes an expansion that represents the exact | ||
| 1476 | * 3x3 determinant of its arguments where the first row | ||
| 1477 | * is 1 1 1. | ||
| 1478 | * \return a reference to an expansion, allocated on the stack. | ||
| 1479 | * \code | ||
| 1480 | * const expansion& a21 = ...; | ||
| 1481 | * ... | ||
| 1482 | * const expansion& a33 = ...; | ||
| 1483 | * expansion& d = expansion_det_111_2x3( | ||
| 1484 | * a21,a22,a23,a31,a32,a33 | ||
| 1485 | * ); | ||
| 1486 | * \endcode | ||
| 1487 | * \warning Do not return or use the returned reference outside the | ||
| 1488 | * calling function. | ||
| 1489 | * \relates GEO::expansion | ||
| 1490 | */ | ||
| 1491 | #define expansion_det_111_2x3(a21, a22, a23, a31, a32, a33) \ | ||
| 1492 | new_expansion_on_stack( \ | ||
| 1493 | expansion::det_111_2x3_capacity(a21, a22, a23, a31, a32, a33) \ | ||
| 1494 | )->assign_det_111_2x3(a21, a22, a23, a31, a32, a33) | ||
| 1495 | |||
| 1496 | // =============== geometric functions ============================== | ||
| 1497 | |||
| 1498 | /** | ||
| 1499 | * \brief Computes an expansion that represents the exact | ||
| 1500 | * squared distance between its argument. | ||
| 1501 | * \param[in] a first point (specified as a const double*) | ||
| 1502 | * \param[in] b second point (specified as a const double*) | ||
| 1503 | * \param[in] dim dimension of the points | ||
| 1504 | * \return a reference to an expansion, allocated on the stack. | ||
| 1505 | * \code | ||
| 1506 | * const double* p1 = ...; | ||
| 1507 | * const double* p2 = ...; | ||
| 1508 | * expansion& d12 = expansion_sq_dist(p1,p2,3); | ||
| 1509 | * \endcode | ||
| 1510 | * \warning Do not return or use the returned reference outside the | ||
| 1511 | * calling function. | ||
| 1512 | * \relates GEO::expansion | ||
| 1513 | */ | ||
| 1514 | #define expansion_sq_dist(a, b, dim) \ | ||
| 1515 | new_expansion_on_stack( \ | ||
| 1516 | expansion::sq_dist_capacity(dim) \ | ||
| 1517 | )->assign_sq_dist(a, b, dim) | ||
| 1518 | |||
| 1519 | /** | ||
| 1520 | * \brief Computes an expansion that represents the exact | ||
| 1521 | * dot product dot(a-c,b-c) | ||
| 1522 | * \param[in] a first point (specified as a const double*) | ||
| 1523 | * \param[in] b second point (specified as a const double*) | ||
| 1524 | * \param[in] c third point (specified as a const double*) | ||
| 1525 | * \param[in] dim dimension of the points | ||
| 1526 | * \return a reference to an expansion, allocated on the stack. | ||
| 1527 | * \code | ||
| 1528 | * const double* p1 = ...; | ||
| 1529 | * const double* p2 = ...; | ||
| 1530 | * const double* p0 = ...; | ||
| 1531 | * expansion& dot12 = expansion_dot_at(p1,p2,p0,3); | ||
| 1532 | * \endcode | ||
| 1533 | * \warning Do not return or use the returned reference outside the | ||
| 1534 | * calling function. | ||
| 1535 | * \relates GEO::expansion | ||
| 1536 | */ | ||
| 1537 | #define expansion_dot_at(a, b, c, dim) \ | ||
| 1538 | new_expansion_on_stack( \ | ||
| 1539 | expansion::dot_at_capacity(dim) \ | ||
| 1540 | )->assign_dot_at(a, b, c, dim) | ||
| 1541 | |||
| 1542 | |||
| 1543 | /** | ||
| 1544 | * \brief Computes an expansion that represents the exact | ||
| 1545 | * squared length of a 3d vector | ||
| 1546 | * \param[in] x,y,z coordinates of the vector (specified as expansion) | ||
| 1547 | * \return a reference to an expansion, allocated on the stack. | ||
| 1548 | * \code | ||
| 1549 | * const expansion& x = ...; | ||
| 1550 | * const expansion& y = ...; | ||
| 1551 | * const expansion& z = ...; | ||
| 1552 | * expansion& l = expansion_length2(x,y,z); | ||
| 1553 | * \endcode | ||
| 1554 | * \warning Do not return or use the returned reference outside the | ||
| 1555 | * calling function. | ||
| 1556 | * \relates GEO::expansion | ||
| 1557 | */ | ||
| 1558 | #define expansion_length2(x,y,z) \ | ||
| 1559 | new_expansion_on_stack( \ | ||
| 1560 | expansion::length2_capacity(x,y,z) \ | ||
| 1561 | )->assign_length2(x,y,z) | ||
| 1562 | |||
| 1563 | /************************************************************************/ | ||
| 1564 | |||
| 1565 | /** | ||
| 1566 | * \brief Computes the sign of a 2x2 determinant | ||
| 1567 | * \details Specialization using the low-evel API for expansions. | ||
| 1568 | * This gains some performance as compared to using CGAL's | ||
| 1569 | * determinant template with expansion_nt. | ||
| 1570 | */ | ||
| 1571 | Sign GEOGRAM_API sign_of_expansion_determinant( | ||
| 1572 | const expansion& a00,const expansion& a01, | ||
| 1573 | const expansion& a10,const expansion& a11 | ||
| 1574 | ); | ||
| 1575 | |||
| 1576 | /** | ||
| 1577 | * \brief Computes the sign of a 3x3 determinant | ||
| 1578 | * \details Specialization using the low-evel API for expansions. | ||
| 1579 | * This gains some performance as compared to using CGAL's determinant | ||
| 1580 | * template with expansion_nt. | ||
| 1581 | */ | ||
| 1582 | Sign GEOGRAM_API sign_of_expansion_determinant( | ||
| 1583 | const expansion& a00,const expansion& a01,const expansion& a02, | ||
| 1584 | const expansion& a10,const expansion& a11,const expansion& a12, | ||
| 1585 | const expansion& a20,const expansion& a21,const expansion& a22 | ||
| 1586 | ); | ||
| 1587 | |||
| 1588 | /** | ||
| 1589 | * \brief Computes the sign of a 4x4 determinant | ||
| 1590 | * \details Specialization using the low-evel API for expansions. | ||
| 1591 | * This gains some performance as compared to using CGAL's determinant | ||
| 1592 | * template with expansion_nt. | ||
| 1593 | */ | ||
| 1594 | Sign GEOGRAM_API sign_of_expansion_determinant( | ||
| 1595 | const expansion& a00,const expansion& a01, | ||
| 1596 | const expansion& a02,const expansion& a03, | ||
| 1597 | const expansion& a10,const expansion& a11, | ||
| 1598 | const expansion& a12,const expansion& a13, | ||
| 1599 | const expansion& a20,const expansion& a21, | ||
| 1600 | const expansion& a22,const expansion& a23, | ||
| 1601 | const expansion& a30,const expansion& a31, | ||
| 1602 | const expansion& a32,const expansion& a33 | ||
| 1603 | ); | ||
| 1604 | |||
| 1605 | /************************************************************************/ | ||
| 1606 | |||
| 1607 | /** | ||
| 1608 | * \brief Adds a scalar to an expansion, eliminating zero components | ||
| 1609 | * from the output expansion. | ||
| 1610 | * \param[in] e first expansion | ||
| 1611 | * \param[in] b double to be added to \p e | ||
| 1612 | * \param[out] h the result \p e + \p b | ||
| 1613 | * \details Sets \p h = (\p e + \p b). \p e and \p h can be the same. | ||
| 1614 | * This function is adapted from Jonathan Shewchuk's code. | ||
| 1615 | * See the long version of his paper for details. | ||
| 1616 | * Maintains the nonoverlapping property. If round-to-even is used (as | ||
| 1617 | * with IEEE 754), maintains the strongly nonoverlapping and nonadjacent | ||
| 1618 | * properties as well. (That is, if e has one of these properties, so | ||
| 1619 | * will h.) | ||
| 1620 | */ | ||
| 1621 | void GEOGRAM_API grow_expansion_zeroelim( | ||
| 1622 | const expansion& e, double b, expansion& h | ||
| 1623 | ); | ||
| 1624 | |||
| 1625 | /** | ||
| 1626 | * \brief Multiplies an expansion by a scalar, | ||
| 1627 | * eliminating zero components from the | ||
| 1628 | * output expansion. | ||
| 1629 | * \param[in] e an expansion | ||
| 1630 | * \param[in] b the double to be multiplied by \p e | ||
| 1631 | * \param[out] h the result \p b * \p e | ||
| 1632 | * \details (sets \p h = \p b * \p e). \p e and \p h cannot be the same. | ||
| 1633 | * This function is adapted from Jonathan Shewchuk's code. | ||
| 1634 | * See either version of his paper for details. | ||
| 1635 | * Maintains the nonoverlapping property. If round-to-even is used (as | ||
| 1636 | * with IEEE 754), maintains the strongly nonoverlapping and nonadjacent | ||
| 1637 | * properties as well. (That is, if e has one of these properties, so | ||
| 1638 | * will h.) | ||
| 1639 | */ | ||
| 1640 | void GEOGRAM_API scale_expansion_zeroelim( | ||
| 1641 | const expansion& e, double b, expansion& h | ||
| 1642 | ); | ||
| 1643 | |||
| 1644 | /** | ||
| 1645 | * \brief Sums two expansions, eliminating zero | ||
| 1646 | * components from the output expansion (sets \p h = \p e + \p f). | ||
| 1647 | * \param[in] e the first expansion | ||
| 1648 | * \param[in] f the second expansion | ||
| 1649 | * \param[out] h the result \p e + \p f | ||
| 1650 | * \details h cannot be e or f. | ||
| 1651 | * This function is adapted from Jonathan Shewchuk's code. | ||
| 1652 | * See the long version of his paper for details. | ||
| 1653 | * If round-to-even is used (as with IEEE 754), maintains the strongly | ||
| 1654 | * nonoverlapping property. (That is, if e is strongly nonoverlapping, h | ||
| 1655 | * will be also.) Does NOT maintain the nonoverlapping or nonadjacent | ||
| 1656 | * properties. | ||
| 1657 | * | ||
| 1658 | */ | ||
| 1659 | void GEOGRAM_API fast_expansion_sum_zeroelim( | ||
| 1660 | const expansion& e, const expansion& f, expansion& h | ||
| 1661 | ); | ||
| 1662 | |||
| 1663 | |||
| 1664 | /** | ||
| 1665 | * \brief Computes the difference of two expansions, eliminating zero | ||
| 1666 | * components from the output expansion | ||
| 1667 | * \param[in] e first expansion | ||
| 1668 | * \param[in] f second expansion to be subtracted from e | ||
| 1669 | * \param[out] h the result \p e - \p f | ||
| 1670 | * \details Sets \p h = (\p e - \p f). \p h cannot be \p e or \p f. | ||
| 1671 | * This function is adapted from Jonathan Shewchuk's code. | ||
| 1672 | * See the long version of his paper for details. | ||
| 1673 | * If round-to-even is used (as with IEEE 754), maintains the strongly | ||
| 1674 | * nonoverlapping property. (That is, if e is strongly nonoverlapping, h | ||
| 1675 | * will be also.) Does NOT maintain the nonoverlapping or nonadjacent | ||
| 1676 | * properties. | ||
| 1677 | */ | ||
| 1678 | void GEOGRAM_API fast_expansion_diff_zeroelim( | ||
| 1679 | const expansion& e, const expansion& f, expansion& h | ||
| 1680 | ); | ||
| 1681 | |||
| 1682 | /************************************************************************/ | ||
| 1683 | } | ||
| 1684 | |||
| 1685 | #endif | ||
| 1686 |