Back to index

glibc  2.9
mpn2ldbl.c
Go to the documentation of this file.
00001 /* Copyright (C) 1995, 1996, 1997, 1998, 1999, 2002, 2003, 2004, 2006, 2007
00002        Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004 
00005    The GNU C Library is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU Lesser General Public
00007    License as published by the Free Software Foundation; either
00008    version 2.1 of the License, or (at your option) any later version.
00009 
00010    The GNU C Library is distributed in the hope that it will be useful,
00011    but WITHOUT ANY WARRANTY; without even the implied warranty of
00012    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013    Lesser General Public License for more details.
00014 
00015    You should have received a copy of the GNU Lesser General Public
00016    License along with the GNU C Library; if not, write to the Free
00017    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00018    02111-1307 USA.  */
00019 
00020 #include "gmp.h"
00021 #include "gmp-impl.h"
00022 #include <ieee754.h>
00023 #include <float.h>
00024 #include <math.h>
00025 
00026 /* Convert a multi-precision integer of the needed number of bits (106
00027    for long double) and an integral power of two to a `long double' in
00028    IBM extended format.  */
00029 
00030 long double
00031 __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
00032 {
00033   union ibm_extended_long_double u;
00034   unsigned long lzcount;
00035   unsigned long long hi, lo;
00036   int exponent2;
00037 
00038   u.ieee.negative = sign;
00039   u.ieee.negative2 = sign;
00040   u.ieee.exponent = expt + IBM_EXTENDED_LONG_DOUBLE_BIAS;
00041   u.ieee.exponent2 = 0;
00042   exponent2 = expt - 53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
00043 
00044 #if BITS_PER_MP_LIMB == 32
00045   /* The low order 53 bits (52 + hidden) go into the lower double */
00046   lo = frac_ptr[0];
00047   lo |= (frac_ptr[1] & ((1LL << (53 - 32)) - 1)) << 32;
00048   /* The high order 53 bits (52 + hidden) go into the upper double */
00049   hi = (frac_ptr[1] >> (53 - 32)) & ((1 << 11) - 1);
00050   hi |= ((unsigned long long) frac_ptr[2]) << 11;
00051   hi |= ((unsigned long long) frac_ptr[3]) << (32 + 11);
00052 #elif BITS_PER_MP_LIMB == 64
00053   /* The low order 53 bits (52 + hidden) go into the lower double */
00054   lo = frac_ptr[0] & (((mp_limb_t) 1 << 53) - 1);
00055   /* The high order 53 bits (52 + hidden) go into the upper double */
00056   hi = (frac_ptr[0] >> 53) & (((mp_limb_t) 1 << 11) - 1);
00057   hi |= (frac_ptr[1] << 11);
00058 #else
00059   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
00060 #endif
00061 
00062   if ((hi & (1LL << 52)) == 0 && (hi | lo) != 0)
00063     {
00064       /* denormal number  */
00065       unsigned long long val = hi ? hi : lo;
00066 
00067       if (sizeof (val) == sizeof (long))
00068        lzcount = __builtin_clzl (val);
00069       else if ((val >> 32) != 0)
00070        lzcount = __builtin_clzl ((long) (val >> 32));
00071       else
00072        lzcount = __builtin_clzl ((long) val) + 32;
00073       if (hi)
00074        lzcount = lzcount - 11;
00075       else
00076        lzcount = lzcount + 42;
00077 
00078       if (lzcount > u.ieee.exponent)
00079        {
00080          lzcount = u.ieee.exponent;
00081          u.ieee.exponent = 0;
00082          exponent2 -= lzcount;
00083        }
00084       else
00085        {
00086          u.ieee.exponent -= (lzcount - 1);
00087          exponent2 -= (lzcount - 1);
00088        }
00089 
00090       if (lzcount <= 53)
00091        {
00092          hi = (hi << lzcount) | (lo >> (53 - lzcount));
00093          lo = (lo << lzcount) & ((1LL << 53) - 1);
00094        }
00095       else
00096        {
00097          hi = lo << (lzcount - 53);
00098          lo = 0;
00099        }
00100     }
00101 
00102   if (lo != 0L)
00103     {
00104       /* hidden2 bit of low double controls rounding of the high double.
00105         If hidden2 is '1' and either the explicit mantissa is non-zero
00106         or hi is odd, then round up hi and adjust lo (2nd mantissa)
00107         plus change the sign of the low double to compensate.  */
00108       if ((lo & (1LL << 52)) != 0
00109          && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1))))
00110        {
00111          hi++;
00112          if ((hi & ((1LL << 52) - 1)) == 0)
00113            {
00114              if ((hi & (1LL << 53)) != 0)
00115               hi -= 1LL << 52;
00116              u.ieee.exponent++;
00117            }
00118          u.ieee.negative2 = !sign;
00119          lo = (1LL << 53) - lo;
00120        }
00121 
00122       /* The hidden bit of the lo mantissa is zero so we need to normalize
00123         it for the low double.  Shift it left until the hidden bit is '1'
00124         then adjust the 2nd exponent accordingly.  */
00125 
00126       if (sizeof (lo) == sizeof (long))
00127        lzcount = __builtin_clzl (lo);
00128       else if ((lo >> 32) != 0)
00129        lzcount = __builtin_clzl ((long) (lo >> 32));
00130       else
00131        lzcount = __builtin_clzl ((long) lo) + 32;
00132       lzcount = lzcount - 11;
00133       if (lzcount > 0)
00134        {
00135          lo = lo << lzcount;
00136          exponent2 = exponent2 - lzcount;
00137        }
00138       if (exponent2 > 0)
00139        u.ieee.exponent2 = exponent2;
00140       else
00141        lo >>= 1 - exponent2;
00142     }
00143   else
00144     u.ieee.negative2 = 0;
00145 
00146   u.ieee.mantissa3 = lo & 0xffffffffLL;
00147   u.ieee.mantissa2 = (lo >> 32) & 0xfffff;
00148   u.ieee.mantissa1 = hi & 0xffffffffLL;
00149   u.ieee.mantissa0 = (hi >> 32) & ((1LL << (LDBL_MANT_DIG - 86)) - 1);
00150 
00151   return u.d;
00152 }