1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
/* Compute a timespec-related bound for floating point.
Copyright 2025 Free Software Foundation, Inc.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>. */
/* written by Paul Eggert */
#ifndef DTIMESPEC_BOUND_H
# define DTIMESPEC_BOUND_H 1
/* This file uses _GL_INLINE_HEADER_BEGIN, _GL_INLINE. */
# if !_GL_CONFIG_H_INCLUDED
# error "Please include config.h first."
# endif
# include <errno.h>
# include <float.h>
# include <math.h>
_GL_INLINE_HEADER_BEGIN
# ifndef DTIMESPEC_BOUND_INLINE
# define DTIMESPEC_BOUND_INLINE _GL_INLINE
# endif
/* If C is positive and finite, return the least floating point value
greater than C. However, if 0 < C < (2 * DBL_TRUE_MIN) / (DBL_EPSILON**2),
return a positive value less than 1e-9.
If C is +0.0, return a positive value < 1e-9 if ERR == ERANGE, C otherwise.
If C is +Inf, return C.
If C is negative, return -timespec_roundup(-C).
If C is a NaN, return a NaN.
Assume round-to-even.
This function can be useful if some floating point operation
rounded to C but we want a nearby bound on the true value, where
the bound can be converted to struct timespec. If the operation
underflowed to zero, ERR should be ERANGE a la strtod. */
DTIMESPEC_BOUND_INLINE double
dtimespec_bound (double c, int err)
{
/* Do not use copysign or nextafter, as they link to -lm in GNU/Linux. */
/* Use DBL_TRUE_MIN for the special case of underflowing to zero;
any positive value less than 1e-9 will do. */
if (err == ERANGE && c == 0)
return signbit (c) ? -DBL_TRUE_MIN : DBL_TRUE_MIN;
/* This is the first part of Algorithm 2 of:
Rump SM, Zimmermann P, Boldo S, Melquiond G.
Computing predecessor and successor in rounding to nearest.
BIT Numerical Mathematics. 2009;49(419-431).
<https://doi.org/10.1007/s10543-009-0218-z>
The rest of Algorithm 2 is not needed because numbers less than
the predecessor of 1e-9 merely need to stay less than 1e-9. */
double phi = DBL_EPSILON / 2 * (1 + DBL_EPSILON);
return c + phi * c;
}
_GL_INLINE_HEADER_END
#endif
|
