76 lines
2.6 KiB
C
76 lines
2.6 KiB
C
/* 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
|