libgcc/soft-fp/floatbitintdd.c

   1 /* Software floating-point emulation.
   2    Convert a _BitInt to _Decimal64.
   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 _Decimal64 __bid_floatbitintdd (const UBILtype *, SItype);
  32
  33 _Decimal64
  34 __bid_floatbitintdd (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 mantissa;
  42   SItype exponent = 0;
  43   UBILtype inexact = 0;
  44   union { _Decimal64 d; UDItype u; } 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 (_BitInt(2048)) 9999999999999999e+369DD.  */
  65   if (aiprec > 1279 + (iprec < 0))
  66     {
  67     ovf:
  68       if (iprec < 0)
  69         u.d = -9000000000000000e+369DD;
  70       else
  71         u.d = 9000000000000000e+369DD;
  72       __asm ("" : "+g" (u.d));
  73       u.d += u.d;
  74       __asm ("" : "+g" (u.d));
  75       goto done;
  76     }
  77   /* Bit precision of 9999999999999999uwb.  */
  78   if (aiprec >= 54)
  79     {
  80       USItype pow10_limbs, q_limbs, q2_limbs, j;
  81       USItype exp_bits = 0, e;
  82       UDItype m;
  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 >= 54 + 4 + 10)
  88         {
  89           exp_bits = (aiprec - 54 - 4) / 10;
  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       /* 72 is the highest number of quotient bits needed on
  96          aiprec range of [68, 1279].  E.g. if aiprec is 1277,
  97          exponent will be 363 and exp_bits 1206.  1277 - 1206 + 1
  98          is 72.  Unfortunately that means the result doesn't fit into
  99          UDItype...  */
 100       q_limbs = (72 + BIL_TYPE_SIZE - 1) / BIL_TYPE_SIZE;
 101       q2_limbs = 64 / BIL_TYPE_SIZE;
 102       buf = __builtin_alloca ((q_limbs + pow10_limbs * 2 + q2_limbs + 2)
 103                               * sizeof (UBILtype));
 104       if (exponent)
 105         {
 106           __bid_pow10bitint (buf + q_limbs, exp_bits, exponent);
 107           __divmodbitint4 (buf, q_limbs * BIL_TYPE_SIZE,
 108                            buf + q_limbs + pow10_limbs,
 109                            pow10_limbs * BIL_TYPE_SIZE,
 110                            i, iprec < 0 ? -aiprec : aiprec,
 111                            buf + q_limbs, exp_bits);
 112           if (iprec < 0)
 113             bitint_negate (buf + BITINT_END (q_limbs - 1, 0),
 114                            buf + BITINT_END (q_limbs - 1, 0), q_limbs);
 115           inexact = buf[q_limbs + pow10_limbs];
 116           for (j = 1; j < pow10_limbs; ++j)
 117             inexact |= buf[q_limbs + pow10_limbs + 1];
 118         }
 119       else
 120         {
 121           __builtin_memcpy (buf + BITINT_END (q_limbs - in + 1, 0), i,
 122                             (in - 1) * sizeof (UBILtype));
 123           buf[BITINT_END (q_limbs - in, in - 1)] = msb;
 124           if (iprec < 0)
 125             bitint_negate (buf + BITINT_END (q_limbs - 1, 0),
 126                            buf + BITINT_END (q_limbs - 1, 0), in);
 127           if (q_limbs > in)
 128             __builtin_memset (buf + BITINT_END (0, in), '\0',
 129                               (q_limbs - in) * sizeof (UBILtype));
 130         }
 131       e = 0;
 132 #if BIL_TYPE_SIZE == 64
 133       m = buf[BITINT_END (1, 0)];
 134 #elif BIL_TYPE_SIZE == 32
 135       m = ((UDItype) buf[1] << 32) | buf[BITINT_END (2, 0)];
 136 #else
 137 # error Unsupported BIL_TYPE_SIZE
 138 #endif
 139       if (buf[BITINT_END (0, q_limbs - 1)])
 140         {
 141           if (buf[BITINT_END (0, q_limbs - 1)] > 0x5)
 142             {
 143               /* 1000000000000000000000wb */
 144               if (buf[BITINT_END (0, q_limbs - 1)] > 0x36
 145                   || (buf[BITINT_END (0, q_limbs - 1)] == 0x36
 146                       && m >= (UDItype) 0x35c9adc5dea00000))
 147                 e = 6;
 148               else
 149                 e = 5;
 150             }
 151           /* 100000000000000000000wb */
 152           else if (buf[BITINT_END (0, q_limbs - 1)] == 0x5
 153                    && m >= (UDItype) 0x6bc75e2d63100000)
 154             e = 5;
 155           else
 156             e = 4;
 157         }
 158       else if (m >= (UDItype) 1000000000000000000)
 159         {
 160           if (m >= (UDItype) 10000000000000000000ULL)
 161             e = 4;
 162           else
 163             e = 3;
 164         }
 165       else if (m >= (UDItype) 100000000000000000)
 166         e = 2;
 167       else if (m >= (UDItype) 10000000000000000)
 168         e = 1;
 169       exponent += e;
 170       if (exponent > 369)
 171         goto ovf;
 172       if (e)
 173         {
 174           UBILtype rem, half;
 175           __bid_pow10bitint (buf + q_limbs + pow10_limbs * 2,
 176                              BIL_TYPE_SIZE, e);
 177           __divmodbitint4 (buf + q_limbs + pow10_limbs * 2 + 1,
 178                            q2_limbs * BIL_TYPE_SIZE,
 179                            buf + q_limbs + pow10_limbs * 2 + 1 + q2_limbs,
 180                            BIL_TYPE_SIZE,
 181                            buf, q_limbs * BIL_TYPE_SIZE,
 182                            buf + q_limbs + pow10_limbs * 2, BIL_TYPE_SIZE);
 183           half = buf[q_limbs + pow10_limbs * 2] / 2;
 184           rem = buf[q_limbs + pow10_limbs * 2 + 1 + q2_limbs];
 185           if (inexact)
 186             {
 187               /* If first division discovered some non-0 digits
 188                  and this second division is by 10, e.g.
 189                  for XXXXXX5499999999999 or XXXXXX5000000000001
 190                  if first division is by 10^12 and second by 10^1,
 191                  doing rem |= 1 wouldn't change the 5.  Similarly
 192                  for rem 4 doing rem |= 1 would change it to 5,
 193                  but we don't want to change it in that case.  */
 194               if (e == 1)
 195                 {
 196                   if (rem == 5)
 197                     rem = 6;
 198                   else if (rem != 4)
 199                     rem |= 1;
 200                 }
 201               else
 202                 rem |= 1;
 203             }
 204           /* Set inexact to 0, 1, 2, 3 depending on if remainder
 205              of the divisions is exact 0, smaller than 10^exponent / 2,
 206              exactly 10^exponent / 2 or greater than that.  */
 207           if (rem >= half)
 208             inexact = 2 + (rem > half);
 209           else
 210             inexact = (rem != 0);
 211 #if BIL_TYPE_SIZE == 64
 212           mantissa = buf[q_limbs + pow10_limbs * 2 + 1];
 213 #else
 214           mantissa
 215             = ((buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (0, 1)] << 32)
 216                | buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (1, 0)]);
 217 #endif
 218         }
 219       else
 220 #if BIL_TYPE_SIZE == 64
 221         mantissa = buf[BITINT_END (1, 0)];
 222 #else
 223         mantissa
 224           = ((buf[1] << 32) | buf[BITINT_END (2, 0)]);
 225 #endif
 226     }
 227   else
 228     {
 229 #if BIL_TYPE_SIZE == 64
 230       mantissa = msb;
 231 #else
 232       if (in == 1)
 233         mantissa = iprec < 0 ? (UDItype) (BILtype) msb : (UDItype) msb;
 234       else
 235         mantissa = ((msb << 32) | i[BITINT_END (1, 0)]);
 236 #endif
 237       if (iprec < 0)
 238         mantissa = -mantissa;
 239     }
 240
 241   exponent += 398;
 242   if (mantissa >= (UDItype) 0x20000000000000)
 243     u.u = (((((iprec < 0) << 2) | (UDItype) 3) << 61)
 244            | (((UDItype) exponent) << 51)
 245            | (mantissa ^ (UDItype) 0x20000000000000));
 246   else
 247     u.u = ((((UDItype) (iprec < 0)) << 63)
 248            | (((UDItype) exponent) << 53)
 249            | mantissa);
 250   if (inexact)
 251     {
 252       ui.u = ((((UDItype) (iprec < 0)) << 63)
 253               | (((UDItype) (exponent - 1)) << 53)
 254               | (inexact + 3));
 255       __asm ("" : "+g" (u.d));
 256       __asm ("" : "+g" (ui.d));
 257       u.d += ui.d;
 258       __asm ("" : "+g" (u.d));
 259     }
 260
 261 done:
 262   return u.d;
 263 }
 264 #endif