sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c

   1 /* Copyright (C) 1995-2016 Free Software Foundation, Inc.
   2    This file is part of the GNU C Library.
   3
   4    The GNU C Library is free software; you can redistribute it and/or
   5    modify it under the terms of the GNU Lesser General Public
   6    License as published by the Free Software Foundation; either
   7    version 2.1 of the License, or (at your option) any later version.
   8
   9    The GNU C Library is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12    Lesser General Public License for more details.
  13
  14    You should have received a copy of the GNU Lesser General Public
  15    License along with the GNU C Library; if not, see
  16    <http://www.gnu.org/licenses/>.  */
  17
  18 #include <ieee754.h>
  19 #include <errno.h>
  20 #include <float.h>
  21 #include <math.h>
  22
  23 /* Need to set this when including gmp headers after system headers.  */
  24 #define HAVE_ALLOCA 1
  25
  26 #include "gmp.h"
  27 #include "gmp-impl.h"
  28
  29 /* Convert a multi-precision integer of the needed number of bits (106
  30    for long double) and an integral power of two to a `long double' in
  31    IBM extended format.  */
  32
  33 long double
  34 __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
  35 {
  36   union ibm_extended_long_double u;
  37   unsigned long lzcount;
  38   unsigned long long hi, lo;
  39   int exponent2;
  40
  41   u.d[0].ieee.negative = sign;
  42   u.d[1].ieee.negative = sign;
  43   u.d[0].ieee.exponent = expt + IEEE754_DOUBLE_BIAS;
  44   u.d[1].ieee.exponent = 0;
  45   exponent2 = expt - 53 + IEEE754_DOUBLE_BIAS;
  46
  47 #if BITS_PER_MP_LIMB == 32
  48   /* The low order 53 bits (52 + hidden) go into the lower double */
  49   lo = frac_ptr[0];
  50   lo |= (frac_ptr[1] & ((1LL << (53 - 32)) - 1)) << 32;
  51   /* The high order 53 bits (52 + hidden) go into the upper double */
  52   hi = (frac_ptr[1] >> (53 - 32)) & ((1 << 11) - 1);
  53   hi |= ((unsigned long long) frac_ptr[2]) << 11;
  54   hi |= ((unsigned long long) frac_ptr[3]) << (32 + 11);
  55 #elif BITS_PER_MP_LIMB == 64
  56   /* The low order 53 bits (52 + hidden) go into the lower double */
  57   lo = frac_ptr[0] & (((mp_limb_t) 1 << 53) - 1);
  58   /* The high order 53 bits (52 + hidden) go into the upper double */
  59   hi = (frac_ptr[0] >> 53) & (((mp_limb_t) 1 << 11) - 1);
  60   hi |= (frac_ptr[1] << 11);
  61 #else
  62   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
  63 #endif
  64
  65   if ((hi & (1LL << 52)) == 0 && (hi | lo) != 0)
  66     {
  67       /* denormal number  */
  68       unsigned long long val = hi ? hi : lo;
  69
  70       if (sizeof (val) == sizeof (long))
  71         lzcount = __builtin_clzl (val);
  72       else if ((val >> 32) != 0)
  73         lzcount = __builtin_clzl ((long) (val >> 32));
  74       else
  75         lzcount = __builtin_clzl ((long) val) + 32;
  76       if (hi)
  77         lzcount = lzcount - (64 - 53);
  78       else
  79         lzcount = lzcount + 53 - (64 - 53);
  80
  81       if (lzcount > u.d[0].ieee.exponent)
  82         {
  83           lzcount = u.d[0].ieee.exponent;
  84           u.d[0].ieee.exponent = 0;
  85           exponent2 -= lzcount;
  86         }
  87       else
  88         {
  89           u.d[0].ieee.exponent -= (lzcount - 1);
  90           exponent2 -= (lzcount - 1);
  91         }
  92
  93       if (lzcount <= 53)
  94         {
  95           hi = (hi << lzcount) | (lo >> (53 - lzcount));
  96           lo = (lo << lzcount) & ((1LL << 53) - 1);
  97         }
  98       else
  99         {
 100           hi = lo << (lzcount - 53);
 101           lo = 0;
 102         }
 103     }
 104
 105   if (lo != 0)
 106     {
 107       /* hidden bit of low double controls rounding of the high double.
 108          If hidden is '1' and either the explicit mantissa is non-zero
 109          or hi is odd, then round up hi and adjust lo (2nd mantissa)
 110          plus change the sign of the low double to compensate.  */
 111       if ((lo & (1LL << 52)) != 0
 112           && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1)) != 0))
 113         {
 114           hi++;
 115           if ((hi & (1LL << 53)) != 0)
 116             {
 117               hi >>= 1;
 118               u.d[0].ieee.exponent++;
 119               if (u.d[0].ieee.exponent == IEEE754_DOUBLE_BIAS + DBL_MAX_EXP)
 120                 {
 121                   /* Overflow.  The appropriate overflowed result must
 122                      be produced (if an infinity, that means the low
 123                      part must be zero).  */
 124                   __set_errno (ERANGE);
 125                   return (sign ? -LDBL_MAX : LDBL_MAX) * LDBL_MAX;
 126                 }
 127             }
 128           u.d[1].ieee.negative = !sign;
 129           lo = (1LL << 53) - lo;
 130         }
 131
 132       /* Normalize the low double.  Shift the mantissa left until
 133          the hidden bit is '1' and adjust the exponent accordingly.  */
 134
 135       if (sizeof (lo) == sizeof (long))
 136         lzcount = __builtin_clzl (lo);
 137       else if ((lo >> 32) != 0)
 138         lzcount = __builtin_clzl ((long) (lo >> 32));
 139       else
 140         lzcount = __builtin_clzl ((long) lo) + 32;
 141       lzcount = lzcount - (64 - 53);
 142       lo <<= lzcount;
 143       exponent2 -= lzcount;
 144
 145       if (exponent2 > 0)
 146         u.d[1].ieee.exponent = exponent2;
 147       else if (exponent2 > -53)
 148         lo >>= 1 - exponent2;
 149       else
 150         lo = 0;
 151     }
 152   else
 153     u.d[1].ieee.negative = 0;
 154
 155   u.d[1].ieee.mantissa1 = lo;
 156   u.d[1].ieee.mantissa0 = lo >> 32;
 157   u.d[0].ieee.mantissa1 = hi;
 158   u.d[0].ieee.mantissa0 = hi >> 32;
 159
 160   return u.ld;
 161 }