1e87fcf200
The ia32 _BitInt support revealed a bug in floatbitint?d.c. As can be even guessed from how the code is written in the loop, the intention was to set inexact to non-zero whenever the remainder after division wasn't zero, but I've ended up just checking whether the 2 least significant limbs of the remainder were non-zero. Now, in the dfp/bitint-4.c test in one case the remainder happens to have least significant 64 bits zero and then the higher limbs are non-zero; with 32-bit limbs that means 2 least significant limbs are zero and so the code acted as if it was exactly divisible. Fixed thusly. 2024-02-10 Jakub Jelinek <jakub@redhat.com> * soft-fp/floatbitintdd.c (__bid_floatbitintdd): Or in all remainder limbs into inexact rather than just first two. * soft-fp/floatbitintsd.c (__bid_floatbitintsd): Likewise. * soft-fp/floatbitinttd.c (__bid_floatbitinttd): Likewise.
275 lines
8.3 KiB
C
275 lines
8.3 KiB
C
/* Software floating-point emulation.
|
|
Convert a _BitInt to _Decimal128.
|
|
|
|
Copyright (C) 2023 Free Software Foundation, Inc.
|
|
|
|
This file is part of GCC.
|
|
|
|
GCC is free software; you can redistribute it and/or modify it under
|
|
the terms of the GNU General Public License as published by the Free
|
|
Software Foundation; either version 3, or (at your option) any later
|
|
version.
|
|
|
|
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
Under Section 7 of GPL version 3, you are granted additional
|
|
permissions described in the GCC Runtime Library Exception, version
|
|
3.1, as published by the Free Software Foundation.
|
|
|
|
You should have received a copy of the GNU General Public License and
|
|
a copy of the GCC Runtime Library Exception along with this program;
|
|
see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
<http://www.gnu.org/licenses/>. */
|
|
|
|
#include "soft-fp.h"
|
|
#include "bitint.h"
|
|
|
|
#ifdef __BITINT_MAXWIDTH__
|
|
extern _Decimal128 __bid_floatbitinttd (const UBILtype *, SItype);
|
|
|
|
_Decimal128
|
|
__bid_floatbitinttd (const UBILtype *i, SItype iprec)
|
|
{
|
|
iprec = bitint_reduce_prec (&i, iprec);
|
|
USItype aiprec = iprec < 0 ? -iprec : iprec;
|
|
USItype in = (aiprec + BIL_TYPE_SIZE - 1) / BIL_TYPE_SIZE;
|
|
USItype idx = BITINT_END (0, in - 1);
|
|
UBILtype msb = i[idx];
|
|
UDItype mantissahi, mantissalo;
|
|
SItype exponent = 0;
|
|
UBILtype inexact = 0;
|
|
union { _Decimal128 d; UDItype u[2]; } u, ui;
|
|
if (aiprec % BIL_TYPE_SIZE)
|
|
{
|
|
if (iprec > 0)
|
|
msb &= ((UBILtype) 1 << (aiprec % BIL_TYPE_SIZE)) - 1;
|
|
else
|
|
msb |= (UBILtype) -1 << (aiprec % BIL_TYPE_SIZE);
|
|
}
|
|
if (iprec < 0)
|
|
{
|
|
SItype n = sizeof (0ULL) * __CHAR_BIT__ + 1 - __builtin_clzll (~msb);
|
|
aiprec = (in - 1) * BIL_TYPE_SIZE + n;
|
|
}
|
|
else if (msb == 0)
|
|
aiprec = 1;
|
|
else
|
|
{
|
|
SItype n = sizeof (0ULL) * __CHAR_BIT__ - __builtin_clzll (msb);
|
|
aiprec = (in - 1) * BIL_TYPE_SIZE + n;
|
|
}
|
|
/* Number of bits in
|
|
(_BitInt(32768)) 9999999999999999999999999999999999e+6111DL. */
|
|
if (aiprec > 20414 + (iprec < 0))
|
|
{
|
|
ovf:
|
|
if (iprec < 0)
|
|
u.d = -9000000000000000000000000000000000e+6111DL;
|
|
else
|
|
u.d = 9000000000000000000000000000000000e+6111DL;
|
|
__asm ("" : "+g" (u.d));
|
|
u.d += u.d;
|
|
__asm ("" : "+g" (u.d));
|
|
goto done;
|
|
}
|
|
/* Bit precision of 9999999999999999999999999999999999uwb. */
|
|
if (aiprec >= 113)
|
|
{
|
|
USItype pow10_limbs, q_limbs, q2_limbs, j, k;
|
|
USItype exp_bits = 0, e;
|
|
UBILtype *buf;
|
|
/* First do a possibly large divide smaller enough such that
|
|
we only need to check remainder for 0 or non-0 and then
|
|
we'll do further division. */
|
|
if (aiprec >= 113 + 4 + 10)
|
|
{
|
|
exp_bits = ((aiprec - 113 - 4) * (UDItype) 30) / 299;
|
|
exponent = exp_bits * 3;
|
|
/* Upper estimate for pow10 (exponent) bits. */
|
|
exp_bits = exp_bits * 10 - exp_bits / 30;
|
|
}
|
|
pow10_limbs = (exp_bits + BIL_TYPE_SIZE - 1) / BIL_TYPE_SIZE;
|
|
/* 127 is the highest number of quotient bits needed on
|
|
aiprec range of [127, 20414]. E.g. if aiprec is 20409,
|
|
exponent will be 6105 and exp_bits 20283. 20409 - 20283 + 1
|
|
is 127. */
|
|
q_limbs = (127 + BIL_TYPE_SIZE - 1) / BIL_TYPE_SIZE;
|
|
q2_limbs = 128 / BIL_TYPE_SIZE;
|
|
buf = __builtin_alloca ((q_limbs + pow10_limbs * 2 + q2_limbs + 2)
|
|
* sizeof (UBILtype));
|
|
if (exponent)
|
|
{
|
|
__bid_pow10bitint (buf + q_limbs, exp_bits, exponent);
|
|
__divmodbitint4 (buf, q_limbs * BIL_TYPE_SIZE,
|
|
buf + q_limbs + pow10_limbs,
|
|
pow10_limbs * BIL_TYPE_SIZE,
|
|
i, iprec < 0 ? -aiprec : aiprec,
|
|
buf + q_limbs, exp_bits);
|
|
if (iprec < 0)
|
|
bitint_negate (buf + BITINT_END (q_limbs - 1, 0),
|
|
buf + BITINT_END (q_limbs - 1, 0), q_limbs);
|
|
inexact = buf[q_limbs + pow10_limbs];
|
|
for (j = 1; j < pow10_limbs; ++j)
|
|
inexact |= buf[q_limbs + pow10_limbs + j];
|
|
}
|
|
else
|
|
{
|
|
__builtin_memcpy (buf + BITINT_END (q_limbs - in + 1, 0), i,
|
|
(in - 1) * sizeof (UBILtype));
|
|
buf[BITINT_END (q_limbs - in, in - 1)] = msb;
|
|
if (iprec < 0)
|
|
bitint_negate (buf + BITINT_END (q_limbs - 1, 0),
|
|
buf + BITINT_END (q_limbs - 1, 0), in);
|
|
if (q_limbs > in)
|
|
__builtin_memset (buf + BITINT_END (0, in), '\0',
|
|
(q_limbs - in) * sizeof (UBILtype));
|
|
}
|
|
e = 0;
|
|
for (j = 3; j; )
|
|
{
|
|
USItype eprev = e;
|
|
__bid_pow10bitint (buf + q_limbs + pow10_limbs * 2 + 1,
|
|
128, 33 + e + j);
|
|
for (k = BITINT_END (0, q_limbs - 1);
|
|
k != BITINT_END (q_limbs, (USItype) -1); k -= BITINT_INC)
|
|
if (buf[k] > buf[q_limbs + pow10_limbs * 2 + 1 + k])
|
|
{
|
|
e += j;
|
|
break;
|
|
}
|
|
else if (buf[k] < buf[q_limbs + pow10_limbs * 2 + 1 + k])
|
|
break;
|
|
if (k == BITINT_END (q_limbs, (USItype) -1))
|
|
e += j;
|
|
if (j == 2 && e != eprev)
|
|
break;
|
|
else
|
|
--j;
|
|
}
|
|
exponent += e;
|
|
if (exponent > 6111)
|
|
goto ovf;
|
|
if (e)
|
|
{
|
|
UBILtype rem, half;
|
|
__bid_pow10bitint (buf + q_limbs + pow10_limbs * 2,
|
|
BIL_TYPE_SIZE, e);
|
|
__divmodbitint4 (buf + q_limbs + pow10_limbs * 2 + 1,
|
|
q2_limbs * BIL_TYPE_SIZE,
|
|
buf + q_limbs + pow10_limbs * 2 + 1 + q2_limbs,
|
|
BIL_TYPE_SIZE,
|
|
buf, q_limbs * BIL_TYPE_SIZE,
|
|
buf + q_limbs + pow10_limbs * 2, BIL_TYPE_SIZE);
|
|
half = buf[q_limbs + pow10_limbs * 2] / 2;
|
|
rem = buf[q_limbs + pow10_limbs * 2 + 1 + q2_limbs];
|
|
if (inexact)
|
|
{
|
|
/* If first division discovered some non-0 digits
|
|
and this second division is by 10, e.g.
|
|
for XXXXXX5499999999999 or XXXXXX5000000000001
|
|
if first division is by 10^12 and second by 10^1,
|
|
doing rem |= 1 wouldn't change the 5. Similarly
|
|
for rem 4 doing rem |= 1 would change it to 5,
|
|
but we don't want to change it in that case. */
|
|
if (e == 1)
|
|
{
|
|
if (rem == 5)
|
|
rem = 6;
|
|
else if (rem != 4)
|
|
rem |= 1;
|
|
}
|
|
else
|
|
rem |= 1;
|
|
}
|
|
/* Set inexact to 0, 1, 2, 3 depending on if remainder
|
|
of the divisions is exact 0, smaller than 10^exponent / 2,
|
|
exactly 10^exponent / 2 or greater than that. */
|
|
if (rem >= half)
|
|
inexact = 2 + (rem > half);
|
|
else
|
|
inexact = (rem != 0);
|
|
#if BIL_TYPE_SIZE == 64
|
|
mantissahi = buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (0, 1)];
|
|
mantissalo = buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (1, 0)];
|
|
#else
|
|
mantissahi
|
|
= ((UDItype)
|
|
buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (0, 3)] << 32
|
|
| buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (1, 2)]);
|
|
mantissalo
|
|
= ((UDItype)
|
|
buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (2, 1)] << 32
|
|
| buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (3, 0)]);
|
|
#endif
|
|
}
|
|
else
|
|
{
|
|
#if BIL_TYPE_SIZE == 64
|
|
mantissahi = buf[BITINT_END (0, 1)];
|
|
mantissalo = buf[BITINT_END (1, 0)];
|
|
#else
|
|
mantissahi = ((UDItype) buf[BITINT_END (0, 3)] << 32
|
|
| buf[BITINT_END (1, 2)]);
|
|
mantissalo = ((UDItype) buf[BITINT_END (2, 1)] << 32
|
|
| buf[BITINT_END (3, 0)]);
|
|
#endif
|
|
}
|
|
}
|
|
else
|
|
{
|
|
mantissahi = iprec < 0 ? -1 : 0;
|
|
#if BIL_TYPE_SIZE == 64
|
|
if (in == 1)
|
|
mantissalo = msb;
|
|
else
|
|
{
|
|
mantissahi = msb;
|
|
mantissalo = i[BITINT_END (1, 0)];
|
|
}
|
|
#else
|
|
if (in <= 2)
|
|
{
|
|
if (in == 1)
|
|
mantissalo = iprec < 0 ? (UDItype) (BILtype) msb : (UDItype) msb;
|
|
else
|
|
mantissalo = (UDItype) msb << 32 | i[BITINT_END (1, 0)];
|
|
}
|
|
else
|
|
{
|
|
if (in == 3)
|
|
mantissahi = iprec < 0 ? (UDItype) (BILtype) msb : (UDItype) msb;
|
|
else
|
|
mantissahi = (UDItype) msb << 32 | i[BITINT_END (1, 2)];
|
|
mantissalo = ((UDItype) i[BITINT_END (in - 2, 1)] << 32
|
|
| i[BITINT_END (in - 1, 0)]);
|
|
}
|
|
#endif
|
|
if (iprec < 0)
|
|
mantissahi
|
|
= ~mantissahi + __builtin_add_overflow (~mantissalo, 1, &mantissalo);
|
|
}
|
|
|
|
exponent += 6176;
|
|
u.u[__FLOAT_WORD_ORDER__ != __ORDER_BIG_ENDIAN__]
|
|
= ((((UDItype) (iprec < 0)) << 63)
|
|
| (((UDItype) exponent) << 49)
|
|
| mantissahi);
|
|
u.u[__FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__] = mantissalo;
|
|
if (inexact)
|
|
{
|
|
ui.u[__FLOAT_WORD_ORDER__ != __ORDER_BIG_ENDIAN__]
|
|
= (((UDItype) (iprec < 0)) << 63) | (((UDItype) exponent - 1) << 49);
|
|
ui.u[__FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__] = inexact + 3;
|
|
__asm ("" : "+g" (u.d));
|
|
__asm ("" : "+g" (ui.d));
|
|
u.d += ui.d;
|
|
__asm ("" : "+g" (u.d));
|
|
}
|
|
|
|
done:
|
|
return u.d;
|
|
}
|
|
#endif
|