1/* mpfr_nextabove, mpfr_nextbelow, mpfr_nexttoward -- next representable
2floating-point number
3
4Copyright 1999, 2001-2023 Free Software Foundation, Inc.
5Contributed by the AriC and Caramba projects, INRIA.
6
7This file is part of the GNU MPFR Library.
8
9The GNU MPFR Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14The GNU MPFR Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
21https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2251 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
46void
47mpfr_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
88void
89mpfr_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
120void
121mpfr_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
134void
135mpfr_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
149void
150mpfr_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