2ce182e258
This patch adds the library helpers for multiplication, division + modulo
and casts from and to floating point (both binary and decimal).
As described in the intro, the first step is try to reduce further the
passed in precision by skipping over most significant limbs with just zeros
or sign bit copies. For multiplication and division I've implemented
a simple algorithm, using something smarter like Karatsuba or Toom N-Way
might be faster for very large _BitInts (which we don't support right now
anyway), but could mean more code in libgcc, which maybe isn't what people
are willing to accept.
For the to/from floating point conversions the patch uses soft-fp, because
it already has tons of handy macros which can be used for that. In theory
it could be implemented using {,unsigned} long long or {,unsigned} __int128
to/from floating point conversions with some frexp before/after, but at that
point we already need to force it into integer registers and analyze it
anyway. Plus, for 32-bit arches there is no __int128 that could be used
for XF/TF mode stuff.
I know that soft-fp is owned by glibc and I think the op-common.h change
should be propagated there, but the bitint stuff is really GCC specific
and IMHO doesn't belong into the glibc copy.
2023-09-06 Jakub Jelinek <jakub@redhat.com>
PR c/102989
libgcc/
* config/aarch64/t-softfp (softfp_extras): Use += rather than :=.
* config/i386/64/t-softfp (softfp_extras): Likewise.
* config/i386/libgcc-glibc.ver (GCC_14.0.0): Export _BitInt support
routines.
* config/i386/t-softfp (softfp_extras): Add fixxfbitint and
bf, hf and xf mode floatbitint.
(CFLAGS-floatbitintbf.c, CFLAGS-floatbitinthf.c): Add -msse2.
* config/riscv/t-softfp32 (softfp_extras): Use += rather than :=.
* config/rs6000/t-e500v1-fp (softfp_extras): Likewise.
* config/rs6000/t-e500v2-fp (softfp_extras): Likewise.
* config/t-softfp (softfp_floatbitint_funcs): New.
(softfp_bid_list): New.
(softfp_func_list): Add sf and df mode from and to _BitInt libcalls.
(softfp_bid_file_list): New.
(LIB2ADD_ST): Add $(softfp_bid_file_list).
* config/t-softfp-sfdftf (softfp_extras): Add fixtfbitint and
floatbitinttf.
* config/t-softfp-tf (softfp_extras): Likewise.
* libgcc2.c (bitint_reduce_prec): New inline function.
(BITINT_INC, BITINT_END): Define.
(bitint_mul_1, bitint_addmul_1): New helper functions.
(__mulbitint3): New function.
(bitint_negate, bitint_submul_1): New helper functions.
(__divmodbitint4): New function.
* libgcc2.h (LIBGCC2_UNITS_PER_WORD): When building _BitInt support
libcalls, redefine depending on __LIBGCC_BITINT_LIMB_WIDTH__.
(__mulbitint3, __divmodbitint4): Declare.
* libgcc-std.ver.in (GCC_14.0.0): Export _BitInt support routines.
* Makefile.in (lib2funcs): Add _mulbitint3.
(LIB2_DIVMOD_FUNCS): Add _divmodbitint4.
* soft-fp/bitint.h: New file.
* soft-fp/fixdfbitint.c: New file.
* soft-fp/fixsfbitint.c: New file.
* soft-fp/fixtfbitint.c: New file.
* soft-fp/fixxfbitint.c: New file.
* soft-fp/floatbitintbf.c: New file.
* soft-fp/floatbitintdf.c: New file.
* soft-fp/floatbitinthf.c: New file.
* soft-fp/floatbitintsf.c: New file.
* soft-fp/floatbitinttf.c: New file.
* soft-fp/floatbitintxf.c: New file.
* soft-fp/op-common.h (_FP_FROM_INT): Add support for rsize up to
4 * _FP_W_TYPE_SIZE rather than just 2 * _FP_W_TYPE_SIZE.
* soft-fp/bitintpow10.c: New file.
* soft-fp/fixsdbitint.c: New file.
* soft-fp/fixddbitint.c: New file.
* soft-fp/fixtdbitint.c: New file.
* soft-fp/floatbitintsd.c: New file.
* soft-fp/floatbitintdd.c: New file.
* soft-fp/floatbitinttd.c: New file.
271 lines
8.2 KiB
C
271 lines
8.2 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 + 1];
|
|
}
|
|
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
|
|
= ((buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (0, 3)] << 32)
|
|
| buf[q_limbs + pow10_limbs * 2 + 1 + BITINT_END (1, 2)]);
|
|
mantissalo
|
|
= ((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 = (buf[BITINT_END (0, 3)] << 32) | buf[BITINT_END (1, 2)];
|
|
mantissalo = (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 = (msb << 32) | i[BITINT_END (1, 0)];
|
|
}
|
|
else
|
|
{
|
|
if (in == 3)
|
|
mantissahi = iprec < 0 ? (UDItype) (BILtype) msb : (UDItype) msb;
|
|
else
|
|
mantissahi = (msb << 32) | i[BITINT_END (1, 2)];
|
|
mantissalo = ((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
|