| 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_BASIC_NUMERIC | ||
| 41 | #define GEOGRAM_BASIC_NUMERIC | ||
| 42 | |||
| 43 | #include <geogram/basic/common.h> | ||
| 44 | #include <geogram/basic/assert.h> | ||
| 45 | #include <cmath> | ||
| 46 | #include <float.h> | ||
| 47 | #include <limits.h> | ||
| 48 | #include <algorithm> // for std::min / std::max | ||
| 49 | #include <stdint.h> | ||
| 50 | #include <limits> | ||
| 51 | #include <type_traits> | ||
| 52 | #include <iostream> | ||
| 53 | #include <cstdlib> | ||
| 54 | |||
| 55 | #ifndef M_PI | ||
| 56 | /** | ||
| 57 | * \brief Value of the constant PI if not defined by the system | ||
| 58 | */ | ||
| 59 | #define M_PI 3.14159265358979323846 | ||
| 60 | #endif | ||
| 61 | |||
| 62 | /** | ||
| 63 | * \file geogram/basic/numeric.h | ||
| 64 | * \brief Types and functions for numbers manipulation | ||
| 65 | */ | ||
| 66 | |||
| 67 | namespace GEO { | ||
| 68 | |||
| 69 | /** | ||
| 70 | * \brief Integer constants that represent the sign of a value | ||
| 71 | */ | ||
| 72 | enum Sign { | ||
| 73 | /** Value is negative */ | ||
| 74 | NEGATIVE = -1, | ||
| 75 | /** Value is zero */ | ||
| 76 | ZERO = 0, | ||
| 77 | /** Value is positive */ | ||
| 78 | POSITIVE = 1 | ||
| 79 | }; | ||
| 80 | |||
| 81 | |||
| 82 | /** | ||
| 83 | * \brief Compares two values | ||
| 84 | * \param[in] a , b the two values to compare | ||
| 85 | * \tparam T the type of the value | ||
| 86 | * \retval POSITIVE if \p a is greater than \p b | ||
| 87 | * \retval ZERO if \p a is equal to \p b | ||
| 88 | * \retval NEGATIVE if \p a is smaller than \p b | ||
| 89 | * \see Sign | ||
| 90 | */ | ||
| 91 | template <class T> | ||
| 92 | 172496177 | inline Sign geo_cmp(const T& a, const T& b) { | |
| 93 | 172496177 | return Sign((a > b) - (a < b)); | |
| 94 | } | ||
| 95 | |||
| 96 | |||
| 97 | /** | ||
| 98 | * \brief Gets the sign of a value | ||
| 99 | * \details Returns -1, 0, or 1 whether value \p x is resp. negative, zero | ||
| 100 | * or positive. The function uses operator<() and operator>() to compare | ||
| 101 | * the value to 0 (zero). The integer constant zero must make | ||
| 102 | * senses for the type of the value, or T must be constructible from | ||
| 103 | * integer constant zero. | ||
| 104 | * \param[in] x the value to test | ||
| 105 | * \tparam T the type of the value | ||
| 106 | * \return the sign of the value | ||
| 107 | * \see Sign | ||
| 108 | */ | ||
| 109 | template <class T> | ||
| 110 | 172418651 | inline Sign geo_sgn(const T& x) { | |
| 111 | 172418651 | return geo_cmp(x, T(0)); | |
| 112 | } | ||
| 113 | |||
| 114 | /** | ||
| 115 | * \brief Defines numeric types used in Vorpaline. | ||
| 116 | * \details | ||
| 117 | * These types names have the form (u)int<size> or float<size>, | ||
| 118 | * where the (optional) u denotes an unsigned type, | ||
| 119 | * and the size is in bits. | ||
| 120 | */ | ||
| 121 | namespace Numeric { | ||
| 122 | |||
| 123 | /** Generic pointer type */ | ||
| 124 | typedef void* pointer; | ||
| 125 | |||
| 126 | /** Integer type with a width of 8 bits */ | ||
| 127 | typedef int8_t int8; | ||
| 128 | |||
| 129 | /** Integer type with a width of 16 bits */ | ||
| 130 | typedef int16_t int16; | ||
| 131 | |||
| 132 | /** Integer type with a width of 32 bits */ | ||
| 133 | typedef int32_t int32; | ||
| 134 | |||
| 135 | /** Integer type with a width of 64 bits */ | ||
| 136 | typedef int64_t int64; | ||
| 137 | |||
| 138 | /** Unsigned integer type with a width of 8 bits */ | ||
| 139 | typedef uint8_t uint8; | ||
| 140 | |||
| 141 | /** Unsigned integer type with a width of 16 bits */ | ||
| 142 | typedef uint16_t uint16; | ||
| 143 | |||
| 144 | /** Unsigned integer type with a width of 32 bits */ | ||
| 145 | typedef uint32_t uint32; | ||
| 146 | |||
| 147 | /** Unsigned integer type with a width of 64 bits */ | ||
| 148 | typedef uint64_t uint64; | ||
| 149 | |||
| 150 | /** Floating point type with a width of 32 bits */ | ||
| 151 | typedef float float32; | ||
| 152 | |||
| 153 | /** Floating point type with a width of 64 bits */ | ||
| 154 | typedef double float64; | ||
| 155 | |||
| 156 | /** | ||
| 157 | * \brief Gets 32 bits float maximum positive value | ||
| 158 | */ | ||
| 159 | ✗ | inline constexpr float32 max_float32() { | |
| 160 | ✗ | return std::numeric_limits<float32>::max(); | |
| 161 | } | ||
| 162 | |||
| 163 | /** | ||
| 164 | * \brief Gets 32 bits float minimum negative value | ||
| 165 | */ | ||
| 166 | ✗ | inline constexpr float32 min_float32() { | |
| 167 | // Note: numeric_limits<>::min() is not | ||
| 168 | // what we want (it returns the smallest | ||
| 169 | // positive non-denormal). | ||
| 170 | ✗ | return -max_float32(); | |
| 171 | } | ||
| 172 | |||
| 173 | /** | ||
| 174 | * \brief Gets 64 bits float maximum positive value | ||
| 175 | */ | ||
| 176 | 93933800 | inline constexpr float64 max_float64() { | |
| 177 | 93933800 | return std::numeric_limits<float64>::max(); | |
| 178 | } | ||
| 179 | |||
| 180 | /** | ||
| 181 | * \brief Gets 64 bits float minimum negative value | ||
| 182 | */ | ||
| 183 | 679495 | inline constexpr float64 min_float64() { | |
| 184 | // Note: numeric_limits<>::min() is not | ||
| 185 | // what we want (it returns the smallest | ||
| 186 | // positive non-denormal). | ||
| 187 | 679495 | return -max_float64(); | |
| 188 | } | ||
| 189 | |||
| 190 | /** | ||
| 191 | * \brief Checks whether a 32 bits float is "not a number" | ||
| 192 | */ | ||
| 193 | bool GEOGRAM_API is_nan(float32 x); | ||
| 194 | |||
| 195 | /** | ||
| 196 | * \brief Checks whether a 64 bits float is "not a number" | ||
| 197 | */ | ||
| 198 | bool GEOGRAM_API is_nan(float64 x); | ||
| 199 | |||
| 200 | /** | ||
| 201 | * \brief Resets the random number generator. | ||
| 202 | * \details Uses "algo:random_seed" | ||
| 203 | */ | ||
| 204 | void GEOGRAM_API random_reset(); | ||
| 205 | |||
| 206 | /** | ||
| 207 | * \brief Resets the random number generator. | ||
| 208 | * \param[in] seed the random seed or -1 to use default | ||
| 209 | */ | ||
| 210 | void GEOGRAM_API random_reset(int seed); | ||
| 211 | |||
| 212 | /** | ||
| 213 | * \brief Returns a 32 bits integer between 0 and RAND_MAX | ||
| 214 | */ | ||
| 215 | int32 GEOGRAM_API random_int32(); | ||
| 216 | |||
| 217 | /** | ||
| 218 | * \brief Returns a 32 bits float between 0 and 1 | ||
| 219 | */ | ||
| 220 | float32 GEOGRAM_API random_float32(); | ||
| 221 | |||
| 222 | /** | ||
| 223 | * \brief Returns a 64 bits float between 0 and 1 | ||
| 224 | */ | ||
| 225 | float64 GEOGRAM_API random_float64(); | ||
| 226 | |||
| 227 | /** | ||
| 228 | * \brief Limits helper class that extends std::numeric_limits | ||
| 229 | * \details LimitsHelper extends std::numeric_limits to provide | ||
| 230 | * additional information about numeric types \p T. | ||
| 231 | * Template parameter \p is_numeric receives the value \c | ||
| 232 | * std::numeric_limits<T>::is_specialized which is \c true for all | ||
| 233 | * numeric types and \c false for the other types. The template is | ||
| 234 | * specialized for \p is_numeric == \c true to define additional | ||
| 235 | * information. For non-numeric types, the default template does not | ||
| 236 | * define anything. | ||
| 237 | * \tparam T an object type | ||
| 238 | * \tparam is_numeric is true if type \p T is a numeric type, false | ||
| 239 | * otherwise. | ||
| 240 | */ | ||
| 241 | template <class T, bool is_numeric> | ||
| 242 | struct LimitsHelper : std::numeric_limits<T> { | ||
| 243 | }; | ||
| 244 | |||
| 245 | /** | ||
| 246 | * \brief Specialization of LimitsHelper for numeric types | ||
| 247 | * \details This specialization defines the following values: | ||
| 248 | * - size - the size of the numeric type in bytes | ||
| 249 | * - numbits - the size of the numeric type in bits | ||
| 250 | * \tparam T a numeric type | ||
| 251 | */ | ||
| 252 | template <class T> | ||
| 253 | struct LimitsHelper<T, true> : std::numeric_limits<T> { | ||
| 254 | /** The size of the numeric type in bytes */ | ||
| 255 | static const size_t size = sizeof(T); | ||
| 256 | /** The size of the numeric type in bits */ | ||
| 257 | static const size_t numbits = 8 * sizeof(T); | ||
| 258 | }; | ||
| 259 | |||
| 260 | /** | ||
| 261 | * \brief Extends std::numeric_limits with additional information | ||
| 262 | * \details Limits provides additional information about numeric types | ||
| 263 | * that are not available in std::numeric_limits: | ||
| 264 | * - size: the size of the numeric type in bytes | ||
| 265 | * - numbits: the size of the numeric type in bits | ||
| 266 | * These types are defined in the helper class LimitsHelper for | ||
| 267 | * numeric types only. They are not defined for non-numeric types. | ||
| 268 | */ | ||
| 269 | template <class T> | ||
| 270 | struct Limits : | ||
| 271 | LimitsHelper<T, std::numeric_limits<T>::is_specialized> { | ||
| 272 | }; | ||
| 273 | |||
| 274 | /** | ||
| 275 | * \brief place holder for optimizing internal number representation | ||
| 276 | * \details there are specializations for expansion_nt, rational_nt | ||
| 277 | */ | ||
| 278 | template <class T> inline void optimize_number_representation(T& x) { | ||
| 279 | geo_argused(x); | ||
| 280 | } | ||
| 281 | |||
| 282 | /** | ||
| 283 | * \brief Compares two rational numbers given as separate | ||
| 284 | * numerators and denominators. | ||
| 285 | * \param[in] a_num , a_denom defines a = \p a_num / \p a_denom | ||
| 286 | * \param[in] b_num , b_denom defines b = \p b_num / \p b_denom | ||
| 287 | * \return the sign of a - b | ||
| 288 | */ | ||
| 289 | template <class T> inline Sign ratio_compare( | ||
| 290 | const T& a_num, const T& a_denom, const T& b_num, const T& b_denom | ||
| 291 | ) { | ||
| 292 | if(a_denom == b_denom) { | ||
| 293 | return Sign(geo_cmp(a_num,b_num)*geo_sgn(a_denom)); | ||
| 294 | } | ||
| 295 | return Sign( | ||
| 296 | geo_cmp(a_num*b_denom, b_num*a_denom) * | ||
| 297 | geo_sgn(a_denom) * geo_sgn(b_denom) | ||
| 298 | ); | ||
| 299 | } | ||
| 300 | } | ||
| 301 | |||
| 302 | /************************************************************************/ | ||
| 303 | |||
| 304 | |||
| 305 | /** | ||
| 306 | * \brief Gets the square value of a value | ||
| 307 | * \param[in] x a value of type \p T | ||
| 308 | * \tparam T the type of the value | ||
| 309 | * \return the square value of \p x | ||
| 310 | */ | ||
| 311 | template <class T> | ||
| 312 | 3483165777 | inline T geo_sqr(T x) { | |
| 313 | 3483165777 | return x * x; | |
| 314 | } | ||
| 315 | |||
| 316 | /** | ||
| 317 | * \brief Clamps a value to a range | ||
| 318 | * \details Clamps the value \p x to a range defined by \p min and \p max. | ||
| 319 | * This modifies the value of \p x directly. | ||
| 320 | * \param[in,out] x a value of type \p T | ||
| 321 | * \param[in] min the lower bound of the clamping range | ||
| 322 | * \param[in] max the upper bound of the clamping range | ||
| 323 | */ | ||
| 324 | template <class T> | ||
| 325 | ✗ | inline void geo_clamp(T& x, T min, T max) { | |
| 326 | ✗ | if(x < min) { | |
| 327 | ✗ | x = min; | |
| 328 | ✗ | } else if(x > max) { | |
| 329 | ✗ | x = max; | |
| 330 | } | ||
| 331 | ✗ | } | |
| 332 | |||
| 333 | /** | ||
| 334 | * \brief The type for storing and manipulating indices. | ||
| 335 | * \internal | ||
| 336 | * Vorpaline uses 32 bit indices (can be changed to 64 bits | ||
| 337 | * if need be, but this will double memory consumption of | ||
| 338 | * all combinatorial data structures). | ||
| 339 | */ | ||
| 340 | typedef geo_index_t index_t; | ||
| 341 | |||
| 342 | /** | ||
| 343 | * \brief Gets the maximum positive value of type index_t. | ||
| 344 | */ | ||
| 345 | 1437 | inline constexpr index_t max_index_t() { | |
| 346 | 1437 | return std::numeric_limits<index_t>::max(); | |
| 347 | } | ||
| 348 | |||
| 349 | /** | ||
| 350 | * \brief The type for storing and manipulating indices differences. | ||
| 351 | * \details Can be negative (for instance to indicate special values like | ||
| 352 | * borders). | ||
| 353 | */ | ||
| 354 | typedef geo_signed_index_t signed_index_t; | ||
| 355 | |||
| 356 | /** | ||
| 357 | * \brief Gets the maximum positive value of type signed_index_t. | ||
| 358 | */ | ||
| 359 | inline constexpr signed_index_t max_signed_index_t() { | ||
| 360 | return std::numeric_limits<signed_index_t>::max(); | ||
| 361 | } | ||
| 362 | |||
| 363 | /** | ||
| 364 | * \brief Gets the minimum negative value of type signed_index_t. | ||
| 365 | */ | ||
| 366 | inline signed_index_t min_signed_index_t() { | ||
| 367 | return std::numeric_limits<signed_index_t>::min(); | ||
| 368 | } | ||
| 369 | |||
| 370 | /** | ||
| 371 | * \brief The type for storing coordinate indices, and iterating on | ||
| 372 | * the coordinates of a point. | ||
| 373 | */ | ||
| 374 | typedef geo_coord_index_t coord_index_t; | ||
| 375 | |||
| 376 | /** | ||
| 377 | * \TODOC | ||
| 378 | */ | ||
| 379 | 24 | inline double round(double x) { | |
| 380 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
|
24 | return ((x - floor(x)) > 0.5 ? ceil(x) : floor(x)); |
| 381 | } | ||
| 382 | |||
| 383 | /************************************************************************/ | ||
| 384 | |||
| 385 | /** | ||
| 386 | * \brief The dummy index value. | ||
| 387 | * \details Used for instance on the border of a surface, where adjacent | ||
| 388 | * facets are set to NO_INDEX. | ||
| 389 | */ | ||
| 390 | static constexpr index_t NO_INDEX = index_t(-1); | ||
| 391 | |||
| 392 | /************************************************************************/ | ||
| 393 | |||
| 394 | /** | ||
| 395 | * \brief type traits for scalars | ||
| 396 | * \details all basic C++ arithmetic types plus expansion_nt, interval_nt, | ||
| 397 | * rational_nt and geogramplus exact_nt. | ||
| 398 | * Used to avoid ambiguous declarations in scalar * vector products. | ||
| 399 | */ | ||
| 400 | template <class T> struct is_scalar { | ||
| 401 | typedef typename std::is_arithmetic<T>::type type; | ||
| 402 | static constexpr bool value = std::is_arithmetic<T>::value; | ||
| 403 | }; | ||
| 404 | |||
| 405 | /************************************************************************/ | ||
| 406 | } | ||
| 407 | |||
| 408 | |||
| 409 | /************* Disable floating point contraction **************************/ | ||
| 410 | |||
| 411 | /** | ||
| 412 | * \brief add GEO_FP_CONTRACT_OFF at the beginning of files or functions | ||
| 413 | * where FMA operations should not be generated in place of a*b+c | ||
| 414 | * \details it is important for instance for geometric predicates, that | ||
| 415 | * rely on strict IEEE754 implementation of product and addition. | ||
| 416 | */ | ||
| 417 | |||
| 418 | #if defined(GOMGEN) | ||
| 419 | # define GEO_FP_CONTRACT_OFF(x) | ||
| 420 | #elif defined(__clang__) | ||
| 421 | # define GEO_FP_CONTRACT_OFF _Pragma("clang fp contract(off)") | ||
| 422 | #elif defined(_MSC_VER) | ||
| 423 | # define GEO_FP_CONTRACT_OFF _Pragma("fp_contract(off)") | ||
| 424 | #elif defined(__GNUC__) | ||
| 425 | |||
| 426 | // GCC does not have any pragma to deactivate FMA generation, | ||
| 427 | // so instead we check that they are deactivated (by the command-line | ||
| 428 | // option -ffp-contract=off) and fire an assertion fail if it was not | ||
| 429 | // the case. | ||
| 430 | struct GeoAssertNoFpContract { | ||
| 431 | 504 | GeoAssertNoFpContract() { | |
| 432 | #ifdef GEOGRAM_PSM | ||
| 433 | if(fp_contraction_enabled()) { | ||
| 434 | std::cerr << "Needs to be compiled with -ffp-contract-off" | ||
| 435 | << std::endl; | ||
| 436 | abort(); | ||
| 437 | } | ||
| 438 | #else | ||
| 439 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 504 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
504 | geo_assert(!fp_contraction_enabled()); |
| 440 | #endif | ||
| 441 | 504 | } | |
| 442 | 504 | static bool fp_contraction_enabled() { | |
| 443 | 504 | return (a2plusb(0x1.0000002p0, -0x1.0000004p0) != 0.0); | |
| 444 | } | ||
| 445 | 504 | __attribute__((noipa)) static double a2plusb(double a, double b) { | |
| 446 | 504 | return a * a + b; | |
| 447 | } | ||
| 448 | }; | ||
| 449 | # define GEO_FP_CONTRACT_OFF \ | ||
| 450 | static GeoAssertNoFpContract CPP_CONCAT(assert_no_fp_contract_,__LINE__); | ||
| 451 | #else | ||
| 452 | # define GEO_FP_CONTRACT_OFF _Pragma("STDC FP_CONTRACT OFF") | ||
| 453 | #endif | ||
| 454 | |||
| 455 | #endif | ||
| 456 |