libgcc/soft-fp/floatbitinttd.c

   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 + 1];
 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             = ((buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (0, 3)] << 32)
 200                | buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (1, 2)]);
 201           mantissalo
 202             = ((buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (2, 1)] << 32)
 203                | buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (3, 0)]);
 204 #endif
 205         }
 206       else
 207         {
 208 #if BIL_TYPE_SIZE == 64
 209           mantissahi = buf[BITINT_END (0, 1)];
 210           mantissalo = buf[BITINT_END (1, 0)];
 211 #else
 212           mantissahi = (buf[BITINT_END (0, 3)] << 32) | buf[BITINT_END (1, 2)];
 213           mantissalo = (buf[BITINT_END (2, 1)] << 32) | buf[BITINT_END (3, 0)];
 214 #endif
 215         }
 216     }
 217   else
 218     {
 219       mantissahi = iprec < 0 ? -1 : 0;
 220 #if BIL_TYPE_SIZE == 64
 221       if (in == 1)
 222         mantissalo = msb;
 223       else
 224         {
 225           mantissahi = msb;
 226           mantissalo = i[BITINT_END (1, 0)];
 227         }
 228 #else
 229       if (in <= 2)
 230         {
 231           if (in == 1)
 232             mantissalo = iprec < 0 ? (UDItype) (BILtype) msb : (UDItype) msb;
 233           else
 234             mantissalo = (msb << 32) | i[BITINT_END (1, 0)];
 235         }
 236       else
 237         {
 238           if (in == 3)
 239             mantissahi = iprec < 0 ? (UDItype) (BILtype) msb : (UDItype) msb;
 240           else
 241             mantissahi = (msb << 32) | i[BITINT_END (1, 2)];
 242           mantissalo = ((i[BITINT_END (in - 2, 1)] << 32)
 243                         | i[BITINT_END (in - 1, 0)]);
 244         }
 245 #endif
 246       if (iprec < 0)
 247         mantissahi
 248           = ~mantissahi + __builtin_add_overflow (~mantissalo, 1, &mantissalo);
 249     }
 250
 251   exponent += 6176;
 252   u.u[__FLOAT_WORD_ORDER__ != __ORDER_BIG_ENDIAN__]
 253     = ((((UDItype) (iprec < 0)) << 63)
 254        | (((UDItype) exponent) << 49)
 255        | mantissahi);
 256   u.u[__FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__] = mantissalo;
 257   if (inexact)
 258     {
 259       ui.u[__FLOAT_WORD_ORDER__ != __ORDER_BIG_ENDIAN__]
 260         = (((UDItype) (iprec < 0)) << 63) | (((UDItype) exponent - 1) << 49);
 261       ui.u[__FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__] = inexact + 3;
 262       __asm ("" : "+g" (u.d));
 263       __asm ("" : "+g" (ui.d));
 264       u.d += ui.d;
 265       __asm ("" : "+g" (u.d));
 266     }
 267
 268 done:
 269   return u.d;
 270 }
 271 #endif