Home | History | Annotate | Line # | Download | only in printf 
      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.mantissa3; /* Low-order 32 bits of fraction.  */
     43   res_ptr[1] = u.ieee.mantissa2;
     44   res_ptr[2] = u.ieee.mantissa1;
     45   res_ptr[3] = u.ieee.mantissa0; /* High-order 32 bits.  */
     46   #define N 4
     47 #elif BITS_PER_MP_LIMB == 64
     48   res_ptr[0] = ((mp_limb_t) u.ieee.mantissa2 << 32) | u.ieee.mantissa3;
     49   res_ptr[1] = ((mp_limb_t) u.ieee.mantissa0 << 32) | u.ieee.mantissa1;
     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 }
    137