libquadmath/printf/flt1282mpn.c

   1 /* Copyright (C) 1995,1996,1997,1998,1999,2002,2003
   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 <float.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include "gmp-impl.h"
  24
  25 /* Convert a `__float128' in IEEE854 quad-precision format to a
  26    multi-precision integer representing the significand scaled up by its
  27    number of bits (113 for long double) and an integral power of two
  28    (MPN frexpl). */
  29
  30 mp_size_t
  31 mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size,
  32                       int *expt, int *is_neg,
  33                       __float128 value)
  34 {
  35   ieee854_float128 u;
  36   u.value = value;
  37
  38   *is_neg = u.ieee.negative;
  39   *expt = (int) u.ieee.exponent - IEEE854_FLOAT128_BIAS;
  40
  41 #if BITS_PER_MP_LIMB == 32
  42   res_ptr[0] = u.ieee.mant_low; /* Low-order 32 bits of fraction.  */
  43   res_ptr[1] = (u.ieee.mant_low >> 32);
  44   res_ptr[2] = u.ieee.mant_high;
  45   res_ptr[3] = (u.ieee.mant_high >> 32); /* High-order 32 bits.  */
  46   #define N 4
  47 #elif BITS_PER_MP_LIMB == 64
  48   res_ptr[0] = u.ieee.mant_low;
  49   res_ptr[1] = u.ieee.mant_high;
  50   #define N 2
  51 #else
  52   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
  53 #endif
  54 /* The format does not fill the last limb.  There are some zeros.  */
  55 #define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
  56                            - (FLT128_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
  57
  58   if (u.ieee.exponent == 0)
  59     {
  60       /* A biased exponent of zero is a special case.
  61          Either it is a zero or it is a denormal number.  */
  62       if (res_ptr[0] == 0 && res_ptr[1] == 0
  63           && res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4.  */
  64         /* It's zero.  */
  65         *expt = 0;
  66       else
  67         {
  68           /* It is a denormal number, meaning it has no implicit leading
  69              one bit, and its exponent is in fact the format minimum.  */
  70           int cnt;
  71
  72 #if N == 2
  73           if (res_ptr[N - 1] != 0)
  74             {
  75               count_leading_zeros (cnt, res_ptr[N - 1]);
  76               cnt -= NUM_LEADING_ZEROS;
  77               res_ptr[N - 1] = res_ptr[N - 1] << cnt
  78                                | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
  79               res_ptr[0] <<= cnt;
  80               *expt = FLT128_MIN_EXP - 1 - cnt;
  81             }
  82           else
  83             {
  84               count_leading_zeros (cnt, res_ptr[0]);
  85               if (cnt >= NUM_LEADING_ZEROS)
  86                 {
  87                   res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS);
  88                   res_ptr[0] = 0;
  89                 }
  90               else
  91                 {
  92                   res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt);
  93                   res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt);
  94                 }
  95               *expt = FLT128_MIN_EXP - 1
  96                 - (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt;
  97             }
  98 #else
  99           int j, k, l;
 100
 101           for (j = N - 1; j > 0; j--)
 102             if (res_ptr[j] != 0)
 103               break;
 104
 105           count_leading_zeros (cnt, res_ptr[j]);
 106           cnt -= NUM_LEADING_ZEROS;
 107           l = N - 1 - j;
 108           if (cnt < 0)
 109             {
 110               cnt += BITS_PER_MP_LIMB;
 111               l--;
 112             }
 113           if (!cnt)
 114             for (k = N - 1; k >= l; k--)
 115               res_ptr[k] = res_ptr[k-l];
 116           else
 117             {
 118               for (k = N - 1; k > l; k--)
 119                 res_ptr[k] = res_ptr[k-l] << cnt
 120                              | res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt);
 121               res_ptr[k--] = res_ptr[0] << cnt;
 122             }
 123
 124           for (; k >= 0; k--)
 125             res_ptr[k] = 0;
 126           *expt = FLT128_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt;
 127 #endif
 128         }
 129     }
 130   else
 131     /* Add the implicit leading one bit for a normalized number.  */
 132     res_ptr[N - 1] |= (mp_limb_t) 1 << (FLT128_MANT_DIG - 1
 133                                         - ((N - 1) * BITS_PER_MP_LIMB));
 134
 135   return N;
 136 }