Home | History | Annotate | Line # | Download | only in src 
      1 /* mpfr_get_float128 -- convert a multiple precision floating-point
      2                         number to a _Float128 number
      3 
      4 Copyright 2012-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 #include "mpfr-impl.h"
     25 
     26 #ifdef MPFR_WANT_FLOAT128
     27 
     28 /* generic code */
     29 _Float128
     30 mpfr_get_float128 (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
     31 {
     32 
     33   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     34     return (_Float128) mpfr_get_d (x, rnd_mode);
     35   else /* now x is a normal non-zero number */
     36     {
     37       _Float128 r; /* result */
     38       _Float128 m;
     39       mpfr_exp_t e;  /* exponent of x (before rounding) */
     40       mpfr_exp_t sh; /* exponent shift, so that x/2^sh is in the double range */
     41       const int emin = -16381;
     42       const int esub = emin - IEEE_FLOAT128_MANT_DIG;
     43       int sign;
     44 
     45       sign = MPFR_SIGN (x);
     46       e = MPFR_GET_EXP (x);
     47 
     48       if (MPFR_UNLIKELY (e <= esub))
     49         {
     50           if (MPFR_IS_LIKE_RNDZ (rnd_mode, sign < 0) ||
     51               (rnd_mode == MPFR_RNDN && (e < esub || mpfr_powerof2_raw (x))))
     52             return sign < 0 ? -0.0 : 0.0;
     53           r = 1.0;
     54           sh = esub;
     55         }
     56       else
     57         {
     58           mpfr_t y;
     59           mp_limb_t *yp;
     60           int prec, i;  /* small enough to fit in an int */
     61           MPFR_SAVE_EXPO_DECL (expo);
     62 
     63           MPFR_SAVE_EXPO_MARK (expo);
     64 
     65           /* First round x to the target _Float128 precision, taking the
     66              reduced precision of the subnormals into account, so that all
     67              subsequent operations are exact (this avoids double rounding
     68              problems). */
     69           prec = e < emin ? e - esub : IEEE_FLOAT128_MANT_DIG;
     70           MPFR_ASSERTD (prec >= MPFR_PREC_MIN);
     71           mpfr_init2 (y, prec);
     72 
     73           mpfr_set (y, x, rnd_mode);
     74           sh = MPFR_GET_EXP (y);
     75           MPFR_SET_EXP (y, 0);
     76           MPFR_SET_POS (y);
     77           yp = MPFR_MANT (y);
     78 
     79           r = 0.0;
     80           for (i = 0; i < MPFR_LIMB_SIZE (y); i++)
     81             {
     82               /* Note: MPFR_LIMB_MAX is avoided below as it might not
     83                  always work if GMP_NUMB_BITS > IEEE_FLOAT128_MANT_DIG.
     84                  MPFR_LIMB_HIGHBIT has the advantage to fit on 1 bit. */
     85               r += yp[i];
     86               r *= 1 / (2 * (_Float128) MPFR_LIMB_HIGHBIT);
     87             }
     88 
     89           mpfr_clear (y);
     90 
     91           MPFR_SAVE_EXPO_FREE (expo);
     92         }
     93 
     94       /* we now have to multiply r by 2^sh */
     95       MPFR_ASSERTD (r > 0);
     96       if (sh != 0)
     97         {
     98           /* An overflow may occur (example: 0.5*2^1024) */
     99           while (r < 1.0)
    100             {
    101               r += r;
    102               sh--;
    103             }
    104 
    105           if (sh > 0)
    106             m = 2.0;
    107           else
    108             {
    109               m = 0.5;
    110               sh = -sh;
    111             }
    112 
    113           for (;;)
    114             {
    115               if (sh % 2)
    116                 r = r * m;
    117               sh >>= 1;
    118               if (sh == 0)
    119                 break;
    120               m = m * m;
    121             }
    122         }
    123       if (sign < 0)
    124         r = -r;
    125       return r;
    126     }
    127 }
    128 
    129 #endif /* MPFR_WANT_FLOAT128 */
    130