Home | History | Annotate | Line # | Download | only in src 
      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
     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
     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
    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
    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
    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