GCC Code Coverage Report


Directory: ./
File: lib/geogram/basic/numeric.h
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 10 14 71.4%
Functions: 3 3 100.0%
Branches: 62 134 46.3%

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 inline Sign geo_cmp(const T& a, const T& b) {
93
51/78
✓ Branch 0 taken 821960 times.
✓ Branch 1 taken 9970737 times.
✓ Branch 2 taken 7157014 times.
✓ Branch 3 taken 12171350 times.
✓ Branch 4 taken 322626 times.
✓ Branch 5 taken 70443 times.
✓ Branch 6 taken 16002121 times.
✓ Branch 7 taken 33203951 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 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 180822 times.
✓ Branch 40 taken 926 times.
✓ Branch 41 taken 89340 times.
✓ Branch 42 taken 30306 times.
✓ Branch 43 taken 303 times.
✓ Branch 44 taken 30345 times.
✓ Branch 45 taken 29984 times.
✓ Branch 46 taken 30367 times.
✓ Branch 47 taken 59637 times.
✗ 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 not taken.
✗ Branch 58 not taken.
✓ Branch 59 taken 80655 times.
✓ Branch 60 taken 72918 times.
✓ Branch 61 taken 7737 times.
✓ Branch 62 taken 15601 times.
✓ Branch 63 taken 20314 times.
✓ Branch 64 taken 10393 times.
✓ Branch 65 taken 141 times.
✓ Branch 66 taken 35155 times.
✓ Branch 67 taken 10677 times.
✗ Branch 68 not taken.
✓ Branch 69 taken 279520 times.
✓ Branch 70 taken 271431 times.
✓ Branch 71 taken 8089 times.
✓ Branch 72 taken 96912 times.
✓ Branch 73 taken 32953 times.
✓ Branch 74 taken 114073 times.
✓ Branch 75 taken 20264 times.
✓ Branch 76 taken 204580 times.
✓ Branch 77 taken 2450 times.
148026577 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 inline Sign geo_sgn(const T& x) {
111 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 inline constexpr float64 max_float64() {
177 return std::numeric_limits<float64>::max();
178 }
179
180 /**
181 * \brief Gets 64 bits float minimum negative value
182 */
183 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 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 inline T geo_sqr(T x) {
313
9/30
✓ Branch 0 taken 15822252 times.
✓ Branch 1 taken 1665852 times.
✓ Branch 2 taken 15250638 times.
✓ Branch 3 taken 1565330 times.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 292417 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 371631 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 371631 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 79295 times.
✗ 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.
2827940641 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 inline constexpr index_t max_index_t() {
346 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 0 not taken.
✓ Branch 1 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 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