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

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