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