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

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