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

   1 /* Copyright (C) 1995-2014 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 "gmp.h"
  19 #include "gmp-impl.h"
  20 #include <ieee754.h>
  21 #include <float.h>
  22 #include <math.h>
  23
  24 /* Convert a multi-precision integer of the needed number of bits (106
  25    for long double) and an integral power of two to a `long double' in
  26    IBM extended format.  */
  27
  28 long double
  29 __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
  30 {
  31   union ibm_extended_long_double u;
  32   unsigned long lzcount;
  33   unsigned long long hi, lo;
  34   int exponent2;
  35
  36   u.d[0].ieee.negative = sign;
  37   u.d[1].ieee.negative = sign;
  38   u.d[0].ieee.exponent = expt + IEEE754_DOUBLE_BIAS;
  39   u.d[1].ieee.exponent = 0;
  40   exponent2 = expt - 53 + IEEE754_DOUBLE_BIAS;
  41
  42 #if BITS_PER_MP_LIMB == 32
  43   /* The low order 53 bits (52 + hidden) go into the lower double */
  44   lo = frac_ptr[0];
  45   lo |= (frac_ptr[1] & ((1LL << (53 - 32)) - 1)) << 32;
  46   /* The high order 53 bits (52 + hidden) go into the upper double */
  47   hi = (frac_ptr[1] >> (53 - 32)) & ((1 << 11) - 1);
  48   hi |= ((unsigned long long) frac_ptr[2]) << 11;
  49   hi |= ((unsigned long long) frac_ptr[3]) << (32 + 11);
  50 #elif BITS_PER_MP_LIMB == 64
  51   /* The low order 53 bits (52 + hidden) go into the lower double */
  52   lo = frac_ptr[0] & (((mp_limb_t) 1 << 53) - 1);
  53   /* The high order 53 bits (52 + hidden) go into the upper double */
  54   hi = (frac_ptr[0] >> 53) & (((mp_limb_t) 1 << 11) - 1);
  55   hi |= (frac_ptr[1] << 11);
  56 #else
  57   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
  58 #endif
  59
  60   if ((hi & (1LL << 52)) == 0 && (hi | lo) != 0)
  61     {
  62       /* denormal number  */
  63       unsigned long long val = hi ? hi : lo;
  64
  65       if (sizeof (val) == sizeof (long))
  66         lzcount = __builtin_clzl (val);
  67       else if ((val >> 32) != 0)
  68         lzcount = __builtin_clzl ((long) (val >> 32));
  69       else
  70         lzcount = __builtin_clzl ((long) val) + 32;
  71       if (hi)
  72         lzcount = lzcount - (64 - 53);
  73       else
  74         lzcount = lzcount + 53 - (64 - 53);
  75
  76       if (lzcount > u.d[0].ieee.exponent)
  77         {
  78           lzcount = u.d[0].ieee.exponent;
  79           u.d[0].ieee.exponent = 0;
  80           exponent2 -= lzcount;
  81         }
  82       else
  83         {
  84           u.d[0].ieee.exponent -= (lzcount - 1);
  85           exponent2 -= (lzcount - 1);
  86         }
  87
  88       if (lzcount <= 53)
  89         {
  90           hi = (hi << lzcount) | (lo >> (53 - lzcount));
  91           lo = (lo << lzcount) & ((1LL << 53) - 1);
  92         }
  93       else
  94         {
  95           hi = lo << (lzcount - 53);
  96           lo = 0;
  97         }
  98     }
  99
 100   if (lo != 0)
 101     {
 102       /* hidden bit of low double controls rounding of the high double.
 103          If hidden is '1' and either the explicit mantissa is non-zero
 104          or hi is odd, then round up hi and adjust lo (2nd mantissa)
 105          plus change the sign of the low double to compensate.  */
 106       if ((lo & (1LL << 52)) != 0
 107           && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1)) != 0))
 108         {
 109           hi++;
 110           if ((hi & (1LL << 53)) != 0)
 111             {
 112               hi >>= 1;
 113               u.d[0].ieee.exponent++;
 114             }
 115           u.d[1].ieee.negative = !sign;
 116           lo = (1LL << 53) - lo;
 117         }
 118
 119       /* Normalize the low double.  Shift the mantissa left until
 120          the hidden bit is '1' and adjust the exponent accordingly.  */
 121
 122       if (sizeof (lo) == sizeof (long))
 123         lzcount = __builtin_clzl (lo);
 124       else if ((lo >> 32) != 0)
 125         lzcount = __builtin_clzl ((long) (lo >> 32));
 126       else
 127         lzcount = __builtin_clzl ((long) lo) + 32;
 128       lzcount = lzcount - (64 - 53);
 129       lo <<= lzcount;
 130       exponent2 -= lzcount;
 131
 132       if (exponent2 > 0)
 133         u.d[1].ieee.exponent = exponent2;
 134       else if (exponent2 > -53)
 135         lo >>= 1 - exponent2;
 136       else
 137         lo = 0;
 138     }
 139   else
 140     u.d[1].ieee.negative = 0;
 141
 142   u.d[1].ieee.mantissa1 = lo;
 143   u.d[1].ieee.mantissa0 = lo >> 32;
 144   u.d[0].ieee.mantissa1 = hi;
 145   u.d[0].ieee.mantissa0 = hi >> 32;
 146
 147   return u.ld;
 148 }