1 /* mpfr_nextabove, mpfr_nextbelow, mpfr_nexttoward -- next representable
2 floating-point number
3
4 Copyright 1999, 2001-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba projects, INRIA.
6
7 This file is part of the GNU MPFR Library.
8
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23
24 /* Note concerning the exceptions: In case of NaN result, the NaN flag is
25 * set as usual. No underflow or overflow is generated (this contradicts
26 * the obsolete IEEE 754-1985 standard for nextafter, but conforms to the
27 * IEEE 754-2008/2019 standards, where the nextUp and nextDown operations
28 * replaced nextafter).
29 *
30 * For mpfr_nexttoward where x and y are zeros of different signs, the
31 * behavior is the same as the nextafter function from IEEE 754-1985
32 * (x is returned), but differs from the ISO C99 nextafter and nexttoward
33 * functions (where y is returned).
34 *
35 * In short:
36 * - mpfr_nextabove and mpfr_nextbelow are similar to nextUp and nextDown
37 * respectively (IEEE 2008+, ISO C2x), but with the usual MPFR rule for
38 * the NaN flag (because MPFR has a single kind of NaN);
39 * - mpfr_nexttoward is also similar to these functions concerning the
40 * exceptions and the sign of 0, and for the particular case x = y, it
41 * follows the old nextafter function from IEEE 754-1985.
42 */
43
44 #include "mpfr-impl.h"
45
46 void
mpfr_nexttozero(mpfr_ptr x)47 mpfr_nexttozero (mpfr_ptr x)
48 {
49 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
50 {
51 if (MPFR_IS_INF (x))
52 {
53 mpfr_setmax (x, __gmpfr_emax);
54 }
55 else
56 {
57 MPFR_ASSERTN (MPFR_IS_ZERO (x));
58 MPFR_CHANGE_SIGN(x);
59 mpfr_setmin (x, __gmpfr_emin);
60 }
61 }
62 else
63 {
64 mp_size_t xn;
65 int sh;
66 mp_limb_t *xp;
67
68 xn = MPFR_LIMB_SIZE (x);
69 MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC(x));
70 xp = MPFR_MANT(x);
71 mpn_sub_1 (xp, xp, xn, MPFR_LIMB_ONE << sh);
72 if (MPFR_UNLIKELY (MPFR_LIMB_MSB (xp[xn-1]) == 0))
73 { /* was an exact power of two: not normalized any more,
74 thus do not use MPFR_GET_EXP. */
75 mpfr_exp_t exp = MPFR_EXP (x);
76 if (MPFR_UNLIKELY (exp == __gmpfr_emin))
77 MPFR_SET_ZERO(x);
78 else
79 {
80 MPFR_SET_EXP (x, exp - 1);
81 /* The following is valid whether xn = 1 or xn > 1. */
82 xp[xn-1] |= MPFR_LIMB_HIGHBIT;
83 }
84 }
85 }
86 }
87
88 void
mpfr_nexttoinf(mpfr_ptr x)89 mpfr_nexttoinf (mpfr_ptr x)
90 {
91 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
92 {
93 if (MPFR_IS_ZERO (x))
94 mpfr_setmin (x, __gmpfr_emin);
95 }
96 else
97 {
98 mp_size_t xn;
99 int sh;
100 mp_limb_t *xp;
101
102 xn = MPFR_LIMB_SIZE (x);
103 MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC(x));
104 xp = MPFR_MANT(x);
105 if (MPFR_UNLIKELY (mpn_add_1 (xp, xp, xn, MPFR_LIMB_ONE << sh)))
106 /* got 1.0000... */
107 {
108 mpfr_exp_t exp = MPFR_EXP (x);
109 if (MPFR_UNLIKELY (exp == __gmpfr_emax))
110 MPFR_SET_INF (x);
111 else
112 {
113 MPFR_SET_EXP (x, exp + 1);
114 xp[xn-1] = MPFR_LIMB_HIGHBIT;
115 }
116 }
117 }
118 }
119
120 void
mpfr_nextabove(mpfr_ptr x)121 mpfr_nextabove (mpfr_ptr x)
122 {
123 if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
124 {
125 __gmpfr_flags |= MPFR_FLAGS_NAN;
126 return;
127 }
128 if (MPFR_IS_NEG(x))
129 mpfr_nexttozero (x);
130 else
131 mpfr_nexttoinf (x);
132 }
133
134 void
mpfr_nextbelow(mpfr_ptr x)135 mpfr_nextbelow (mpfr_ptr x)
136 {
137 if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
138 {
139 __gmpfr_flags |= MPFR_FLAGS_NAN;
140 return;
141 }
142
143 if (MPFR_IS_NEG(x))
144 mpfr_nexttoinf (x);
145 else
146 mpfr_nexttozero (x);
147 }
148
149 void
mpfr_nexttoward(mpfr_ptr x,mpfr_srcptr y)150 mpfr_nexttoward (mpfr_ptr x, mpfr_srcptr y)
151 {
152 int s;
153
154 if (MPFR_UNLIKELY(MPFR_IS_NAN(x)))
155 {
156 __gmpfr_flags |= MPFR_FLAGS_NAN;
157 return;
158 }
159 else if (MPFR_UNLIKELY(MPFR_IS_NAN(x) || MPFR_IS_NAN(y)))
160 {
161 MPFR_SET_NAN(x);
162 __gmpfr_flags |= MPFR_FLAGS_NAN;
163 return;
164 }
165
166 s = mpfr_cmp (x, y);
167 if (s == 0)
168 return;
169 else if (s < 0)
170 mpfr_nextabove (x);
171 else
172 mpfr_nextbelow (x);
173 }
174