diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/s_nextupl.c')
-rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_nextupl.c | 78 |
1 files changed, 78 insertions, 0 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_nextupl.c b/sysdeps/ieee754/ldbl-128ibm/s_nextupl.c new file mode 100644 index 0000000..2a4fae7 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/s_nextupl.c @@ -0,0 +1,78 @@ +/* Return the least floating-point number greater than X. + Copyright (C) 2016 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library 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 + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <math_ldbl_opt.h> + +/* Return the least floating-point number greater than X. */ +long double +__nextupl (long double x) +{ + int64_t hx, ihx, lx; + double xhi, xlo, yhi; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + ihx = hx & 0x7fffffffffffffffLL; + + if (ihx > 0x7ff0000000000000LL) /* x is nan. */ + return x + x; /* Signal the nan. */ + if (ihx == 0) + return LDBL_TRUE_MIN; + + long double u; + if ((hx == 0x7fefffffffffffffLL) && (lx == 0x7c8ffffffffffffeLL)) + return INFINITY; + if ((uint64_t) hx >= 0xfff0000000000000ULL) + { + u = -0x1.fffffffffffff7ffffffffffff8p+1023L; + return u; + } + if (ihx <= 0x0360000000000000LL) + { /* x <= LDBL_MIN. */ + x += LDBL_TRUE_MIN; + if (x == 0.0L) /* Handle negative LDBL_TRUE_MIN case. */ + x = -0.0L; + return x; + } + /* If the high double is an exact power of two and the low + double is the opposite sign, then 1ulp is one less than + what we might determine from the high double. Similarly + if X is an exact power of two, and negative, because + making it a little larger will result in the exponent + decreasing by one and normalisation of the mantissa. */ + if ((hx & 0x000fffffffffffffLL) == 0 + && ((lx != 0 && lx != 0x8000000000000000LL && (hx ^ lx) < 0) + || ((lx == 0 || lx == 0x8000000000000000LL) && hx < 0))) + ihx -= 1LL << 52; + if (ihx < (106LL << 52)) + { /* ulp will denormal. */ + INSERT_WORDS64 (yhi, ihx & (0x7ffLL << 52)); + u = yhi * 0x1p-105; + } + else + { + INSERT_WORDS64 (yhi, (ihx & (0x7ffLL << 52)) - (105LL << 52)); + u = yhi; + } + return x + u; +} + +weak_alias (__nextupl, nextupl) |