| 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 GEO_INTERVAL_NT | ||
| 41 | #define GEO_INTERVAL_NT | ||
| 42 | |||
| 43 | #include <geogram/numerics/expansion_nt.h> | ||
| 44 | #include <geogram/basic/vechg.h> | ||
| 45 | #include <iomanip> | ||
| 46 | #include <limits> | ||
| 47 | #include <cmath> | ||
| 48 | |||
| 49 | #ifdef __SSE__ | ||
| 50 | #include <xmmintrin.h> | ||
| 51 | #else | ||
| 52 | #include <fenv.h> | ||
| 53 | #endif | ||
| 54 | |||
| 55 | // Uncomment to activate checks (keeps an arbitrary precision | ||
| 56 | // representation of the number and checks that the interval | ||
| 57 | // contains it). | ||
| 58 | // #define INTERVAL_CHECK | ||
| 59 | |||
| 60 | // Uncomment to deactivate all interval filters in geogram | ||
| 61 | // (used for debugging and performance tests) | ||
| 62 | // #define GEO_NO_INTERVALS | ||
| 63 | |||
| 64 | namespace GEO { | ||
| 65 | |||
| 66 | /** | ||
| 67 | * \brief Base class for interval arithmetics. | ||
| 68 | * \details Has the verification mechanics | ||
| 69 | * (define INTERVAL_CHECK to activate it). Verification | ||
| 70 | * stores an additional expansion_nt and checks | ||
| 71 | * that this expansion_nt is contained by the | ||
| 72 | * interval. | ||
| 73 | */ | ||
| 74 | class intervalBase { | ||
| 75 | public: | ||
| 76 | |||
| 77 | static void set_FPU_round_to_nearest() { | ||
| 78 | #ifdef __SSE__ | ||
| 79 | _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST); | ||
| 80 | #else | ||
| 81 | fesetround(FE_TONEAREST); | ||
| 82 | #endif | ||
| 83 | } | ||
| 84 | |||
| 85 | static void set_FPU_round_to_upper() { | ||
| 86 | #ifdef __SSE__ | ||
| 87 | _MM_SET_ROUNDING_MODE(_MM_ROUND_UP); | ||
| 88 | #else | ||
| 89 | fesetround(FE_UPWARD); | ||
| 90 | #endif | ||
| 91 | } | ||
| 92 | |||
| 93 | enum Sign2 { | ||
| 94 | SIGN2_ERROR = -1, | ||
| 95 | SIGN2_ZERO = 0, | ||
| 96 | SIGN2_NP, | ||
| 97 | SIGN2_PP, | ||
| 98 | SIGN2_ZP, | ||
| 99 | SIGN2_NN, | ||
| 100 | SIGN2_NZ | ||
| 101 | }; | ||
| 102 | |||
| 103 | 12600724 | static bool sign_is_determined(Sign2 s) { | |
| 104 | return | ||
| 105 |
2/2✓ Branch 0 taken 7234575 times.
✓ Branch 1 taken 5366149 times.
|
12600724 | s == SIGN2_ZERO || |
| 106 |
3/4✓ Branch 0 taken 12600724 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5746559 times.
✓ Branch 3 taken 1488016 times.
|
25201448 | s == SIGN2_NN || |
| 107 | 12600724 | s == SIGN2_PP ; | |
| 108 | } | ||
| 109 | |||
| 110 | 1467220 | static bool sign_is_non_zero(Sign2 s) { | |
| 111 | return | ||
| 112 |
4/4✓ Branch 0 taken 881506 times.
✓ Branch 1 taken 585714 times.
✓ Branch 2 taken 837723 times.
✓ Branch 3 taken 43783 times.
|
1467220 | s == SIGN2_NN || |
| 113 | 1467220 | s == SIGN2_PP ; | |
| 114 | } | ||
| 115 | |||
| 116 | 6131369 | static Sign convert_sign(Sign2 s) { | |
| 117 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 6131369 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
6131369 | geo_assert(sign_is_determined(s)); |
| 118 |
2/2✓ Branch 0 taken 2923871 times.
✓ Branch 1 taken 3207498 times.
|
6131369 | if(s == SIGN2_NN) { |
| 119 | 2923871 | return NEGATIVE; | |
| 120 | } | ||
| 121 |
1/2✓ Branch 0 taken 3207498 times.
✗ Branch 1 not taken.
|
3207498 | if(s == SIGN2_PP) { |
| 122 | 3207498 | return POSITIVE; | |
| 123 | } | ||
| 124 | ✗ | return ZERO; | |
| 125 | } | ||
| 126 | |||
| 127 | 49409314 | intervalBase() { | |
| 128 | 49409314 | control_set(0); | |
| 129 | 49409314 | } | |
| 130 | |||
| 131 | 7891308 | intervalBase(double x) { | |
| 132 | 7891308 | control_set(x); | |
| 133 | 7891308 | } | |
| 134 | |||
| 135 | intervalBase(const intervalBase& rhs) = default; | ||
| 136 | |||
| 137 | intervalBase& operator=(const intervalBase& rhs) = default; | ||
| 138 | |||
| 139 | protected: | ||
| 140 | #ifdef INTERVAL_CHECK | ||
| 141 | void control_set(const expansion_nt& x) { | ||
| 142 | control_ = x; | ||
| 143 | } | ||
| 144 | void control_set(const intervalBase& x) { | ||
| 145 | control_ = x.control_; | ||
| 146 | } | ||
| 147 | void control_set(double x) { | ||
| 148 | control_ = expansion_nt(x); | ||
| 149 | } | ||
| 150 | void control_negate() { | ||
| 151 | control_.rep().negate(); | ||
| 152 | } | ||
| 153 | void control_add(const intervalBase& x) { | ||
| 154 | control_ += x.control_; | ||
| 155 | } | ||
| 156 | void control_sub(const intervalBase& x) { | ||
| 157 | control_ -= x.control_; | ||
| 158 | } | ||
| 159 | void control_mul(const intervalBase& x) { | ||
| 160 | control_ *= x.control_; | ||
| 161 | } | ||
| 162 | void control_check(double inf, double sup) { | ||
| 163 | typedef std::numeric_limits< double > dbl; | ||
| 164 | if(inf > sup) { | ||
| 165 | std::cerr.precision(dbl::max_digits10); | ||
| 166 | std::cerr << "inf() > sup() !!" << std::endl; | ||
| 167 | std::cerr << "inf()=" << inf << std::endl; | ||
| 168 | std::cerr << "sup()=" << sup << std::endl; | ||
| 169 | geo_assert_not_reached; | ||
| 170 | } | ||
| 171 | if(control_ < inf || control_ > sup) { | ||
| 172 | std::cerr.precision(dbl::max_digits10); | ||
| 173 | std::cerr << "[" << inf << "," << sup << "]" | ||
| 174 | << " " << control_.estimate() << ":" | ||
| 175 | << control_.rep().length() | ||
| 176 | << std::endl; | ||
| 177 | geo_assert_not_reached; | ||
| 178 | } | ||
| 179 | } | ||
| 180 | expansion_nt control_; /**< exact represented value, for tests */ | ||
| 181 | #else | ||
| 182 | 57300622 | void control_set(double x) { | |
| 183 | 57300622 | geo_argused(x); | |
| 184 | 57300622 | } | |
| 185 | 49409314 | void control_set(const expansion_nt& x) { | |
| 186 | 49409314 | geo_argused(x); | |
| 187 | 49409314 | } | |
| 188 | void control_set(const intervalBase& x) { | ||
| 189 | geo_argused(x); | ||
| 190 | } | ||
| 191 | 212819 | void control_negate() { | |
| 192 | 212819 | } | |
| 193 | 5357103 | void control_add(const intervalBase& x) { | |
| 194 | 5357103 | geo_argused(x); | |
| 195 | 5357103 | } | |
| 196 | 31640458 | void control_sub(const intervalBase& x) { | |
| 197 | 31640458 | geo_argused(x); | |
| 198 | 31640458 | } | |
| 199 | 59233686 | void control_mul(const intervalBase& x) { | |
| 200 | 59233686 | geo_argused(x); | |
| 201 | 59233686 | } | |
| 202 | void control_check(double inf, double sup) { | ||
| 203 | geo_argused(inf); | ||
| 204 | geo_argused(sup); | ||
| 205 | } | ||
| 206 | #endif | ||
| 207 | }; | ||
| 208 | |||
| 209 | |||
| 210 | /*******************************************************************/ | ||
| 211 | |||
| 212 | /** | ||
| 213 | * \brief Interval arithmetics in round to upper (RU) mode. | ||
| 214 | */ | ||
| 215 | class intervalRU : public intervalBase { | ||
| 216 | public: | ||
| 217 | |||
| 218 | /** | ||
| 219 | * \brief Sets FPU rounding mode for using this type of interval. | ||
| 220 | * \details Declare one instance at the beginning of a function or | ||
| 221 | * bloc of code that uses intervalRU. Do not mix other floating | ||
| 222 | * point operations in the bloc. | ||
| 223 | * \code | ||
| 224 | * { | ||
| 225 | * intervalRU::Rounding rounding; | ||
| 226 | * ... do some computations with intervalRU | ||
| 227 | * } | ||
| 228 | * \endcode | ||
| 229 | */ | ||
| 230 | struct Rounding { | ||
| 231 | Rounding() { | ||
| 232 | set_FPU_round_to_upper(); | ||
| 233 | } | ||
| 234 | ~Rounding() { | ||
| 235 | set_FPU_round_to_nearest(); | ||
| 236 | } | ||
| 237 | }; | ||
| 238 | |||
| 239 | intervalRU() : | ||
| 240 | intervalBase(), | ||
| 241 | ln_(0.0), | ||
| 242 | u_(0.0) | ||
| 243 | { | ||
| 244 | control_check(); | ||
| 245 | } | ||
| 246 | |||
| 247 | intervalRU(double x) : | ||
| 248 | intervalBase(x), | ||
| 249 | ln_(-x), | ||
| 250 | u_(x) | ||
| 251 | { | ||
| 252 | control_check(); | ||
| 253 | } | ||
| 254 | |||
| 255 | intervalRU(double l, double u) : ln_(-l), u_(u) { | ||
| 256 | // note: we cannot control check here. | ||
| 257 | } | ||
| 258 | |||
| 259 | intervalRU(const intervalRU& rhs) = default; | ||
| 260 | |||
| 261 | intervalRU(const expansion_nt& rhs) { | ||
| 262 | *this = rhs; | ||
| 263 | } | ||
| 264 | |||
| 265 | intervalRU& operator=(const intervalRU& rhs) = default; | ||
| 266 | |||
| 267 | intervalRU& operator=(double rhs) { | ||
| 268 | ln_ = -rhs; | ||
| 269 | u_ = rhs; | ||
| 270 | control_set(rhs); | ||
| 271 | control_check(); | ||
| 272 | return *this; | ||
| 273 | } | ||
| 274 | |||
| 275 | intervalRU& operator=(const expansion_nt& rhs) { | ||
| 276 | |||
| 277 | // Optimized expansion-to-interval conversion: | ||
| 278 | // | ||
| 279 | // Add components starting from the one of largest magnitude | ||
| 280 | // Stop as soon as next component is smaller than ulp (and then | ||
| 281 | // expand interval by ulp). | ||
| 282 | |||
| 283 | index_t l = rhs.length(); | ||
| 284 | ln_ = -rhs.component(l-1); | ||
| 285 | u_ = rhs.component(l-1); | ||
| 286 | |||
| 287 | for(int comp_idx=int(l)-2; comp_idx>=0; --comp_idx) { | ||
| 288 | double comp = rhs.component(index_t(comp_idx)); | ||
| 289 | u_ += comp; | ||
| 290 | ln_ -= comp; | ||
| 291 | if(comp > 0) { | ||
| 292 | double new_u = u_ + comp; | ||
| 293 | if(new_u == u_) { | ||
| 294 | u_ = std::nextafter( | ||
| 295 | u_, std::numeric_limits<double>::infinity() | ||
| 296 | ); | ||
| 297 | break; | ||
| 298 | } else { | ||
| 299 | u_ = new_u; | ||
| 300 | } | ||
| 301 | } else { | ||
| 302 | // If we stored l, we would write: | ||
| 303 | // new_l = l + comp | ||
| 304 | // But we store ln = -l, so we write: | ||
| 305 | // -new_ln = -ln + comp | ||
| 306 | // Which means: | ||
| 307 | // new_ln = ln - comp | ||
| 308 | |||
| 309 | double new_ln = ln_ - comp; | ||
| 310 | if(new_ln == ln_) { | ||
| 311 | ln_ = std::nextafter( | ||
| 312 | ln_, std::numeric_limits<double>::infinity() | ||
| 313 | ); | ||
| 314 | break; | ||
| 315 | } else { | ||
| 316 | ln_ = new_ln; | ||
| 317 | } | ||
| 318 | } | ||
| 319 | } | ||
| 320 | |||
| 321 | control_set(rhs); | ||
| 322 | control_check(); | ||
| 323 | return *this; | ||
| 324 | } | ||
| 325 | |||
| 326 | double inf() const { | ||
| 327 | return -ln_; | ||
| 328 | } | ||
| 329 | |||
| 330 | double sup() const { | ||
| 331 | return u_; | ||
| 332 | } | ||
| 333 | |||
| 334 | double estimate() const { | ||
| 335 | // 0.5*(lb+ub) -> | ||
| 336 | return 0.5*(-ln_+u_); | ||
| 337 | } | ||
| 338 | |||
| 339 | bool is_nan() const { | ||
| 340 | return !(ln_==ln_) || !(u_==u_); | ||
| 341 | } | ||
| 342 | |||
| 343 | Sign2 sign() const { | ||
| 344 | // Branchless | ||
| 345 | int lz = int(ln_ == 0); | ||
| 346 | int ln = int(ln_ > 0); // inverted, it is ln_ !!! | ||
| 347 | int lp = int(ln_ < 0); // inverted, it is ln_ !!! | ||
| 348 | int uz = int(u_ == 0); | ||
| 349 | int un = int(u_ < 0); | ||
| 350 | int up = int(u_ > 0); | ||
| 351 | Sign2 result = Sign2( | ||
| 352 | ln*up*SIGN2_NP+ | ||
| 353 | lp*up*SIGN2_PP+ | ||
| 354 | lz*up*SIGN2_ZP+ | ||
| 355 | ln*un*SIGN2_NN+ | ||
| 356 | ln*uz*SIGN2_NZ | ||
| 357 | ); | ||
| 358 | result = Sign2( | ||
| 359 | int(result) + | ||
| 360 | int(result==SIGN2_ZERO && !(lz&&uz)) * SIGN2_ERROR | ||
| 361 | ); | ||
| 362 | return result; | ||
| 363 | } | ||
| 364 | |||
| 365 | intervalRU& negate() { | ||
| 366 | std::swap(ln_, u_); | ||
| 367 | control_negate(); | ||
| 368 | control_check(); | ||
| 369 | return *this; | ||
| 370 | } | ||
| 371 | |||
| 372 | intervalRU& operator+=(const intervalRU &x) { | ||
| 373 | // lb += x.lb -> -lbn += -x.lbn -> lbn += x.lbn | ||
| 374 | ln_ += x.ln_; | ||
| 375 | u_ += x.u_; | ||
| 376 | control_add(x); | ||
| 377 | control_check(); | ||
| 378 | return *this; | ||
| 379 | } | ||
| 380 | |||
| 381 | intervalRU& operator-=(const intervalRU &x) { | ||
| 382 | // +=(x.negate()) -> | ||
| 383 | ln_ += x.u_; | ||
| 384 | u_ += x.ln_; | ||
| 385 | control_sub(x); | ||
| 386 | control_check(); | ||
| 387 | return *this; | ||
| 388 | } | ||
| 389 | |||
| 390 | |||
| 391 | intervalRU& operator*=(const intervalRU &b) { | ||
| 392 | // get bounds of both operands | ||
| 393 | double aln = ln_; | ||
| 394 | double au = u_; | ||
| 395 | double bln = b.ln_; | ||
| 396 | double bu = b.u_; | ||
| 397 | |||
| 398 | // negated bounds round to upper | ||
| 399 | // (equivalent to bounds round to lower) | ||
| 400 | double lln = (-aln)*bln; | ||
| 401 | double lun = aln*bu; | ||
| 402 | double uln = au*bln; | ||
| 403 | double uun = (-au)*bu; | ||
| 404 | |||
| 405 | // bounds round to upper | ||
| 406 | double ll = aln*bln; | ||
| 407 | double lu = (-aln)*bu; | ||
| 408 | double ul = au*(-bln); | ||
| 409 | double uu = au*bu; | ||
| 410 | |||
| 411 | ln_ = std::max(std::max(lln,lun),std::max(uln,uun)); | ||
| 412 | u_ = std::max(std::max(ll,lu),std::max(ul,uu)); | ||
| 413 | |||
| 414 | control_mul(b); | ||
| 415 | control_check(); | ||
| 416 | return *this; | ||
| 417 | } | ||
| 418 | |||
| 419 | protected: | ||
| 420 | |||
| 421 | #ifdef INTERVAL_CHECK | ||
| 422 | void control_check() { | ||
| 423 | // expansion_nt used in control_check() operates | ||
| 424 | // in round to nearest mode !! | ||
| 425 | set_FPU_round_to_nearest(); | ||
| 426 | intervalBase::control_check(inf(),sup()); | ||
| 427 | set_FPU_round_to_upper(); | ||
| 428 | } | ||
| 429 | #else | ||
| 430 | void control_check() { | ||
| 431 | } | ||
| 432 | #endif | ||
| 433 | double ln_; | ||
| 434 | double u_; | ||
| 435 | }; | ||
| 436 | |||
| 437 | |||
| 438 | inline intervalRU operator+(const intervalRU& a, const intervalRU& b) { | ||
| 439 | intervalRU result = a; | ||
| 440 | return result += b; | ||
| 441 | } | ||
| 442 | |||
| 443 | inline intervalRU operator-(const intervalRU& a, const intervalRU& b) { | ||
| 444 | intervalRU result = a; | ||
| 445 | return result -= b; | ||
| 446 | } | ||
| 447 | |||
| 448 | inline intervalRU operator*(const intervalRU& a, const intervalRU& b) { | ||
| 449 | intervalRU result = a; | ||
| 450 | return result *= b; | ||
| 451 | } | ||
| 452 | |||
| 453 | /*************************************************************************/ | ||
| 454 | |||
| 455 | /** | ||
| 456 | * \brief Number type for interval arithmetics | ||
| 457 | * \details Interval class in "round to nearest" mode, by Richard Harris: | ||
| 458 | * https://accu.org/journals/overload/19/103/harris_1974/ | ||
| 459 | * Propagates proportional errors at a rate of 1+/-0.5eps | ||
| 460 | * Handles denormals properly (as a special case). | ||
| 461 | */ | ||
| 462 | class intervalRN : public intervalBase { | ||
| 463 | public: | ||
| 464 | |||
| 465 | // operates in default rounding mode | ||
| 466 | // (so Rounding subclass does nothing) | ||
| 467 | struct Rounding { | ||
| 468 | 6081443 | Rounding() { | |
| 469 | 6081443 | } | |
| 470 | 6081443 | ~Rounding() { | |
| 471 | 6081443 | } | |
| 472 | }; | ||
| 473 | |||
| 474 | intervalRN() : | ||
| 475 | intervalBase(), | ||
| 476 | lb_(0.0), | ||
| 477 | ub_(0.0) | ||
| 478 | { | ||
| 479 | control_check(); | ||
| 480 | } | ||
| 481 | |||
| 482 | 7891308 | intervalRN(double x) : | |
| 483 | intervalBase(x), | ||
| 484 | 7891308 | lb_(x), | |
| 485 | 7891308 | ub_(x) | |
| 486 | { | ||
| 487 | 7891308 | control_check(); | |
| 488 | 7891308 | } | |
| 489 | |||
| 490 | intervalRN(double l, double u) : lb_(l), ub_(u) { | ||
| 491 | // note: we cannot control check here. | ||
| 492 | } | ||
| 493 | |||
| 494 | |||
| 495 | intervalRN(const intervalRN& rhs) = default; | ||
| 496 | |||
| 497 | 49409314 | intervalRN(const expansion_nt& rhs) { | |
| 498 | 49409314 | *this = rhs; | |
| 499 | 49409314 | } | |
| 500 | |||
| 501 | intervalRN& operator=(const intervalRN& rhs) = default; | ||
| 502 | |||
| 503 | intervalRN& operator=(double rhs) { | ||
| 504 | lb_ = rhs; | ||
| 505 | ub_ = rhs; | ||
| 506 | control_set(rhs); | ||
| 507 | control_check(); | ||
| 508 | return *this; | ||
| 509 | } | ||
| 510 | |||
| 511 | 49409314 | intervalRN& operator=(const expansion_nt& rhs) { | |
| 512 | |||
| 513 | // Optimized expansion-to-interval conversion: | ||
| 514 | // | ||
| 515 | // Add components starting from the one of largest magnitude | ||
| 516 | // Stop as soon as next component is smaller than ulp (and then | ||
| 517 | // expand interval by ulp). | ||
| 518 | |||
| 519 | 49409314 | index_t l = rhs.length(); | |
| 520 | 49409314 | lb_ = rhs.component(l-1); | |
| 521 | 49409314 | ub_ = rhs.component(l-1); | |
| 522 | |||
| 523 |
2/2✓ Branch 0 taken 36204237 times.
✓ Branch 1 taken 13205077 times.
|
49409314 | for(int comp_idx=int(l)-2; comp_idx>=0; --comp_idx) { |
| 524 | 36204237 | double comp = rhs.component(index_t(comp_idx)); | |
| 525 |
2/2✓ Branch 0 taken 18075053 times.
✓ Branch 1 taken 18129184 times.
|
36204237 | if(comp > 0) { |
| 526 | 18075053 | double nub = ub_ + comp; | |
| 527 |
1/2✓ Branch 0 taken 18075053 times.
✗ Branch 1 not taken.
|
18075053 | if(nub == ub_) { |
| 528 | 18075053 | ub_ = std::nextafter( | |
| 529 | ub_, std::numeric_limits<double>::infinity() | ||
| 530 | ); | ||
| 531 | 18075053 | break; | |
| 532 | } else { | ||
| 533 | ✗ | ub_ = nub; | |
| 534 | ✗ | adjust(); | |
| 535 | } | ||
| 536 | } else { | ||
| 537 | 18129184 | double nlb = lb_ + comp; | |
| 538 |
1/2✓ Branch 0 taken 18129184 times.
✗ Branch 1 not taken.
|
18129184 | if(nlb == lb_) { |
| 539 | 18129184 | lb_ = std::nextafter( | |
| 540 | 18129184 | lb_, -std::numeric_limits<double>::infinity() | |
| 541 | ); | ||
| 542 | 18129184 | break; | |
| 543 | } else { | ||
| 544 | ✗ | lb_ = nlb; | |
| 545 | ✗ | adjust(); | |
| 546 | } | ||
| 547 | } | ||
| 548 | } | ||
| 549 | 49409314 | control_set(rhs); | |
| 550 | 49409314 | control_check(); | |
| 551 | 49409314 | return *this; | |
| 552 | } | ||
| 553 | |||
| 554 | double inf() const { | ||
| 555 | return lb_; | ||
| 556 | } | ||
| 557 | |||
| 558 | double sup() const { | ||
| 559 | return ub_; | ||
| 560 | } | ||
| 561 | |||
| 562 | double estimate() const { | ||
| 563 | return 0.5*(lb_ + ub_); | ||
| 564 | } | ||
| 565 | |||
| 566 | 127489651 | bool is_nan() const { | |
| 567 |
2/4✓ Branch 0 taken 127489651 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 127489651 times.
|
127489651 | return !(lb_==lb_) || !(ub_==ub_); |
| 568 | } | ||
| 569 | |||
| 570 | 9022279 | Sign2 sign() const { | |
| 571 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 9022279 times.
|
9022279 | if(is_nan()) { |
| 572 | ✗ | return SIGN2_ERROR; | |
| 573 | } | ||
| 574 | // Branchless (not sure it is super though...) | ||
| 575 | 9022279 | int lz = int(lb_ == 0); | |
| 576 | 9022279 | int ln = int(lb_ < 0); | |
| 577 | 9022279 | int lp = int(lb_ > 0); | |
| 578 | 9022279 | int uz = int(ub_ == 0); | |
| 579 | 9022279 | int un = int(ub_ < 0); | |
| 580 | 9022279 | int up = int(ub_ > 0); | |
| 581 | 9022279 | Sign2 result = Sign2( | |
| 582 | 9022279 | ln*up*SIGN2_NP+ | |
| 583 | 9022279 | lp*up*SIGN2_PP+ | |
| 584 | 9022279 | lz*up*SIGN2_ZP+ | |
| 585 | 9022279 | ln*un*SIGN2_NN+ | |
| 586 | 9022279 | ln*uz*SIGN2_NZ | |
| 587 | ); | ||
| 588 | 9022279 | result = Sign2( | |
| 589 | int(result) + | ||
| 590 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 9022279 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
9022279 | int(result==SIGN2_ZERO && !(lz&&uz)) * SIGN2_ERROR |
| 591 | ); | ||
| 592 | 9022279 | return result; | |
| 593 | } | ||
| 594 | |||
| 595 | 212819 | intervalRN& negate() { | |
| 596 | 212819 | lb_ = -lb_; | |
| 597 | 212819 | ub_ = -ub_; | |
| 598 | 212819 | std::swap(lb_, ub_); | |
| 599 | 212819 | control_negate(); | |
| 600 | 212819 | control_check(); | |
| 601 | 212819 | return *this; | |
| 602 | } | ||
| 603 | |||
| 604 | 5357103 | intervalRN& operator+=(const intervalRN &x) { | |
| 605 | 5357103 | lb_ += x.lb_; | |
| 606 | 5357103 | ub_ += x.ub_; | |
| 607 | 5357103 | adjust(); | |
| 608 | 5357103 | control_add(x); | |
| 609 | 5357103 | control_check(); | |
| 610 | 5357103 | return *this; | |
| 611 | } | ||
| 612 | |||
| 613 | 31640458 | intervalRN& operator-=(const intervalRN &x) { | |
| 614 | 31640458 | lb_ -= x.ub_; | |
| 615 | 31640458 | ub_ -= x.lb_; | |
| 616 | 31640458 | adjust(); | |
| 617 | 31640458 | control_sub(x); | |
| 618 | 31640458 | control_check(); | |
| 619 | 31640458 | return *this; | |
| 620 | } | ||
| 621 | |||
| 622 | 59233686 | intervalRN& operator*=(const intervalRN &x) { | |
| 623 |
3/6✓ Branch 1 taken 59233686 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 59233686 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 59233686 times.
✗ Branch 7 not taken.
|
59233686 | if(!is_nan() && !x.is_nan()) { |
| 624 | 59233686 | double ll = lb_*x.lb_; | |
| 625 | 59233686 | double lu = lb_*x.ub_; | |
| 626 | 59233686 | double ul = ub_*x.lb_; | |
| 627 | 59233686 | double uu = ub_*x.ub_; | |
| 628 | |||
| 629 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 59233686 times.
|
59233686 | if(!(ll==ll)) ll = 0.0; |
| 630 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 59233686 times.
|
59233686 | if(!(lu==lu)) lu = 0.0; |
| 631 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 59233686 times.
|
59233686 | if(!(ul==ul)) ul = 0.0; |
| 632 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 59233686 times.
|
59233686 | if(!(uu==uu)) uu = 0.0; |
| 633 | |||
| 634 |
2/2✓ Branch 0 taken 21428845 times.
✓ Branch 1 taken 37804841 times.
|
59233686 | if(lu<ll) std::swap(lu, ll); |
| 635 |
2/2✓ Branch 0 taken 11391051 times.
✓ Branch 1 taken 47842635 times.
|
59233686 | if(ul<ll) std::swap(ul, ll); |
| 636 |
2/2✓ Branch 0 taken 8269456 times.
✓ Branch 1 taken 50964230 times.
|
59233686 | if(uu<ll) std::swap(uu, ll); |
| 637 | |||
| 638 |
2/2✓ Branch 0 taken 26435720 times.
✓ Branch 1 taken 32797966 times.
|
59233686 | if(lu>uu) uu = lu; |
| 639 |
2/2✓ Branch 0 taken 8529331 times.
✓ Branch 1 taken 50704355 times.
|
59233686 | if(ul>uu) uu = ul; |
| 640 | |||
| 641 | 59233686 | lb_ = ll; | |
| 642 | 59233686 | ub_ = uu; | |
| 643 | |||
| 644 | 59233686 | adjust(); | |
| 645 | } else { | ||
| 646 | ✗ | lb_ = std::numeric_limits<double>::quiet_NaN(); | |
| 647 | ✗ | ub_ = std::numeric_limits<double>::quiet_NaN(); | |
| 648 | } | ||
| 649 | 59233686 | control_mul(x); | |
| 650 | 59233686 | control_check(); | |
| 651 | 59233686 | return *this; | |
| 652 | } | ||
| 653 | |||
| 654 | protected: | ||
| 655 | |||
| 656 | 96231247 | void adjust() { | |
| 657 | static constexpr double i = std::numeric_limits<double>::infinity(); | ||
| 658 | static constexpr double e = std::numeric_limits<double>::epsilon(); | ||
| 659 | // e = nextafter(1.0) - 1.0 | ||
| 660 | static constexpr double m = std::numeric_limits<double>::min(); | ||
| 661 | // m = smallest normalized | ||
| 662 | static constexpr double l = 1.0-e; | ||
| 663 | static constexpr double u = 1.0+e; | ||
| 664 | static constexpr double em = e*m; | ||
| 665 | |||
| 666 |
6/10✓ Branch 0 taken 96231247 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 96231247 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 14232677 times.
✓ Branch 5 taken 81998570 times.
✓ Branch 6 taken 14232677 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 14232677 times.
✗ Branch 9 not taken.
|
96231247 | if(lb_==lb_ && ub_==ub_ && (lb_!=ub_ || (lb_!=i && lb_!=-i))) { |
| 667 | |||
| 668 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 96231247 times.
|
96231247 | if(lb_>ub_) { |
| 669 | ✗ | std::swap(lb_, ub_); | |
| 670 | } | ||
| 671 | |||
| 672 |
2/2✓ Branch 0 taken 43551013 times.
✓ Branch 1 taken 52680234 times.
|
96231247 | if(lb_>m) { |
| 673 | 43551013 | lb_ *= l; | |
| 674 |
2/2✓ Branch 0 taken 46748176 times.
✓ Branch 1 taken 5932058 times.
|
52680234 | } else if(lb_<-m) { |
| 675 | 46748176 | lb_ *= u; | |
| 676 | } else { | ||
| 677 | 5932058 | lb_ -= em; | |
| 678 | } | ||
| 679 | |||
| 680 |
2/2✓ Branch 0 taken 47242974 times.
✓ Branch 1 taken 48988273 times.
|
96231247 | if(ub_>m) { |
| 681 | 47242974 | ub_ *= u; | |
| 682 |
2/2✓ Branch 0 taken 43056262 times.
✓ Branch 1 taken 5932011 times.
|
48988273 | } else if(ub_<-m) { |
| 683 | 43056262 | ub_ *= l; | |
| 684 | } else { | ||
| 685 | 5932011 | ub_ += em; | |
| 686 | } | ||
| 687 | } else { | ||
| 688 | ✗ | lb_ = std::numeric_limits<double>::quiet_NaN(); | |
| 689 | ✗ | ub_ = std::numeric_limits<double>::quiet_NaN(); | |
| 690 | } | ||
| 691 | 96231247 | } | |
| 692 | |||
| 693 | #ifdef INTERVAL_CHECK | ||
| 694 | void control_check() { | ||
| 695 | intervalBase::control_check(inf(),sup()); | ||
| 696 | } | ||
| 697 | #else | ||
| 698 | 153744688 | void control_check() { | |
| 699 | 153744688 | } | |
| 700 | #endif | ||
| 701 | |||
| 702 | private: | ||
| 703 | double lb_; /**< lower bound */ | ||
| 704 | double ub_; /**< upper bound */ | ||
| 705 | }; | ||
| 706 | |||
| 707 | 5357103 | inline intervalRN operator+(const intervalRN& a, const intervalRN& b) { | |
| 708 | 5357103 | intervalRN result = a; | |
| 709 |
1/2✓ Branch 1 taken 5357103 times.
✗ Branch 2 not taken.
|
10714206 | return result += b; |
| 710 | } | ||
| 711 | |||
| 712 | 31640458 | inline intervalRN operator-(const intervalRN& a, const intervalRN& b) { | |
| 713 | 31640458 | intervalRN result = a; | |
| 714 |
1/2✓ Branch 1 taken 31640458 times.
✗ Branch 2 not taken.
|
63280916 | return result -= b; |
| 715 | } | ||
| 716 | |||
| 717 | 59233686 | inline intervalRN operator*(const intervalRN& a, const intervalRN& b) { | |
| 718 | 59233686 | intervalRN result = a; | |
| 719 |
1/2✓ Branch 1 taken 59233686 times.
✗ Branch 2 not taken.
|
118467372 | return result *= b; |
| 720 | } | ||
| 721 | |||
| 722 | /**********************************************************************/ | ||
| 723 | |||
| 724 | /** | ||
| 725 | * \brief Dummy interval class, that does not compute anything and that | ||
| 726 | * always says that sign is undetermined. | ||
| 727 | * \details Used for performance testing and for debugging. | ||
| 728 | */ | ||
| 729 | class intervalDummy : public intervalBase { | ||
| 730 | public: | ||
| 731 | struct Rounding { | ||
| 732 | Rounding() { | ||
| 733 | geo_argused(this); // silence a warning in predicate filters. | ||
| 734 | } | ||
| 735 | }; | ||
| 736 | |||
| 737 | intervalDummy() { | ||
| 738 | } | ||
| 739 | |||
| 740 | intervalDummy(double x) { | ||
| 741 | geo_argused(x); | ||
| 742 | } | ||
| 743 | |||
| 744 | intervalDummy(double l, double u) { | ||
| 745 | geo_argused(l); | ||
| 746 | geo_argused(u); | ||
| 747 | } | ||
| 748 | |||
| 749 | intervalDummy(const intervalDummy& rhs) = default; | ||
| 750 | |||
| 751 | intervalDummy(const expansion_nt& rhs) { | ||
| 752 | geo_argused(rhs); | ||
| 753 | } | ||
| 754 | |||
| 755 | intervalDummy& operator=(const intervalDummy& rhs) = default; | ||
| 756 | |||
| 757 | intervalDummy& operator=(double rhs) { | ||
| 758 | geo_argused(rhs); | ||
| 759 | return *this; | ||
| 760 | } | ||
| 761 | |||
| 762 | intervalDummy& operator=(const expansion_nt& rhs) { | ||
| 763 | geo_argused(rhs); | ||
| 764 | return *this; | ||
| 765 | } | ||
| 766 | |||
| 767 | double inf() const { | ||
| 768 | return -std::numeric_limits<double>::infinity(); | ||
| 769 | } | ||
| 770 | |||
| 771 | double sup() const { | ||
| 772 | return std::numeric_limits<double>::infinity(); | ||
| 773 | } | ||
| 774 | |||
| 775 | double estimate() const { | ||
| 776 | return std::numeric_limits<double>::quiet_NaN(); | ||
| 777 | } | ||
| 778 | |||
| 779 | bool is_nan() const { | ||
| 780 | return true; | ||
| 781 | } | ||
| 782 | |||
| 783 | Sign2 sign() const { | ||
| 784 | return SIGN2_NP; | ||
| 785 | } | ||
| 786 | |||
| 787 | intervalDummy& negate() { | ||
| 788 | return *this; | ||
| 789 | } | ||
| 790 | |||
| 791 | intervalDummy& operator+=(const intervalDummy &x) { | ||
| 792 | geo_argused(x); | ||
| 793 | return *this; | ||
| 794 | } | ||
| 795 | |||
| 796 | intervalDummy& operator-=(const intervalDummy &x) { | ||
| 797 | geo_argused(x); | ||
| 798 | return *this; | ||
| 799 | } | ||
| 800 | |||
| 801 | intervalDummy& operator*=(const intervalDummy &x) { | ||
| 802 | geo_argused(x); | ||
| 803 | return *this; | ||
| 804 | } | ||
| 805 | }; | ||
| 806 | |||
| 807 | inline intervalDummy operator+( | ||
| 808 | const intervalDummy& a, const intervalDummy& b | ||
| 809 | ) { | ||
| 810 | intervalDummy result = a; | ||
| 811 | return result += b; | ||
| 812 | } | ||
| 813 | |||
| 814 | inline intervalDummy operator-( | ||
| 815 | const intervalDummy& a, const intervalDummy& b | ||
| 816 | ) { | ||
| 817 | intervalDummy result = a; | ||
| 818 | return result -= b; | ||
| 819 | } | ||
| 820 | |||
| 821 | inline intervalDummy operator*( | ||
| 822 | const intervalDummy& a, const intervalDummy& b | ||
| 823 | ) { | ||
| 824 | intervalDummy result = a; | ||
| 825 | return result *= b; | ||
| 826 | } | ||
| 827 | |||
| 828 | /**********************************************************************/ | ||
| 829 | |||
| 830 | #ifdef GEO_NO_INTERVALS | ||
| 831 | typedef intervalDummy interval_nt; | ||
| 832 | #else | ||
| 833 | // interval_nt uses round to nearet (RN) | ||
| 834 | typedef intervalRN interval_nt; | ||
| 835 | |||
| 836 | // one could use round to upper (RU) instead, supposed to be slightly faster | ||
| 837 | // but we chose not to use it because: | ||
| 838 | // - does not seem to be significantly faster | ||
| 839 | // - makes the code more complicated, needs to change rounding mode) | ||
| 840 | // (but it is handled properly with the Rounding class provided that | ||
| 841 | // client code does not forget to declare one) | ||
| 842 | // - needs the compilation flag -frounding-math (disable compiler | ||
| 843 | // optimizations that suppose "round to nearest" mode) | ||
| 844 | // typedef intervalRU interval_nt; | ||
| 845 | #endif | ||
| 846 | |||
| 847 | |||
| 848 | /** \brief Specialization of GEO::is_scalar */ | ||
| 849 | template <> struct is_scalar<interval_nt> { | ||
| 850 | typedef interval_nt type; | ||
| 851 | static constexpr bool value = true; | ||
| 852 | }; | ||
| 853 | |||
| 854 | /** | ||
| 855 | * \brief Specialization of 2d homogeneous-coord vector difference | ||
| 856 | * for interval_nt coordinates | ||
| 857 | * \details We cannot check for p1.w and p2.w equality (since they | ||
| 858 | * might be in the same interval without being equal), so we always | ||
| 859 | * reduce to the same denominator. | ||
| 860 | */ | ||
| 861 | 1029186 | template <> inline vec2Hg<interval_nt> operator-( | |
| 862 | const vec2Hg<interval_nt>& p1, const vec2Hg<interval_nt>& p2 | ||
| 863 | ) { | ||
| 864 | return vec2Hg<interval_nt>( | ||
| 865 |
1/2✓ Branch 1 taken 1029186 times.
✗ Branch 2 not taken.
|
1029186 | det2x2(p1.x,p1.w,p2.x,p2.w), |
| 866 |
1/2✓ Branch 1 taken 1029186 times.
✗ Branch 2 not taken.
|
1029186 | det2x2(p1.y,p1.w,p2.y,p2.w), |
| 867 |
1/2✓ Branch 1 taken 1029186 times.
✗ Branch 2 not taken.
|
2058372 | p1.w*p2.w |
| 868 | 2058372 | ); | |
| 869 | } | ||
| 870 | |||
| 871 | |||
| 872 | /** | ||
| 873 | * \brief Specialization of 3d homogeneous-coord vector difference | ||
| 874 | * for interval_nt coordinates | ||
| 875 | * \details We cannot check for p1.w and p2.w equality (since they | ||
| 876 | * might be in the same interval without being equal), so we always | ||
| 877 | * reduce to the same denominator. | ||
| 878 | */ | ||
| 879 | 71229 | template <> inline vec3Hg<interval_nt> operator-( | |
| 880 | const vec3Hg<interval_nt>& p1, const vec3Hg<interval_nt>& p2 | ||
| 881 | ) { | ||
| 882 | return vec3Hg<interval_nt>( | ||
| 883 |
1/2✓ Branch 1 taken 71229 times.
✗ Branch 2 not taken.
|
71229 | det2x2(p1.x,p1.w,p2.x,p2.w), |
| 884 |
1/2✓ Branch 1 taken 71229 times.
✗ Branch 2 not taken.
|
71229 | det2x2(p1.y,p1.w,p2.y,p2.w), |
| 885 |
1/2✓ Branch 1 taken 71229 times.
✗ Branch 2 not taken.
|
71229 | det2x2(p1.z,p1.w,p2.z,p2.w), |
| 886 |
1/2✓ Branch 1 taken 71229 times.
✗ Branch 2 not taken.
|
142458 | p1.w * p2.w |
| 887 | 142458 | ); | |
| 888 | } | ||
| 889 | } | ||
| 890 | |||
| 891 | #endif | ||
| 892 |