Home | History | Annotate | Line # | Download | only in soft-fp
      1 /* Software floating-point emulation.
      2    Convert a _BitInt to _Decimal128.
      3 
      4    Copyright (C) 2023 Free Software Foundation, Inc.
      5 
      6 This file is part of GCC.
      7 
      8 GCC is free software; you can redistribute it and/or modify it under
      9 the terms of the GNU General Public License as published by the Free
     10 Software Foundation; either version 3, or (at your option) any later
     11 version.
     12 
     13 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
     14 WARRANTY; without even the implied warranty of MERCHANTABILITY or
     15 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     16 for more details.
     17 
     18 Under Section 7 of GPL version 3, you are granted additional
     19 permissions described in the GCC Runtime Library Exception, version
     20 3.1, as published by the Free Software Foundation.
     21 
     22 You should have received a copy of the GNU General Public License and
     23 a copy of the GCC Runtime Library Exception along with this program;
     24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
     25 <http://www.gnu.org/licenses/>.  */
     26 
     27 #include "soft-fp.h"
     28 #include "bitint.h"
     29 
     30 #ifdef __BITINT_MAXWIDTH__
     31 extern _Decimal128 __bid_floatbitinttd (const UBILtype *, SItype);
     32 
     33 _Decimal128
     34 __bid_floatbitinttd (const UBILtype *i, SItype iprec)
     35 {
     36   iprec = bitint_reduce_prec (&i, iprec);
     37   USItype aiprec = iprec < 0 ? -iprec : iprec;
     38   USItype in = (aiprec + BIL_TYPE_SIZE - 1) / BIL_TYPE_SIZE;
     39   USItype idx = BITINT_END (0, in - 1);
     40   UBILtype msb = i[idx];
     41   UDItype mantissahi, mantissalo;
     42   SItype exponent = 0;
     43   UBILtype inexact = 0;
     44   union { _Decimal128 d; UDItype u[2]; } u, ui;
     45   if (aiprec % BIL_TYPE_SIZE)
     46     {
     47       if (iprec > 0)
     48 	msb &= ((UBILtype) 1 << (aiprec % BIL_TYPE_SIZE)) - 1;
     49       else
     50 	msb |= (UBILtype) -1 << (aiprec % BIL_TYPE_SIZE);
     51     }
     52   if (iprec < 0)
     53     {
     54       SItype n = sizeof (0ULL) * __CHAR_BIT__ + 1 - __builtin_clzll (~msb);
     55       aiprec = (in - 1) * BIL_TYPE_SIZE + n;
     56     }
     57   else if (msb == 0)
     58     aiprec = 1;
     59   else
     60     {
     61       SItype n = sizeof (0ULL) * __CHAR_BIT__ - __builtin_clzll (msb);
     62       aiprec = (in - 1) * BIL_TYPE_SIZE + n;
     63     }
     64   /* Number of bits in
     65      (_BitInt(32768)) 9999999999999999999999999999999999e+6111DL.  */
     66   if (aiprec > 20414 + (iprec < 0))
     67     {
     68     ovf:
     69       if (iprec < 0)
     70 	u.d = -9000000000000000000000000000000000e+6111DL;
     71       else
     72 	u.d = 9000000000000000000000000000000000e+6111DL;
     73       __asm ("" : "+g" (u.d));
     74       u.d += u.d;
     75       __asm ("" : "+g" (u.d));
     76       goto done;
     77     }
     78   /* Bit precision of 9999999999999999999999999999999999uwb.  */
     79   if (aiprec >= 113)
     80     {
     81       USItype pow10_limbs, q_limbs, q2_limbs, j, k;
     82       USItype exp_bits = 0, e;
     83       UBILtype *buf;
     84       /* First do a possibly large divide smaller enough such that
     85 	 we only need to check remainder for 0 or non-0 and then
     86 	 we'll do further division.  */
     87       if (aiprec >= 113 + 4 + 10)
     88 	{
     89 	  exp_bits = ((aiprec - 113 - 4) * (UDItype) 30) / 299;
     90 	  exponent = exp_bits * 3;
     91 	  /* Upper estimate for pow10 (exponent) bits.  */
     92 	  exp_bits = exp_bits * 10 - exp_bits / 30;
     93 	}
     94       pow10_limbs = (exp_bits + BIL_TYPE_SIZE - 1) / BIL_TYPE_SIZE;
     95       /* 127 is the highest number of quotient bits needed on
     96 	 aiprec range of [127, 20414].  E.g. if aiprec is 20409,
     97 	 exponent will be 6105 and exp_bits 20283.  20409 - 20283 + 1
     98 	 is 127.  */
     99       q_limbs = (127 + BIL_TYPE_SIZE - 1) / BIL_TYPE_SIZE;
    100       q2_limbs = 128 / BIL_TYPE_SIZE;
    101       buf = __builtin_alloca ((q_limbs + pow10_limbs * 2 + q2_limbs + 2)
    102 			      * sizeof (UBILtype));
    103       if (exponent)
    104 	{
    105 	  __bid_pow10bitint (buf + q_limbs, exp_bits, exponent);
    106 	  __divmodbitint4 (buf, q_limbs * BIL_TYPE_SIZE,
    107 			   buf + q_limbs + pow10_limbs,
    108 			   pow10_limbs * BIL_TYPE_SIZE,
    109 			   i, iprec < 0 ? -aiprec : aiprec,
    110 			   buf + q_limbs, exp_bits);
    111 	  if (iprec < 0)
    112 	    bitint_negate (buf + BITINT_END (q_limbs - 1, 0),
    113 			   buf + BITINT_END (q_limbs - 1, 0), q_limbs);
    114 	  inexact = buf[q_limbs + pow10_limbs];
    115 	  for (j = 1; j < pow10_limbs; ++j)
    116 	    inexact |= buf[q_limbs + pow10_limbs + j];
    117 	}
    118       else
    119 	{
    120 	  __builtin_memcpy (buf + BITINT_END (q_limbs - in + 1, 0), i,
    121 			    (in - 1) * sizeof (UBILtype));
    122 	  buf[BITINT_END (q_limbs - in, in - 1)] = msb;
    123 	  if (iprec < 0)
    124 	    bitint_negate (buf + BITINT_END (q_limbs - 1, 0),
    125 			   buf + BITINT_END (q_limbs - 1, 0), in);
    126 	  if (q_limbs > in)
    127 	    __builtin_memset (buf + BITINT_END (0, in), '\0',
    128 			      (q_limbs - in) * sizeof (UBILtype));
    129 	}
    130       e = 0;
    131       for (j = 3; j; )
    132 	{
    133 	  USItype eprev = e;
    134 	  __bid_pow10bitint (buf + q_limbs + pow10_limbs * 2 + 1,
    135 			     128, 33 + e + j);
    136 	  for (k = BITINT_END (0, q_limbs - 1);
    137 	       k != BITINT_END (q_limbs, (USItype) -1); k -= BITINT_INC)
    138 	    if (buf[k] > buf[q_limbs + pow10_limbs * 2 + 1 + k])
    139 	      {
    140 		e += j;
    141 		break;
    142 	      }
    143 	    else if (buf[k] < buf[q_limbs + pow10_limbs * 2 + 1 + k])
    144 	      break;
    145 	  if (k == BITINT_END (q_limbs, (USItype) -1))
    146 	    e += j;
    147 	  if (j == 2 && e != eprev)
    148 	    break;
    149 	  else
    150 	    --j;
    151 	}
    152       exponent += e;
    153       if (exponent > 6111)
    154 	goto ovf;
    155       if (e)
    156 	{
    157 	  UBILtype rem, half;
    158 	  __bid_pow10bitint (buf + q_limbs + pow10_limbs * 2,
    159 			     BIL_TYPE_SIZE, e);
    160 	  __divmodbitint4 (buf + q_limbs + pow10_limbs * 2 + 1,
    161 			   q2_limbs * BIL_TYPE_SIZE,
    162 			   buf + q_limbs + pow10_limbs * 2 + 1 + q2_limbs,
    163 			   BIL_TYPE_SIZE,
    164 			   buf, q_limbs * BIL_TYPE_SIZE,
    165 			   buf + q_limbs + pow10_limbs * 2, BIL_TYPE_SIZE);
    166 	  half = buf[q_limbs + pow10_limbs * 2] / 2;
    167 	  rem = buf[q_limbs + pow10_limbs * 2 + 1 + q2_limbs];
    168 	  if (inexact)
    169 	    {
    170 	      /* If first division discovered some non-0 digits
    171 		 and this second division is by 10, e.g.
    172 		 for XXXXXX5499999999999 or XXXXXX5000000000001
    173 		 if first division is by 10^12 and second by 10^1,
    174 		 doing rem |= 1 wouldn't change the 5.  Similarly
    175 		 for rem 4 doing rem |= 1 would change it to 5,
    176 		 but we don't want to change it in that case.  */
    177 	      if (e == 1)
    178 		{
    179 		  if (rem == 5)
    180 		    rem = 6;
    181 		  else if (rem != 4)
    182 		    rem |= 1;
    183 		}
    184 	      else
    185 		rem |= 1;
    186 	    }
    187 	  /* Set inexact to 0, 1, 2, 3 depending on if remainder
    188 	     of the divisions is exact 0, smaller than 10^exponent / 2,
    189 	     exactly 10^exponent / 2 or greater than that.  */
    190 	  if (rem >= half)
    191 	    inexact = 2 + (rem > half);
    192 	  else
    193 	    inexact = (rem != 0);
    194 #if BIL_TYPE_SIZE == 64
    195 	  mantissahi = buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (0, 1)];
    196 	  mantissalo = buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (1, 0)];
    197 #else
    198 	  mantissahi
    199 	    = ((UDItype)
    200 	       buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (0, 3)] << 32
    201 	       | buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (1, 2)]);
    202 	  mantissalo
    203 	    = ((UDItype)
    204 	       buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (2, 1)] << 32
    205 	       | buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (3, 0)]);
    206 #endif
    207 	}
    208       else
    209 	{
    210 #if BIL_TYPE_SIZE == 64
    211 	  mantissahi = buf[BITINT_END (0, 1)];
    212 	  mantissalo = buf[BITINT_END (1, 0)];
    213 #else
    214 	  mantissahi = ((UDItype) buf[BITINT_END (0, 3)] << 32
    215 			| buf[BITINT_END (1, 2)]);
    216 	  mantissalo = ((UDItype) buf[BITINT_END (2, 1)] << 32
    217 			| buf[BITINT_END (3, 0)]);
    218 #endif
    219 	}
    220     }
    221   else
    222     {
    223       mantissahi = iprec < 0 ? -1 : 0;
    224 #if BIL_TYPE_SIZE == 64
    225       if (in == 1)
    226 	mantissalo = msb;
    227       else
    228 	{
    229 	  mantissahi = msb;
    230 	  mantissalo = i[BITINT_END (1, 0)];
    231 	}
    232 #else
    233       if (in <= 2)
    234 	{
    235 	  if (in == 1)
    236 	    mantissalo = iprec < 0 ? (UDItype) (BILtype) msb : (UDItype) msb;
    237 	  else
    238 	    mantissalo = (UDItype) msb << 32 | i[BITINT_END (1, 0)];
    239 	}
    240       else
    241 	{
    242 	  if (in == 3)
    243 	    mantissahi = iprec < 0 ? (UDItype) (BILtype) msb : (UDItype) msb;
    244 	  else
    245 	    mantissahi = (UDItype) msb << 32 | i[BITINT_END (1, 2)];
    246 	  mantissalo = ((UDItype) i[BITINT_END (in - 2, 1)] << 32
    247 			| i[BITINT_END (in - 1, 0)]);
    248 	}
    249 #endif
    250       if (iprec < 0)
    251 	mantissahi
    252 	  = ~mantissahi + __builtin_add_overflow (~mantissalo, 1, &mantissalo);
    253     }
    254 
    255   exponent += 6176;
    256   u.u[__FLOAT_WORD_ORDER__ != __ORDER_BIG_ENDIAN__]
    257     = ((((UDItype) (iprec < 0)) << 63)
    258        | (((UDItype) exponent) << 49)
    259        | mantissahi);
    260   u.u[__FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__] = mantissalo;
    261   if (inexact)
    262     {
    263       ui.u[__FLOAT_WORD_ORDER__ != __ORDER_BIG_ENDIAN__]
    264 	= (((UDItype) (iprec < 0)) << 63) | (((UDItype) exponent - 1) << 49);
    265       ui.u[__FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__] = inexact + 3;
    266       __asm ("" : "+g" (u.d));
    267       __asm ("" : "+g" (ui.d));
    268       u.d += ui.d;
    269       __asm ("" : "+g" (u.d));
    270     }
    271 
    272 done:
    273   return u.d;
    274 }
    275 #endif
    276