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