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 */