sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.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 "longlong.h"
  21 #include <ieee754.h>
  22 #include <float.h>
  23 #include <math.h>
  24 #include <stdlib.h>
  25
  26 /* Convert a `long double' in IBM extended format to a multi-precision
  27    integer representing the significand scaled up by its number of
  28    bits (106 for long double) and an integral power of two (MPN
  29    frexpl). */
  30
  31 mp_size_t
  32 __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
  33                            int *expt, int *is_neg,
  34                            long double value)
  35 {
  36   union ibm_extended_long_double u;
  37   unsigned long long hi, lo;
  38   int ediff;
  39
  40   u.ld = value;
  41
  42   *is_neg = u.d[0].ieee.negative;
  43   *expt = (int) u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;
  44
  45   lo = ((long long) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
  46   hi = ((long long) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
  47
  48   /* If the lower double is not a denormal or zero then set the hidden
  49      53rd bit.  */
  50   if (u.d[1].ieee.exponent != 0)
  51     lo |= 1ULL << 52;
  52   else
  53     lo = lo << 1;
  54
  55   /* The lower double is normalized separately from the upper.  We may
  56      need to adjust the lower manitissa to reflect this.  */
  57   ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;
  58   if (ediff > 0)
  59     {
  60       if (ediff < 64)
  61         lo = lo >> ediff;
  62       else
  63         lo = 0;
  64     }
  65   else if (ediff < 0)
  66     lo = lo << -ediff;
  67
  68   /* The high double may be rounded and the low double reflects the
  69      difference between the long double and the rounded high double
  70      value.  This is indicated by a differnce between the signs of the
  71      high and low doubles.  */
  72   if (u.d[0].ieee.negative != u.d[1].ieee.negative
  73       && lo != 0)
  74     {
  75       lo = (1ULL << 53) - lo;
  76       if (hi == 0)
  77         {
  78           /* we have a borrow from the hidden bit, so shift left 1.  */
  79           hi = 0x0ffffffffffffeLL | (lo >> 51);
  80           lo = 0x1fffffffffffffLL & (lo << 1);
  81           (*expt)--;
  82         }
  83       else
  84         hi--;
  85     }
  86 #if BITS_PER_MP_LIMB == 32
  87   /* Combine the mantissas to be contiguous.  */
  88   res_ptr[0] = lo;
  89   res_ptr[1] = (hi << (53 - 32)) | (lo >> 32);
  90   res_ptr[2] = hi >> 11;
  91   res_ptr[3] = hi >> (32 + 11);
  92   #define N 4
  93 #elif BITS_PER_MP_LIMB == 64
  94   /* Combine the two mantissas to be contiguous.  */
  95   res_ptr[0] = (hi << 53) | lo;
  96   res_ptr[1] = hi >> 11;
  97   #define N 2
  98 #else
  99   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
 100 #endif
 101 /* The format does not fill the last limb.  There are some zeros.  */
 102 #define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
 103                            - (LDBL_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
 104
 105   if (u.d[0].ieee.exponent == 0)
 106     {
 107       /* A biased exponent of zero is a special case.
 108          Either it is a zero or it is a denormal number.  */
 109       if (res_ptr[0] == 0 && res_ptr[1] == 0
 110           && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4.  */
 111         /* It's zero.  */
 112         *expt = 0;
 113       else
 114         {
 115           /* It is a denormal number, meaning it has no implicit leading
 116              one bit, and its exponent is in fact the format minimum.  We
 117              use DBL_MIN_EXP instead of LDBL_MIN_EXP below because the
 118              latter describes the properties of both parts together, but
 119              the exponent is computed from the high part only.  */
 120           int cnt;
 121
 122 #if N == 2
 123           if (res_ptr[N - 1] != 0)
 124             {
 125               count_leading_zeros (cnt, res_ptr[N - 1]);
 126               cnt -= NUM_LEADING_ZEROS;
 127               res_ptr[N - 1] = res_ptr[N - 1] << cnt
 128                                | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
 129               res_ptr[0] <<= cnt;
 130               *expt = DBL_MIN_EXP - 1 - cnt;
 131             }
 132           else
 133             {
 134               count_leading_zeros (cnt, res_ptr[0]);
 135               if (cnt >= NUM_LEADING_ZEROS)
 136                 {
 137                   res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS);
 138                   res_ptr[0] = 0;
 139                 }
 140               else
 141                 {
 142                   res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt);
 143                   res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt);
 144                 }
 145               *expt = DBL_MIN_EXP - 1
 146                 - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt;
 147             }
 148 #else
 149           int j, k, l;
 150
 151           for (j = N - 1; j > 0; j--)
 152             if (res_ptr[j] != 0)
 153               break;
 154
 155           count_leading_zeros (cnt, res_ptr[j]);
 156           cnt -= NUM_LEADING_ZEROS;
 157           l = N - 1 - j;
 158           if (cnt < 0)
 159             {
 160               cnt += BITS_PER_MP_LIMB;
 161               l--;
 162             }
 163           if (!cnt)
 164             for (k = N - 1; k >= l; k--)
 165               res_ptr[k] = res_ptr[k-l];
 166           else
 167             {
 168               for (k = N - 1; k > l; k--)
 169                 res_ptr[k] = res_ptr[k-l] << cnt
 170                              | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt);
 171               res_ptr[k--] = res_ptr[0] << cnt;
 172             }
 173
 174           for (; k >= 0; k--)
 175             res_ptr[k] = 0;
 176           *expt = DBL_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt;
 177 #endif
 178         }
 179     }
 180   else
 181     /* Add the implicit leading one bit for a normalized number.  */
 182     res_ptr[N - 1] |= (mp_limb_t) 1 << (LDBL_MANT_DIG - 1
 183                                         - ((N - 1) * BITS_PER_MP_LIMB));
 184
 185   return N;
 186 }