diff options
Diffstat (limited to 'libc')
35 files changed, 1507 insertions, 801 deletions
diff --git a/libc/config/baremetal/api.td b/libc/config/baremetal/api.td index 690edbd..25aa06a 100644 --- a/libc/config/baremetal/api.td +++ b/libc/config/baremetal/api.td @@ -57,10 +57,7 @@ def MathAPI : PublicAPI<"math.h"> { } def StdIOAPI : PublicAPI<"stdio.h"> { - let Types = [ - "size_t", - "off_t", - ]; + let Types = ["size_t"]; } def StdlibAPI : PublicAPI<"stdlib.h"> { diff --git a/libc/config/gpu/api.td b/libc/config/gpu/api.td index 523ad49..adaf5bf 100644 --- a/libc/config/gpu/api.td +++ b/libc/config/gpu/api.td @@ -64,11 +64,7 @@ def StdIOAPI : PublicAPI<"stdio.h"> { SimpleMacroDef<"_IOLBF", "1">, SimpleMacroDef<"_IONBF", "2">, ]; - let Types = [ - "FILE", - "off_t", - "size_t", - ]; + let Types = ["size_t", "FILE"]; } def IntTypesAPI : PublicAPI<"inttypes.h"> { diff --git a/libc/config/linux/api.td b/libc/config/linux/api.td index 9964971..eb5ed80 100644 --- a/libc/config/linux/api.td +++ b/libc/config/linux/api.td @@ -49,10 +49,7 @@ def CTypeAPI : PublicAPI<"ctype.h"> { } def FCntlAPI : PublicAPI<"fcntl.h"> { - let Types = [ - "mode_t", - "off_t", - ]; + let Types = ["mode_t"]; } def IntTypesAPI : PublicAPI<"inttypes.h"> { @@ -80,12 +77,7 @@ def StdIOAPI : PublicAPI<"stdio.h"> { SimpleMacroDef<"_IOLBF", "1">, SimpleMacroDef<"_IONBF", "2">, ]; - let Types = [ - "FILE", - "cookie_io_functions_t", - "off_t", - "size_t", - ]; + let Types = ["size_t", "FILE", "cookie_io_functions_t"]; } def StdlibAPI : PublicAPI<"stdlib.h"> { diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index cc7671c..2742c33 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -370,6 +370,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.exp10f libc.src.math.exp2 libc.src.math.exp2f + libc.src.math.exp2m1f libc.src.math.expm1 libc.src.math.expm1f libc.src.math.fabs diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst index 265261b..970a43c 100644 --- a/libc/docs/math/index.rst +++ b/libc/docs/math/index.rst @@ -270,7 +270,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | exp2 | |check| | |check| | | | | 7.12.6.4 | F.10.3.4 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| exp2m1 | | | | | | 7.12.6.5 | F.10.3.5 | +| exp2m1 | |check| | | | | | 7.12.6.5 | F.10.3.5 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | expm1 | |check| | |check| | | | | 7.12.6.6 | F.10.3.6 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/fuzzing/CMakeLists.txt b/libc/fuzzing/CMakeLists.txt index 8248768..816691b 100644 --- a/libc/fuzzing/CMakeLists.txt +++ b/libc/fuzzing/CMakeLists.txt @@ -1,6 +1,7 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fsanitize=fuzzer") add_custom_target(libc-fuzzer) +add_subdirectory(__support) # TODO(#85680): Re-enable math fuzzing after headers are sorted out # add_subdirectory(math) add_subdirectory(stdlib) diff --git a/libc/fuzzing/__support/CMakeLists.txt b/libc/fuzzing/__support/CMakeLists.txt new file mode 100644 index 0000000..278e914 --- /dev/null +++ b/libc/fuzzing/__support/CMakeLists.txt @@ -0,0 +1,7 @@ +add_libc_fuzzer( + uint_fuzz + SRCS + uint_fuzz.cpp + DEPENDS + libc.src.__support.uint +) diff --git a/libc/fuzzing/__support/uint_fuzz.cpp b/libc/fuzzing/__support/uint_fuzz.cpp new file mode 100644 index 0000000..f48f00d --- /dev/null +++ b/libc/fuzzing/__support/uint_fuzz.cpp @@ -0,0 +1,70 @@ +#include "src/__support/CPP/bit.h" +#include "src/__support/UInt.h" +#include "src/string/memory_utils/inline_memcpy.h" + +using namespace LIBC_NAMESPACE; + +// Helper function when using gdb / lldb to set a breakpoint and inspect values. +template <typename T> void debug_and_trap(const char *msg, T a, T b) { + __builtin_trap(); +} + +#define DEBUG_AND_TRAP() + +#define TEST_BINOP(OP) \ + if ((a OP b) != (static_cast<T>(BigInt(a) OP BigInt(b)))) \ + debug_and_trap(#OP, a, b); + +#define TEST_SHIFTOP(OP) \ + if ((a OP b) != (static_cast<T>(BigInt(a) OP b))) \ + debug_and_trap(#OP, a, b); + +#define TEST_FUNCTION(FUN) \ + if (FUN(a) != FUN(BigInt(a))) \ + debug_and_trap(#FUN, a, b); + +// Test that basic arithmetic operations of BigInt behave like their scalar +// counterparts. +template <typename T, typename BigInt> void run_tests(T a, T b) { + TEST_BINOP(+) + TEST_BINOP(-) + TEST_BINOP(*) + if (b != 0) + TEST_BINOP(/) + if (b >= 0 && b < cpp::numeric_limits<T>::digits) { + TEST_SHIFTOP(<<) + TEST_SHIFTOP(>>) + } + if constexpr (!BigInt::SIGNED) { + TEST_FUNCTION(cpp::has_single_bit) + TEST_FUNCTION(cpp::countr_zero) + TEST_FUNCTION(cpp::countl_zero) + TEST_FUNCTION(cpp::countl_one) + TEST_FUNCTION(cpp::countr_one) + } +} + +// Reads a T from libfuzzer data. +template <typename T> T read(const uint8_t *data, size_t &remainder) { + T out = 0; + constexpr size_t T_SIZE = sizeof(T); + const size_t copy_size = remainder < T_SIZE ? remainder : T_SIZE; + inline_memcpy(&out, data, copy_size); + remainder -= copy_size; + return out; +} + +template <typename T, typename BigInt> +void run_tests(const uint8_t *data, size_t size) { + const auto a = read<T>(data, size); + const auto b = read<T>(data, size); + run_tests<T, BigInt>(a, b); +} + +extern "C" int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) { + // unsigned + run_tests<uint64_t, BigInt<64, false, uint16_t>>(data, size); + // signed + run_tests<int64_t, BigInt<64, true, uint16_t>>(data, size); + return 0; +} diff --git a/libc/include/CMakeLists.txt b/libc/include/CMakeLists.txt index 02c7dc8..4203f0b 100644 --- a/libc/include/CMakeLists.txt +++ b/libc/include/CMakeLists.txt @@ -41,10 +41,9 @@ add_gen_header( DEF_FILE fcntl.h.def GEN_HDR fcntl.h DEPENDS + .llvm_libc_common_h .llvm-libc-macros.fcntl_macros .llvm-libc-types.mode_t - .llvm-libc-types.off_t - .llvm_libc_common_h ) add_gen_header( @@ -265,14 +264,13 @@ add_gen_header( DEF_FILE stdio.h.def GEN_HDR stdio.h DEPENDS + .llvm_libc_common_h .llvm-libc-macros.file_seek_macros .llvm-libc-macros.stdio_macros - .llvm-libc-types.FILE - .llvm-libc-types.cookie_io_functions_t - .llvm-libc-types.off_t .llvm-libc-types.size_t .llvm-libc-types.ssize_t - .llvm_libc_common_h + .llvm-libc-types.FILE + .llvm-libc-types.cookie_io_functions_t ) add_gen_header( diff --git a/libc/spec/posix.td b/libc/spec/posix.td index 45f7ecf..cfa8d3a 100644 --- a/libc/spec/posix.td +++ b/libc/spec/posix.td @@ -210,10 +210,7 @@ def POSIX : StandardSpec<"POSIX"> { HeaderSpec FCntl = HeaderSpec< "fcntl.h", [], // Macros - [ - ModeTType, - OffTType, - ], + [ModeTType], [], // Enumerations [ FunctionSpec< @@ -1183,7 +1180,7 @@ def POSIX : StandardSpec<"POSIX"> { HeaderSpec StdIO = HeaderSpec< "stdio.h", [], // Macros - [OffTType], // Types + [], // Types [], // Enumerations [ FunctionSpec< diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td index 719bb9a..bd62870 100644 --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -535,6 +535,8 @@ def StdC : StandardSpec<"stdc"> { FunctionSpec<"exp2", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>, FunctionSpec<"exp2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>, + FunctionSpec<"exp2m1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>, + FunctionSpec<"expm1", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>, FunctionSpec<"expm1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>, diff --git a/libc/src/__support/FPUtil/dyadic_float.h b/libc/src/__support/FPUtil/dyadic_float.h index 73fd738..e0c205f 100644 --- a/libc/src/__support/FPUtil/dyadic_float.h +++ b/libc/src/__support/FPUtil/dyadic_float.h @@ -58,9 +58,9 @@ template <size_t Bits> struct DyadicFloat { // significant bit. LIBC_INLINE constexpr DyadicFloat &normalize() { if (!mantissa.is_zero()) { - int shift_length = static_cast<int>(mantissa.clz()); + int shift_length = cpp::countl_zero(mantissa); exponent -= shift_length; - mantissa.shift_left(static_cast<size_t>(shift_length)); + mantissa <<= static_cast<size_t>(shift_length); } return *this; } @@ -233,7 +233,7 @@ LIBC_INLINE constexpr DyadicFloat<Bits> quick_add(DyadicFloat<Bits> a, result.sign = a.sign; result.exponent = a.exponent; result.mantissa = a.mantissa; - if (result.mantissa.add(b.mantissa)) { + if (result.mantissa.add_overflow(b.mantissa)) { // Mantissa addition overflow. result.shift_right(1); result.mantissa.val[DyadicFloat<Bits>::MantissaType::WORD_COUNT - 1] |= diff --git a/libc/src/__support/RPC/rpc.h b/libc/src/__support/RPC/rpc.h index 5dcae51..05506c0 100644 --- a/libc/src/__support/RPC/rpc.h +++ b/libc/src/__support/RPC/rpc.h @@ -198,12 +198,9 @@ template <bool Invert> struct Process { /// convergent, otherwise the compiler will sink the store and deadlock. [[clang::convergent]] LIBC_INLINE void unlock(uint64_t lane_mask, uint32_t index) { - // Do not move any writes past the unlock + // Do not move any writes past the unlock. atomic_thread_fence(cpp::MemoryOrder::RELEASE); - // Wait for other threads in the warp to finish using the lock - gpu::sync_lane(lane_mask); - // Use exactly one thread to clear the nth bit in the lock array Must // restrict to a single thread to avoid one thread dropping the lock, then // an unrelated warp claiming the lock, then a second thread in this warp @@ -331,6 +328,9 @@ public: LIBC_INLINE uint16_t get_index() const { return index; } LIBC_INLINE void close() { + // Wait for all lanes to finish using the port. + gpu::sync_lane(lane_mask); + // The server is passive, if it own the buffer when it closes we need to // give ownership back to the client. if (owns_buffer && T) diff --git a/libc/src/__support/UInt.h b/libc/src/__support/UInt.h index 282efdb..c1e55ce 100644 --- a/libc/src/__support/UInt.h +++ b/libc/src/__support/UInt.h @@ -14,10 +14,11 @@ #include "src/__support/CPP/limits.h" #include "src/__support/CPP/optional.h" #include "src/__support/CPP/type_traits.h" -#include "src/__support/macros/attributes.h" // LIBC_INLINE -#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY +#include "src/__support/macros/attributes.h" // LIBC_INLINE +#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY +#include "src/__support/macros/properties/compiler.h" // LIBC_COMPILER_IS_CLANG #include "src/__support/macros/properties/types.h" // LIBC_TYPES_HAS_INT128, LIBC_TYPES_HAS_INT64 -#include "src/__support/math_extras.h" // SumCarry, DiffBorrow +#include "src/__support/math_extras.h" // add_with_carry, sub_with_borrow #include "src/__support/number_pair.h" #include <stddef.h> // For size_t @@ -25,71 +26,324 @@ namespace LIBC_NAMESPACE { -namespace internal { -template <typename T> struct half_width; +namespace multiword { -template <> struct half_width<uint64_t> : cpp::type_identity<uint32_t> {}; -template <> struct half_width<uint32_t> : cpp::type_identity<uint16_t> {}; +// A type trait mapping unsigned integers to their half-width unsigned +// counterparts. +template <typename T> struct half_width; template <> struct half_width<uint16_t> : cpp::type_identity<uint8_t> {}; +template <> struct half_width<uint32_t> : cpp::type_identity<uint16_t> {}; +#ifdef LIBC_TYPES_HAS_INT64 +template <> struct half_width<uint64_t> : cpp::type_identity<uint32_t> {}; #ifdef LIBC_TYPES_HAS_INT128 template <> struct half_width<__uint128_t> : cpp::type_identity<uint64_t> {}; #endif // LIBC_TYPES_HAS_INT128 - +#endif // LIBC_TYPES_HAS_INT64 template <typename T> using half_width_t = typename half_width<T>::type; -template <typename T> constexpr NumberPair<T> full_mul(T a, T b) { - NumberPair<T> pa = split(a); - NumberPair<T> pb = split(b); - NumberPair<T> prod; +// An array of two elements that can be used in multiword operations. +template <typename T> struct DoubleWide final : cpp::array<T, 2> { + using UP = cpp::array<T, 2>; + using UP::UP; + LIBC_INLINE constexpr DoubleWide(T lo, T hi) : UP({lo, hi}) {} +}; + +// Converts an unsigned value into a DoubleWide<half_width_t<T>>. +template <typename T> LIBC_INLINE constexpr auto split(T value) { + static_assert(cpp::is_unsigned_v<T>); + using half_type = half_width_t<T>; + return DoubleWide<half_type>( + half_type(value), + half_type(value >> cpp::numeric_limits<half_type>::digits)); +} + +// The low part of a DoubleWide value. +template <typename T> LIBC_INLINE constexpr T lo(const DoubleWide<T> &value) { + return value[0]; +} +// The high part of a DoubleWide value. +template <typename T> LIBC_INLINE constexpr T hi(const DoubleWide<T> &value) { + return value[1]; +} +// The low part of an unsigned value. +template <typename T> LIBC_INLINE constexpr half_width_t<T> lo(T value) { + return lo(split(value)); +} +// The high part of an unsigned value. +template <typename T> LIBC_INLINE constexpr half_width_t<T> hi(T value) { + return hi(split(value)); +} + +// Returns 'a' times 'b' in a DoubleWide<word>. Cannot overflow by construction. +template <typename word> +LIBC_INLINE constexpr DoubleWide<word> mul2(word a, word b) { + if constexpr (cpp::is_same_v<word, uint8_t>) { + return split<uint16_t>(uint16_t(a) * uint16_t(b)); + } else if constexpr (cpp::is_same_v<word, uint16_t>) { + return split<uint32_t>(uint32_t(a) * uint32_t(b)); + } +#ifdef LIBC_TYPES_HAS_INT64 + else if constexpr (cpp::is_same_v<word, uint32_t>) { + return split<uint64_t>(uint64_t(a) * uint64_t(b)); + } +#endif +#ifdef LIBC_TYPES_HAS_INT128 + else if constexpr (cpp::is_same_v<word, uint64_t>) { + return split<__uint128_t>(__uint128_t(a) * __uint128_t(b)); + } +#endif + else { + using half_word = half_width_t<word>; + const auto shiftl = [](word value) -> word { + return value << cpp::numeric_limits<half_word>::digits; + }; + const auto shiftr = [](word value) -> word { + return value >> cpp::numeric_limits<half_word>::digits; + }; + // Here we do a one digit multiplication where 'a' and 'b' are of type + // word. We split 'a' and 'b' into half words and perform the classic long + // multiplication with 'a' and 'b' being two-digit numbers. + + // a a_hi a_lo + // x b => x b_hi b_lo + // ---- ----------- + // c result + // We convert 'lo' and 'hi' from 'half_word' to 'word' so multiplication + // doesn't overflow. + const word a_lo = lo(a); + const word b_lo = lo(b); + const word a_hi = hi(a); + const word b_hi = hi(b); + const word step1 = b_lo * a_lo; // no overflow; + const word step2 = b_lo * a_hi; // no overflow; + const word step3 = b_hi * a_lo; // no overflow; + const word step4 = b_hi * a_hi; // no overflow; + word lo_digit = step1; + word hi_digit = step4; + const word no_carry = 0; + word carry; + word _; // unused carry variable. + lo_digit = add_with_carry<word>(lo_digit, shiftl(step2), no_carry, carry); + hi_digit = add_with_carry<word>(hi_digit, shiftr(step2), carry, _); + lo_digit = add_with_carry<word>(lo_digit, shiftl(step3), no_carry, carry); + hi_digit = add_with_carry<word>(hi_digit, shiftr(step3), carry, _); + return DoubleWide<word>(lo_digit, hi_digit); + } +} + +// In-place 'dst op= rhs' with operation with carry propagation. Returns carry. +template <typename Function, typename word, size_t N, size_t M> +LIBC_INLINE constexpr word inplace_binop(Function op_with_carry, + cpp::array<word, N> &dst, + const cpp::array<word, M> &rhs) { + static_assert(N >= M); + word carry_out = 0; + for (size_t i = 0; i < N; ++i) { + const bool has_rhs_value = i < M; + const word rhs_value = has_rhs_value ? rhs[i] : 0; + const word carry_in = carry_out; + dst[i] = op_with_carry(dst[i], rhs_value, carry_in, carry_out); + // stop early when rhs is over and no carry is to be propagated. + if (!has_rhs_value && carry_out == 0) + break; + } + return carry_out; +} - prod.lo = pa.lo * pb.lo; // exact - prod.hi = pa.hi * pb.hi; // exact - NumberPair<T> lo_hi = split(pa.lo * pb.hi); // exact - NumberPair<T> hi_lo = split(pa.hi * pb.lo); // exact +// In-place addition. Returns carry. +template <typename word, size_t N, size_t M> +LIBC_INLINE constexpr word add_with_carry(cpp::array<word, N> &dst, + const cpp::array<word, M> &rhs) { + return inplace_binop(LIBC_NAMESPACE::add_with_carry<word>, dst, rhs); +} + +// In-place subtraction. Returns borrow. +template <typename word, size_t N, size_t M> +LIBC_INLINE constexpr word sub_with_borrow(cpp::array<word, N> &dst, + const cpp::array<word, M> &rhs) { + return inplace_binop(LIBC_NAMESPACE::sub_with_borrow<word>, dst, rhs); +} + +// In-place multiply-add. Returns carry. +// i.e., 'dst += b * c' +template <typename word, size_t N> +LIBC_INLINE constexpr word mul_add_with_carry(cpp::array<word, N> &dst, word b, + word c) { + return add_with_carry(dst, mul2(b, c)); +} - constexpr size_t HALF_BIT_WIDTH = sizeof(T) * CHAR_BIT / 2; +// An array of two elements serving as an accumulator during multiword +// computations. +template <typename T> struct Accumulator final : cpp::array<T, 2> { + using UP = cpp::array<T, 2>; + LIBC_INLINE constexpr Accumulator() : UP({0, 0}) {} + LIBC_INLINE constexpr T advance(T carry_in) { + auto result = UP::front(); + UP::front() = UP::back(); + UP::back() = carry_in; + return result; + } + LIBC_INLINE constexpr T sum() const { return UP::front(); } + LIBC_INLINE constexpr T carry() const { return UP::back(); } +}; - auto r1 = add_with_carry(prod.lo, lo_hi.lo << HALF_BIT_WIDTH, T(0)); - prod.lo = r1.sum; - prod.hi = add_with_carry(prod.hi, lo_hi.hi, r1.carry).sum; +// In-place multiplication by a single word. Returns carry. +template <typename word, size_t N> +LIBC_INLINE constexpr word scalar_multiply_with_carry(cpp::array<word, N> &dst, + word x) { + Accumulator<word> acc; + for (auto &val : dst) { + const word carry = mul_add_with_carry(acc, val, x); + val = acc.advance(carry); + } + return acc.carry(); +} - auto r2 = add_with_carry(prod.lo, hi_lo.lo << HALF_BIT_WIDTH, T(0)); - prod.lo = r2.sum; - prod.hi = add_with_carry(prod.hi, hi_lo.hi, r2.carry).sum; +// Multiplication of 'lhs' by 'rhs' into 'dst'. Returns carry. +// This function is safe to use for signed numbers. +// https://stackoverflow.com/a/20793834 +// https://pages.cs.wisc.edu/%7Emarkhill/cs354/Fall2008/beyond354/int.mult.html +template <typename word, size_t O, size_t M, size_t N> +LIBC_INLINE constexpr word multiply_with_carry(cpp::array<word, O> &dst, + const cpp::array<word, M> &lhs, + const cpp::array<word, N> &rhs) { + static_assert(O >= M + N); + Accumulator<word> acc; + for (size_t i = 0; i < O; ++i) { + const size_t lower_idx = i < N ? 0 : i - N + 1; + const size_t upper_idx = i < M ? i : M - 1; + word carry = 0; + for (size_t j = lower_idx; j <= upper_idx; ++j) + carry += mul_add_with_carry(acc, lhs[j], rhs[i - j]); + dst[i] = acc.advance(carry); + } + return acc.carry(); +} - return prod; +template <typename word, size_t N> +LIBC_INLINE constexpr void quick_mul_hi(cpp::array<word, N> &dst, + const cpp::array<word, N> &lhs, + const cpp::array<word, N> &rhs) { + Accumulator<word> acc; + word carry = 0; + // First round of accumulation for those at N - 1 in the full product. + for (size_t i = 0; i < N; ++i) + carry += mul_add_with_carry(acc, lhs[i], rhs[N - 1 - i]); + for (size_t i = N; i < 2 * N - 1; ++i) { + acc.advance(carry); + carry = 0; + for (size_t j = i - N + 1; j < N; ++j) + carry += mul_add_with_carry(acc, lhs[j], rhs[i - j]); + dst[i - N] = acc.sum(); + } + dst.back() = acc.carry(); } -template <> -LIBC_INLINE constexpr NumberPair<uint32_t> full_mul<uint32_t>(uint32_t a, - uint32_t b) { - uint64_t prod = uint64_t(a) * uint64_t(b); - NumberPair<uint32_t> result; - result.lo = uint32_t(prod); - result.hi = uint32_t(prod >> 32); - return result; +template <typename word, size_t N> +LIBC_INLINE constexpr bool is_negative(cpp::array<word, N> &array) { + using signed_word = cpp::make_signed_t<word>; + return cpp::bit_cast<signed_word>(array.back()) < 0; } +// An enum for the shift function below. +enum Direction { LEFT, RIGHT }; + +// A bitwise shift on an array of elements. +// TODO: Make the result UB when 'offset' is greater or equal to the number of +// bits in 'array'. This will allow for better code performance. +template <Direction direction, bool is_signed, typename word, size_t N> +LIBC_INLINE constexpr cpp::array<word, N> shift(cpp::array<word, N> array, + size_t offset) { + static_assert(direction == LEFT || direction == RIGHT); + constexpr size_t WORD_BITS = cpp::numeric_limits<word>::digits; + constexpr size_t TOTAL_BITS = N * WORD_BITS; + if (LIBC_UNLIKELY(offset == 0)) + return array; + if (LIBC_UNLIKELY(offset >= TOTAL_BITS)) + return {}; #ifdef LIBC_TYPES_HAS_INT128 -template <> -LIBC_INLINE constexpr NumberPair<uint64_t> full_mul<uint64_t>(uint64_t a, - uint64_t b) { - __uint128_t prod = __uint128_t(a) * __uint128_t(b); - NumberPair<uint64_t> result; - result.lo = uint64_t(prod); - result.hi = uint64_t(prod >> 64); - return result; + if constexpr (TOTAL_BITS == 128) { + using type = cpp::conditional_t<is_signed, __int128_t, __uint128_t>; + auto tmp = cpp::bit_cast<type>(array); + if constexpr (direction == LEFT) + tmp <<= offset; + else + tmp >>= offset; + return cpp::bit_cast<cpp::array<word, N>>(tmp); + } +#endif + const bool is_neg = is_signed && is_negative(array); + constexpr auto at = [](size_t index) -> int { + // reverse iteration when direction == LEFT. + if constexpr (direction == LEFT) + return int(N) - int(index) - 1; + return int(index); + }; + const auto safe_get_at = [&](size_t index) -> word { + // return appropriate value when accessing out of bound elements. + const int i = at(index); + if (i < 0) + return 0; + if (i >= int(N)) + return is_neg ? -1 : 0; + return array[i]; + }; + const size_t index_offset = offset / WORD_BITS; + const size_t bit_offset = offset % WORD_BITS; +#ifdef LIBC_COMPILER_IS_CLANG + __builtin_assume(index_offset < N); +#endif + cpp::array<word, N> out = {}; + for (size_t index = 0; index < N; ++index) { + const word part1 = safe_get_at(index + index_offset); + const word part2 = safe_get_at(index + index_offset + 1); + word &dst = out[at(index)]; + if (bit_offset == 0) + dst = part1; // no crosstalk between parts. + else if constexpr (direction == LEFT) + dst = (part1 << bit_offset) | (part2 >> (WORD_BITS - bit_offset)); + else + dst = (part1 >> bit_offset) | (part2 << (WORD_BITS - bit_offset)); + } + return out; } -#endif // LIBC_TYPES_HAS_INT128 -} // namespace internal +#define DECLARE_COUNTBIT(NAME, INDEX_EXPR) \ + template <typename word, size_t N> \ + LIBC_INLINE constexpr int NAME(const cpp::array<word, N> &val) { \ + int bit_count = 0; \ + for (size_t i = 0; i < N; ++i) { \ + const int word_count = cpp::NAME<word>(val[INDEX_EXPR]); \ + bit_count += word_count; \ + if (word_count != cpp::numeric_limits<word>::digits) \ + break; \ + } \ + return bit_count; \ + } + +DECLARE_COUNTBIT(countr_zero, i) // iterating forward +DECLARE_COUNTBIT(countr_one, i) // iterating forward +DECLARE_COUNTBIT(countl_zero, N - i - 1) // iterating backward +DECLARE_COUNTBIT(countl_one, N - i - 1) // iterating backward + +} // namespace multiword template <size_t Bits, bool Signed, typename WordType = uint64_t> struct BigInt { +private: static_assert(cpp::is_integral_v<WordType> && cpp::is_unsigned_v<WordType>, "WordType must be unsigned integer."); + struct Division { + BigInt quotient; + BigInt remainder; + }; + +public: using word_type = WordType; + using unsigned_type = BigInt<Bits, false, word_type>; + using signed_type = BigInt<Bits, true, word_type>; + LIBC_INLINE_VAR static constexpr bool SIGNED = Signed; LIBC_INLINE_VAR static constexpr size_t BITS = Bits; LIBC_INLINE_VAR @@ -100,10 +354,7 @@ struct BigInt { LIBC_INLINE_VAR static constexpr size_t WORD_COUNT = Bits / WORD_SIZE; - using unsigned_type = BigInt<BITS, false, word_type>; - using signed_type = BigInt<BITS, true, word_type>; - - cpp::array<WordType, WORD_COUNT> val{}; + cpp::array<WordType, WORD_COUNT> val{}; // zero initialized. LIBC_INLINE constexpr BigInt() = default; @@ -112,76 +363,67 @@ struct BigInt { template <size_t OtherBits, bool OtherSigned> LIBC_INLINE constexpr BigInt( const BigInt<OtherBits, OtherSigned, WordType> &other) { - if (OtherBits >= Bits) { + if (OtherBits >= Bits) { // truncate for (size_t i = 0; i < WORD_COUNT; ++i) val[i] = other[i]; - } else { + } else { // zero or sign extend size_t i = 0; for (; i < OtherBits / WORD_SIZE; ++i) val[i] = other[i]; - WordType sign = 0; - if constexpr (Signed && OtherSigned) { - sign = static_cast<WordType>( - -static_cast<cpp::make_signed_t<WordType>>(other.is_neg())); - } - for (; i < WORD_COUNT; ++i) - val[i] = sign; + extend(i, Signed && other.is_neg()); } } // Construct a BigInt from a C array. - template <size_t N, cpp::enable_if_t<N <= WORD_COUNT, int> = 0> - LIBC_INLINE constexpr BigInt(const WordType (&nums)[N]) { - size_t min_wordcount = N < WORD_COUNT ? N : WORD_COUNT; - size_t i = 0; - for (; i < min_wordcount; ++i) + template <size_t N> LIBC_INLINE constexpr BigInt(const WordType (&nums)[N]) { + static_assert(N == WORD_COUNT); + for (size_t i = 0; i < WORD_COUNT; ++i) val[i] = nums[i]; + } - // If nums doesn't completely fill val, then fill the rest with zeroes. - for (; i < WORD_COUNT; ++i) - val[i] = 0; + LIBC_INLINE constexpr explicit BigInt( + const cpp::array<WordType, WORD_COUNT> &words) { + val = words; } // Initialize the first word to |v| and the rest to 0. template <typename T, typename = cpp::enable_if_t<cpp::is_integral_v<T>>> LIBC_INLINE constexpr BigInt(T v) { - val[0] = static_cast<WordType>(v); - - if constexpr (WORD_COUNT == 1) - return; - - if constexpr (Bits < sizeof(T) * CHAR_BIT) { - for (int i = 1; i < WORD_COUNT; ++i) { - v >>= WORD_SIZE; - val[i] = static_cast<WordType>(v); + constexpr size_t T_SIZE = sizeof(T) * CHAR_BIT; + const bool is_neg = Signed && (v < 0); + for (size_t i = 0; i < WORD_COUNT; ++i) { + if (v == 0) { + extend(i, is_neg); + return; } - return; - } - - size_t i = 1; - - if constexpr (WORD_SIZE < sizeof(T) * CHAR_BIT) - for (; i < sizeof(T) * CHAR_BIT / WORD_SIZE; ++i) { + val[i] = static_cast<WordType>(v); + if constexpr (T_SIZE > WORD_SIZE) v >>= WORD_SIZE; - val[i] = static_cast<WordType>(v); - } - - WordType sign = (Signed && (v < 0)) ? ~WordType(0) : WordType(0); - for (; i < WORD_COUNT; ++i) { - val[i] = sign; + else + v = 0; } } + LIBC_INLINE constexpr BigInt &operator=(const BigInt &other) = default; - LIBC_INLINE constexpr explicit BigInt( - const cpp::array<WordType, WORD_COUNT> &words) { - for (size_t i = 0; i < WORD_COUNT; ++i) - val[i] = words[i]; + // constants + LIBC_INLINE static constexpr BigInt zero() { return BigInt(); } + LIBC_INLINE static constexpr BigInt one() { return BigInt(1); } + LIBC_INLINE static constexpr BigInt all_ones() { return ~zero(); } + LIBC_INLINE static constexpr BigInt min() { + BigInt out; + if constexpr (SIGNED) + out.set_msb(); + return out; + } + LIBC_INLINE static constexpr BigInt max() { + BigInt out = all_ones(); + if constexpr (SIGNED) + out.clear_msb(); + return out; } // TODO: Reuse the Sign type. - LIBC_INLINE constexpr bool is_neg() const { - return val.back() >> (WORD_SIZE - 1); - } + LIBC_INLINE constexpr bool is_neg() const { return SIGNED && get_msb(); } template <typename T> LIBC_INLINE constexpr explicit operator T() const { return to<T>(); @@ -191,200 +433,100 @@ struct BigInt { LIBC_INLINE constexpr cpp::enable_if_t< cpp::is_integral_v<T> && !cpp::is_same_v<T, bool>, T> to() const { + constexpr size_t T_SIZE = sizeof(T) * CHAR_BIT; T lo = static_cast<T>(val[0]); - - constexpr size_t T_BITS = sizeof(T) * CHAR_BIT; - - if constexpr (T_BITS <= WORD_SIZE) + if constexpr (T_SIZE <= WORD_SIZE) return lo; - constexpr size_t MAX_COUNT = - T_BITS > Bits ? WORD_COUNT : T_BITS / WORD_SIZE; + T_SIZE > Bits ? WORD_COUNT : T_SIZE / WORD_SIZE; for (size_t i = 1; i < MAX_COUNT; ++i) lo += static_cast<T>(val[i]) << (WORD_SIZE * i); - - if constexpr (Signed && (T_BITS > Bits)) { + if constexpr (Signed && (T_SIZE > Bits)) { // Extend sign for negative numbers. constexpr T MASK = (~T(0) << Bits); if (is_neg()) lo |= MASK; } - return lo; } LIBC_INLINE constexpr explicit operator bool() const { return !is_zero(); } - LIBC_INLINE constexpr BigInt &operator=(const BigInt &other) = default; - LIBC_INLINE constexpr bool is_zero() const { - for (size_t i = 0; i < WORD_COUNT; ++i) { - if (val[i] != 0) + for (auto part : val) + if (part != 0) return false; - } return true; } - // Add x to this number and store the result in this number. + // Add 'rhs' to this number and store the result in this number. // Returns the carry value produced by the addition operation. - LIBC_INLINE constexpr WordType add(const BigInt &x) { - SumCarry<WordType> s{0, 0}; - for (size_t i = 0; i < WORD_COUNT; ++i) { - s = add_with_carry(val[i], x.val[i], s.carry); - val[i] = s.sum; - } - return s.carry; + LIBC_INLINE constexpr WordType add_overflow(const BigInt &rhs) { + return multiword::add_with_carry(val, rhs.val); } LIBC_INLINE constexpr BigInt operator+(const BigInt &other) const { - BigInt result; - SumCarry<WordType> s{0, 0}; - for (size_t i = 0; i < WORD_COUNT; ++i) { - s = add_with_carry(val[i], other.val[i], s.carry); - result.val[i] = s.sum; - } + BigInt result = *this; + result.add_overflow(other); return result; } // This will only apply when initializing a variable from constant values, so // it will always use the constexpr version of add_with_carry. LIBC_INLINE constexpr BigInt operator+(BigInt &&other) const { - BigInt result; - SumCarry<WordType> s{0, 0}; - for (size_t i = 0; i < WORD_COUNT; ++i) { - s = add_with_carry(val[i], other.val[i], s.carry); - result.val[i] = s.sum; - } - return result; + // We use addition commutativity to reuse 'other' and prevent allocation. + other.add_overflow(*this); // Returned carry value is ignored. + return other; } LIBC_INLINE constexpr BigInt &operator+=(const BigInt &other) { - add(other); // Returned carry value is ignored. + add_overflow(other); // Returned carry value is ignored. return *this; } - // Subtract x to this number and store the result in this number. + // Subtract 'rhs' to this number and store the result in this number. // Returns the carry value produced by the subtraction operation. - LIBC_INLINE constexpr WordType sub(const BigInt &x) { - DiffBorrow<WordType> d{0, 0}; - for (size_t i = 0; i < WORD_COUNT; ++i) { - d = sub_with_borrow(val[i], x.val[i], d.borrow); - val[i] = d.diff; - } - return d.borrow; + LIBC_INLINE constexpr WordType sub_overflow(const BigInt &rhs) { + return multiword::sub_with_borrow(val, rhs.val); } LIBC_INLINE constexpr BigInt operator-(const BigInt &other) const { - BigInt result; - DiffBorrow<WordType> d{0, 0}; - for (size_t i = 0; i < WORD_COUNT; ++i) { - d = sub_with_borrow(val[i], other.val[i], d.borrow); - result.val[i] = d.diff; - } + BigInt result = *this; + result.sub_overflow(other); // Returned carry value is ignored. return result; } LIBC_INLINE constexpr BigInt operator-(BigInt &&other) const { - BigInt result; - DiffBorrow<WordType> d{0, 0}; - for (size_t i = 0; i < WORD_COUNT; ++i) { - d = sub_with_borrow(val[i], other.val[i], d.borrow); - result.val[i] = d.diff; - } + BigInt result = *this; + result.sub_overflow(other); // Returned carry value is ignored. return result; } LIBC_INLINE constexpr BigInt &operator-=(const BigInt &other) { // TODO(lntue): Set overflow flag / errno when carry is true. - sub(other); + sub_overflow(other); // Returned carry value is ignored. return *this; } - // Multiply this number with x and store the result in this number. It is - // implemented using the long multiplication algorithm by splitting the - // 64-bit words of this number and |x| in to 32-bit halves but peforming - // the operations using 64-bit numbers. This ensures that we don't lose the - // carry bits. - // Returns the carry value produced by the multiplication operation. + // Multiply this number with x and store the result in this number. LIBC_INLINE constexpr WordType mul(WordType x) { - BigInt<2 * WORD_SIZE, Signed, WordType> partial_sum(0); - for (size_t i = 0; i < WORD_COUNT; ++i) { - NumberPair<WordType> prod = internal::full_mul(val[i], x); - BigInt<2 * WORD_SIZE, Signed, WordType> tmp({prod.lo, prod.hi}); - const WordType carry = partial_sum.add(tmp); - val[i] = partial_sum.val[0]; - partial_sum.val[0] = partial_sum.val[1]; - partial_sum.val[1] = carry; - } - return partial_sum.val[1]; - } - - LIBC_INLINE constexpr BigInt operator*(const BigInt &other) const { - if constexpr (Signed) { - BigInt<Bits, false, WordType> a(*this); - BigInt<Bits, false, WordType> b(other); - const bool a_neg = a.is_neg(); - const bool b_neg = b.is_neg(); - if (a_neg) - a = -a; - if (b_neg) - b = -b; - BigInt<Bits, false, WordType> prod = a * b; - if (a_neg != b_neg) - prod = -prod; - return static_cast<BigInt<Bits, true, WordType>>(prod); - } else { - if constexpr (WORD_COUNT == 1) { - return {val[0] * other.val[0]}; - } else { - BigInt result(0); - BigInt<2 * WORD_SIZE, Signed, WordType> partial_sum(0); - WordType carry = 0; - for (size_t i = 0; i < WORD_COUNT; ++i) { - for (size_t j = 0; j <= i; j++) { - NumberPair<WordType> prod = - internal::full_mul(val[j], other.val[i - j]); - BigInt<2 * WORD_SIZE, Signed, WordType> tmp({prod.lo, prod.hi}); - carry += partial_sum.add(tmp); - } - result.val[i] = partial_sum.val[0]; - partial_sum.val[0] = partial_sum.val[1]; - partial_sum.val[1] = carry; - carry = 0; - } - return result; - } - } + return multiword::scalar_multiply_with_carry(val, x); } - // Return the full product, only unsigned for now. + // Return the full product. template <size_t OtherBits> - LIBC_INLINE constexpr BigInt<Bits + OtherBits, Signed, WordType> + LIBC_INLINE constexpr auto ful_mul(const BigInt<OtherBits, Signed, WordType> &other) const { - BigInt<Bits + OtherBits, Signed, WordType> result(0); - BigInt<2 * WORD_SIZE, Signed, WordType> partial_sum(0); - WordType carry = 0; - constexpr size_t OTHER_WORDCOUNT = - BigInt<OtherBits, Signed, WordType>::WORD_COUNT; - for (size_t i = 0; i <= WORD_COUNT + OTHER_WORDCOUNT - 2; ++i) { - const size_t lower_idx = - i < OTHER_WORDCOUNT ? 0 : i - OTHER_WORDCOUNT + 1; - const size_t upper_idx = i < WORD_COUNT ? i : WORD_COUNT - 1; - for (size_t j = lower_idx; j <= upper_idx; ++j) { - NumberPair<WordType> prod = - internal::full_mul(val[j], other.val[i - j]); - BigInt<2 * WORD_SIZE, Signed, WordType> tmp({prod.lo, prod.hi}); - carry += partial_sum.add(tmp); - } - result.val[i] = partial_sum.val[0]; - partial_sum.val[0] = partial_sum.val[1]; - partial_sum.val[1] = carry; - carry = 0; - } - result.val[WORD_COUNT + OTHER_WORDCOUNT - 1] = partial_sum.val[0]; + BigInt<Bits + OtherBits, Signed, WordType> result; + multiword::multiply_with_carry(result.val, val, other.val); return result; } + LIBC_INLINE constexpr BigInt operator*(const BigInt &other) const { + // Perform full mul and truncate. + return BigInt(ful_mul(other)); + } + // Fast hi part of the full product. The normal product `operator*` returns // `Bits` least significant bits of the full product, while this function will // approximate `Bits` most significant bits of the full product with errors @@ -407,39 +549,17 @@ struct BigInt { // 256 4 16 10 3 // 512 8 64 36 7 LIBC_INLINE constexpr BigInt quick_mul_hi(const BigInt &other) const { - BigInt result(0); - BigInt<2 * WORD_SIZE, Signed, WordType> partial_sum(0); - WordType carry = 0; - // First round of accumulation for those at WORD_COUNT - 1 in the full - // product. - for (size_t i = 0; i < WORD_COUNT; ++i) { - NumberPair<WordType> prod = - internal::full_mul(val[i], other.val[WORD_COUNT - 1 - i]); - BigInt<2 * WORD_SIZE, Signed, WordType> tmp({prod.lo, prod.hi}); - carry += partial_sum.add(tmp); - } - for (size_t i = WORD_COUNT; i < 2 * WORD_COUNT - 1; ++i) { - partial_sum.val[0] = partial_sum.val[1]; - partial_sum.val[1] = carry; - carry = 0; - for (size_t j = i - WORD_COUNT + 1; j < WORD_COUNT; ++j) { - NumberPair<WordType> prod = - internal::full_mul(val[j], other.val[i - j]); - BigInt<2 * WORD_SIZE, Signed, WordType> tmp({prod.lo, prod.hi}); - carry += partial_sum.add(tmp); - } - result.val[i - WORD_COUNT] = partial_sum.val[0]; - } - result.val[WORD_COUNT - 1] = partial_sum.val[1]; + BigInt result; + multiword::quick_mul_hi(result.val, val, other.val); return result; } - // pow takes a power and sets this to its starting value to that power. Zero - // to the zeroth power returns 1. + // BigInt(x).pow_n(n) computes x ^ n. + // Note 0 ^ 0 == 1. LIBC_INLINE constexpr void pow_n(uint64_t power) { - BigInt result = 1; + static_assert(!Signed); + BigInt result = one(); BigInt cur_power = *this; - while (power > 0) { if ((power % 2) > 0) result *= cur_power; @@ -449,38 +569,23 @@ struct BigInt { *this = result; } - // TODO: Make division work correctly for signed integers. - - // div takes another BigInt of the same size and divides this by it. The value - // of this will be set to the quotient, and the return value is the remainder. - LIBC_INLINE constexpr cpp::optional<BigInt> div(const BigInt &other) { - BigInt remainder(0); - if (*this < other) { - remainder = *this; - *this = BigInt(0); - return remainder; - } - if (other == 1) { - return remainder; - } - if (other == 0) { + // Performs inplace signed / unsigned division. Returns remainder if not + // dividing by zero. + // For signed numbers it behaves like C++ signed integer division. + // That is by truncating the fractionnal part + // https://stackoverflow.com/a/3602857 + LIBC_INLINE constexpr cpp::optional<BigInt> div(const BigInt ÷r) { + if (LIBC_UNLIKELY(divider.is_zero())) return cpp::nullopt; - } - - BigInt quotient(0); - BigInt subtractor = other; - int cur_bit = static_cast<int>(subtractor.clz() - this->clz()); - subtractor.shift_left(cur_bit); - - for (; cur_bit >= 0 && *this > 0; --cur_bit, subtractor.shift_right(1)) { - if (*this >= subtractor) { - this->sub(subtractor); - quotient = quotient | (BigInt(1) << cur_bit); - } - } - remainder = *this; - *this = quotient; - return remainder; + if (LIBC_UNLIKELY(divider == BigInt::one())) + return BigInt::zero(); + Division result; + if constexpr (SIGNED) + result = divide_signed(*this, divider); + else + result = divide_unsigned(*this, divider); + *this = result.quotient; + return result.remainder; } // Efficiently perform BigInt / (x * 2^e), where x is a half-word-size @@ -496,19 +601,16 @@ struct BigInt { // computation of each step is now properly contained within WordType. // And finally we perform some extra alignment steps for the remaining bits. LIBC_INLINE constexpr cpp::optional<BigInt> - div_uint_half_times_pow_2(internal::half_width_t<WordType> x, size_t e) { - BigInt remainder(0); - - if (x == 0) { + div_uint_half_times_pow_2(multiword::half_width_t<WordType> x, size_t e) { + BigInt remainder; + if (x == 0) return cpp::nullopt; - } if (e >= Bits) { remainder = *this; - *this = BigInt<Bits, false, WordType>(0); + *this = BigInt<Bits, false, WordType>(); return remainder; } - - BigInt quotient(0); + BigInt quotient; WordType x_word = static_cast<WordType>(x); constexpr size_t LOG2_WORD_SIZE = cpp::bit_width(WORD_SIZE) - 1; constexpr size_t HALF_WORD_SIZE = WORD_SIZE >> 1; @@ -633,189 +735,22 @@ struct BigInt { return *this; } - // TODO: remove and use cpp::countl_zero below. - [[nodiscard]] LIBC_INLINE constexpr int clz() const { - constexpr int word_digits = cpp::numeric_limits<word_type>::digits; - int leading_zeroes = 0; - for (auto i = val.size(); i > 0;) { - --i; - const int zeroes = cpp::countl_zero(val[i]); - leading_zeroes += zeroes; - if (zeroes != word_digits) - break; - } - return leading_zeroes; - } - - // TODO: remove and use cpp::countr_zero below. - [[nodiscard]] LIBC_INLINE constexpr int ctz() const { - constexpr int word_digits = cpp::numeric_limits<word_type>::digits; - int trailing_zeroes = 0; - for (auto word : val) { - const int zeroes = cpp::countr_zero(word); - trailing_zeroes += zeroes; - if (zeroes != word_digits) - break; - } - return trailing_zeroes; - } - - LIBC_INLINE constexpr void shift_left(size_t s) { - if constexpr (Bits == WORD_SIZE) { - // Use native types if possible. - if (s >= WORD_SIZE) { - val[0] = 0; - return; - } - val[0] <<= s; - return; - } - if constexpr ((Bits == 64) && (WORD_SIZE == 32)) { - // Use builtin 64 bits for 32-bit base type if available; - if (s >= 64) { - val[0] = 0; - val[1] = 0; - return; - } - uint64_t tmp = uint64__t(val[0]) + (uint64_t(val[1]) << 62); - tmp <<= s; - val[0] = uint32_t(tmp); - val[1] = uint32_t(tmp >> 32); - return; - } -#ifdef LIBC_TYPES_HAS_INT128 - if constexpr ((Bits == 128) && (WORD_SIZE == 64)) { - // Use builtin 128 bits if available; - if (s >= 128) { - val[0] = 0; - val[1] = 0; - return; - } - __uint128_t tmp = __uint128_t(val[0]) + (__uint128_t(val[1]) << 64); - tmp <<= s; - val[0] = uint64_t(tmp); - val[1] = uint64_t(tmp >> 64); - return; - } -#endif // LIBC_TYPES_HAS_INT128 - if (LIBC_UNLIKELY(s == 0)) - return; - - const size_t drop = s / WORD_SIZE; // Number of words to drop - const size_t shift = s % WORD_SIZE; // Bits to shift in the remaining words. - size_t i = WORD_COUNT; - - if (drop < WORD_COUNT) { - i = WORD_COUNT - 1; - if (shift > 0) { - for (size_t j = WORD_COUNT - 1 - drop; j > 0; --i, --j) { - val[i] = (val[j] << shift) | (val[j - 1] >> (WORD_SIZE - shift)); - } - val[i] = val[0] << shift; - } else { - for (size_t j = WORD_COUNT - 1 - drop; j > 0; --i, --j) { - val[i] = val[j]; - } - val[i] = val[0]; - } - } - - for (size_t j = 0; j < i; ++j) { - val[j] = 0; - } + LIBC_INLINE constexpr BigInt &operator<<=(size_t s) { + val = multiword::shift<multiword::LEFT, SIGNED>(val, s); + return *this; } LIBC_INLINE constexpr BigInt operator<<(size_t s) const { - BigInt result(*this); - result.shift_left(s); - return result; + return BigInt(multiword::shift<multiword::LEFT, SIGNED>(val, s)); } - LIBC_INLINE constexpr BigInt &operator<<=(size_t s) { - shift_left(s); + LIBC_INLINE constexpr BigInt &operator>>=(size_t s) { + val = multiword::shift<multiword::RIGHT, SIGNED>(val, s); return *this; } - LIBC_INLINE constexpr void shift_right(size_t s) { - if constexpr ((Bits == 64) && (WORD_SIZE == 32)) { - // Use builtin 64 bits if available; - if (s >= 64) { - val[0] = 0; - val[1] = 0; - return; - } - uint64_t tmp = uint64_t(val[0]) + (uint64_t(val[1]) << 32); - if constexpr (Signed) { - tmp = static_cast<uint64_t>(static_cast<int64_t>(tmp) >> s); - } else { - tmp >>= s; - } - val[0] = uint32_t(tmp); - val[1] = uint32_t(tmp >> 32); - return; - } -#ifdef LIBC_TYPES_HAS_INT128 - if constexpr ((Bits == 128) && (WORD_SIZE == 64)) { - // Use builtin 128 bits if available; - if (s >= 128) { - val[0] = 0; - val[1] = 0; - return; - } - __uint128_t tmp = __uint128_t(val[0]) + (__uint128_t(val[1]) << 64); - if constexpr (Signed) { - tmp = static_cast<__uint128_t>(static_cast<__int128_t>(tmp) >> s); - } else { - tmp >>= s; - } - val[0] = uint64_t(tmp); - val[1] = uint64_t(tmp >> 64); - return; - } -#endif // LIBC_TYPES_HAS_INT128 - - if (LIBC_UNLIKELY(s == 0)) - return; - const size_t drop = s / WORD_SIZE; // Number of words to drop - const size_t shift = s % WORD_SIZE; // Bit shift in the remaining words. - - size_t i = 0; - WordType sign = Signed ? is_neg() : 0; - - if (drop < WORD_COUNT) { - if (shift > 0) { - for (size_t j = drop; j < WORD_COUNT - 1; ++i, ++j) { - val[i] = (val[j] >> shift) | (val[j + 1] << (WORD_SIZE - shift)); - } - if constexpr (Signed) { - val[i] = static_cast<WordType>( - static_cast<cpp::make_signed_t<WordType>>(val[WORD_COUNT - 1]) >> - shift); - } else { - val[i] = val[WORD_COUNT - 1] >> shift; - } - ++i; - } else { - for (size_t j = drop; j < WORD_COUNT; ++i, ++j) { - val[i] = val[j]; - } - } - } - - for (; i < WORD_COUNT; ++i) { - val[i] = sign; - } - } - LIBC_INLINE constexpr BigInt operator>>(size_t s) const { - BigInt result(*this); - result.shift_right(s); - return result; - } - - LIBC_INLINE constexpr BigInt &operator>>=(size_t s) { - shift_right(s); - return *this; + return BigInt(multiword::shift<multiword::RIGHT, SIGNED>(val, s)); } #define DEFINE_BINOP(OP) \ @@ -833,10 +768,9 @@ struct BigInt { return lhs; \ } - DEFINE_BINOP(&) - DEFINE_BINOP(|) - DEFINE_BINOP(^) - + DEFINE_BINOP(&) // & and &= + DEFINE_BINOP(|) // | and |= + DEFINE_BINOP(^) // ^ and ^= #undef DEFINE_BINOP LIBC_INLINE constexpr BigInt operator~() const { @@ -847,8 +781,8 @@ struct BigInt { } LIBC_INLINE constexpr BigInt operator-() const { - BigInt result = ~(*this); - result.add(BigInt(1)); + BigInt result(*this); + result.negate(); return result; } @@ -865,24 +799,6 @@ struct BigInt { return !(lhs == rhs); } -private: - LIBC_INLINE friend constexpr int cmp(const BigInt &lhs, const BigInt &rhs) { - const auto compare = [](WordType a, WordType b) { - return a == b ? 0 : a > b ? 1 : -1; - }; - if constexpr (Signed) { - const bool lhs_is_neg = lhs.is_neg(); - const bool rhs_is_neg = rhs.is_neg(); - if (lhs_is_neg != rhs_is_neg) - return rhs_is_neg ? 1 : -1; - } - for (size_t i = WORD_COUNT; i-- > 0;) - if (auto cmp = compare(lhs[i], rhs[i]); cmp != 0) - return cmp; - return 0; - } - -public: LIBC_INLINE friend constexpr bool operator>(const BigInt &lhs, const BigInt &rhs) { return cmp(lhs, rhs) > 0; @@ -901,24 +817,24 @@ public: } LIBC_INLINE constexpr BigInt &operator++() { - add(BigInt(1)); + increment(); return *this; } LIBC_INLINE constexpr BigInt operator++(int) { BigInt oldval(*this); - add(BigInt(1)); + increment(); return oldval; } LIBC_INLINE constexpr BigInt &operator--() { - sub(BigInt(1)); + decrement(); return *this; } LIBC_INLINE constexpr BigInt operator--(int) { BigInt oldval(*this); - sub(BigInt(1)); + decrement(); return oldval; } @@ -930,9 +846,117 @@ public: // Return the i-th word of the number. LIBC_INLINE constexpr WordType &operator[](size_t i) { return val[i]; } - LIBC_INLINE WordType *data() { return val; } +private: + LIBC_INLINE friend constexpr int cmp(const BigInt &lhs, const BigInt &rhs) { + constexpr auto compare = [](WordType a, WordType b) { + return a == b ? 0 : a > b ? 1 : -1; + }; + if constexpr (Signed) { + const bool lhs_is_neg = lhs.is_neg(); + const bool rhs_is_neg = rhs.is_neg(); + if (lhs_is_neg != rhs_is_neg) + return rhs_is_neg ? 1 : -1; + } + for (size_t i = WORD_COUNT; i-- > 0;) + if (auto cmp = compare(lhs[i], rhs[i]); cmp != 0) + return cmp; + return 0; + } + + LIBC_INLINE constexpr void bitwise_not() { + for (auto &part : val) + part = ~part; + } + + LIBC_INLINE constexpr void negate() { + bitwise_not(); + increment(); + } - LIBC_INLINE const WordType *data() const { return val; } + LIBC_INLINE constexpr void increment() { + multiword::add_with_carry(val, cpp::array<WordType, 1>{1}); + } + + LIBC_INLINE constexpr void decrement() { + multiword::add_with_carry(val, cpp::array<WordType, 1>{1}); + } + + LIBC_INLINE constexpr void extend(size_t index, bool is_neg) { + const WordType value = is_neg ? cpp::numeric_limits<WordType>::max() + : cpp::numeric_limits<WordType>::min(); + for (size_t i = index; i < WORD_COUNT; ++i) + val[i] = value; + } + + LIBC_INLINE constexpr bool get_msb() const { + return val.back() >> (WORD_SIZE - 1); + } + + LIBC_INLINE constexpr void set_msb() { + val.back() |= mask_leading_ones<WordType, 1>(); + } + + LIBC_INLINE constexpr void clear_msb() { + val.back() &= mask_trailing_ones<WordType, WORD_SIZE - 1>(); + } + + LIBC_INLINE constexpr void set_bit(size_t i) { + const size_t word_index = i / WORD_SIZE; + val[word_index] |= WordType(1) << (i % WORD_SIZE); + } + + LIBC_INLINE constexpr static Division divide_unsigned(const BigInt ÷nd, + const BigInt ÷r) { + BigInt remainder = dividend; + BigInt quotient; + if (remainder >= divider) { + BigInt subtractor = divider; + int cur_bit = multiword::countl_zero(subtractor.val) - + multiword::countl_zero(remainder.val); + subtractor <<= cur_bit; + for (; cur_bit >= 0 && remainder > 0; --cur_bit, subtractor >>= 1) { + if (remainder < subtractor) + continue; + remainder -= subtractor; + quotient.set_bit(cur_bit); + } + } + return Division{quotient, remainder}; + } + + LIBC_INLINE constexpr static Division divide_signed(const BigInt ÷nd, + const BigInt ÷r) { + // Special case because it is not possible to negate the min value of a + // signed integer. + if (dividend == min() && divider == min()) + return Division{one(), zero()}; + // 1. Convert the dividend and divisor to unsigned representation. + unsigned_type udividend(dividend); + unsigned_type udivider(divider); + // 2. Negate the dividend if it's negative, and similarly for the divisor. + const bool dividend_is_neg = dividend.is_neg(); + const bool divider_is_neg = divider.is_neg(); + if (dividend_is_neg) + udividend.negate(); + if (divider_is_neg) + udivider.negate(); + // 3. Use unsigned multiword division algorithm. + const auto unsigned_result = divide_unsigned(udividend, udivider); + // 4. Convert the quotient and remainder to signed representation. + Division result; + result.quotient = signed_type(unsigned_result.quotient); + result.remainder = signed_type(unsigned_result.remainder); + // 5. Negate the quotient if the dividend and divisor had opposite signs. + if (dividend_is_neg != divider_is_neg) + result.quotient.negate(); + // 6. Negate the remainder if the dividend was negative. + if (dividend_is_neg) + result.remainder.negate(); + return result; + } + + friend signed_type; + friend unsigned_type; }; namespace internal { @@ -962,10 +986,8 @@ using Int = BigInt<Bits, true, internal::WordTypeSelectorT<Bits>>; // Provides limits of U/Int<128>. template <> class cpp::numeric_limits<UInt<128>> { public: - LIBC_INLINE static constexpr UInt<128> max() { - return UInt<128>({0xffff'ffff'ffff'ffff, 0xffff'ffff'ffff'ffff}); - } - LIBC_INLINE static constexpr UInt<128> min() { return UInt<128>(0); } + LIBC_INLINE static constexpr UInt<128> max() { return UInt<128>::max(); } + LIBC_INLINE static constexpr UInt<128> min() { return UInt<128>::min(); } // Meant to match std::numeric_limits interface. // NOLINTNEXTLINE(readability-identifier-naming) LIBC_INLINE_VAR static constexpr int digits = 128; @@ -973,12 +995,8 @@ public: template <> class cpp::numeric_limits<Int<128>> { public: - LIBC_INLINE static constexpr Int<128> max() { - return Int<128>({0xffff'ffff'ffff'ffff, 0x7fff'ffff'ffff'ffff}); - } - LIBC_INLINE static constexpr Int<128> min() { - return Int<128>({0, 0x8000'0000'0000'0000}); - } + LIBC_INLINE static constexpr Int<128> max() { return Int<128>::max(); } + LIBC_INLINE static constexpr Int<128> min() { return Int<128>::min(); } // Meant to match std::numeric_limits interface. // NOLINTNEXTLINE(readability-identifier-naming) LIBC_INLINE_VAR static constexpr int digits = 128; @@ -1112,30 +1130,28 @@ has_single_bit(T value) { template <typename T> [[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> countr_zero(const T &value) { - return value.ctz(); + return multiword::countr_zero(value.val); } // Specialization of cpp::countl_zero ('bit.h') for BigInt. template <typename T> [[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> countl_zero(const T &value) { - return value.clz(); + return multiword::countl_zero(value.val); } // Specialization of cpp::countl_one ('bit.h') for BigInt. template <typename T> [[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> countl_one(T value) { - // TODO : Implement a faster version not involving operator~. - return cpp::countl_zero<T>(~value); + return multiword::countl_one(value.val); } // Specialization of cpp::countr_one ('bit.h') for BigInt. template <typename T> [[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> countr_one(T value) { - // TODO : Implement a faster version not involving operator~. - return cpp::countr_zero<T>(~value); + return multiword::countr_one(value.val); } // Specialization of cpp::bit_width ('bit.h') for BigInt. @@ -1182,65 +1198,59 @@ rotr(T value, int rotate) { template <typename T, size_t count> LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, T> mask_trailing_ones() { - static_assert(!T::SIGNED); - if (count == 0) - return T(); - constexpr unsigned T_BITS = CHAR_BIT * sizeof(T); - static_assert(count <= T_BITS && "Invalid bit index"); - using word_type = typename T::word_type; - T out; - constexpr int CHUNK_INDEX_CONTAINING_BIT = - static_cast<int>(count / T::WORD_SIZE); - int index = 0; - for (auto &word : out.val) { - if (index < CHUNK_INDEX_CONTAINING_BIT) - word = -1; - else if (index > CHUNK_INDEX_CONTAINING_BIT) - word = 0; - else - word = mask_trailing_ones<word_type, count % T::WORD_SIZE>(); - ++index; - } + static_assert(!T::SIGNED && count <= T::BITS); + if (count == T::BITS) + return T::all_ones(); + constexpr size_t QUOTIENT = count / T::WORD_SIZE; + constexpr size_t REMAINDER = count % T::WORD_SIZE; + T out; // zero initialized + for (size_t i = 0; i <= QUOTIENT; ++i) + out[i] = i < QUOTIENT + ? -1 + : mask_trailing_ones<typename T::word_type, REMAINDER>(); return out; } // Specialization of mask_leading_ones ('math_extras.h') for BigInt. template <typename T, size_t count> LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, T> mask_leading_ones() { - static_assert(!T::SIGNED); - if (count == 0) - return T(); - constexpr unsigned T_BITS = CHAR_BIT * sizeof(T); - static_assert(count <= T_BITS && "Invalid bit index"); - using word_type = typename T::word_type; - T out; - constexpr int CHUNK_INDEX_CONTAINING_BIT = - static_cast<int>((T::BITS - count - 1ULL) / T::WORD_SIZE); - int index = 0; - for (auto &word : out.val) { - if (index < CHUNK_INDEX_CONTAINING_BIT) - word = 0; - else if (index > CHUNK_INDEX_CONTAINING_BIT) - word = -1; - else - word = mask_leading_ones<word_type, count % T::WORD_SIZE>(); - ++index; - } + static_assert(!T::SIGNED && count <= T::BITS); + if (count == T::BITS) + return T::all_ones(); + constexpr size_t QUOTIENT = (T::BITS - count - 1U) / T::WORD_SIZE; + constexpr size_t REMAINDER = count % T::WORD_SIZE; + T out; // zero initialized + for (size_t i = QUOTIENT; i < T::WORD_COUNT; ++i) + out[i] = i > QUOTIENT + ? -1 + : mask_leading_ones<typename T::word_type, REMAINDER>(); return out; } +// Specialization of mask_trailing_zeros ('math_extras.h') for BigInt. +template <typename T, size_t count> +LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, T> +mask_trailing_zeros() { + return mask_leading_ones<T, T::BITS - count>(); +} + +// Specialization of mask_leading_zeros ('math_extras.h') for BigInt. +template <typename T, size_t count> +LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, T> +mask_leading_zeros() { + return mask_trailing_ones<T, T::BITS - count>(); +} + // Specialization of count_zeros ('math_extras.h') for BigInt. template <typename T> -[[nodiscard]] -LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> +[[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> count_zeros(T value) { return cpp::popcount(~value); } // Specialization of first_leading_zero ('math_extras.h') for BigInt. template <typename T> -[[nodiscard]] -LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> +[[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> first_leading_zero(T value) { return value == cpp::numeric_limits<T>::max() ? 0 : cpp::countl_one(value) + 1; @@ -1248,16 +1258,14 @@ first_leading_zero(T value) { // Specialization of first_leading_one ('math_extras.h') for BigInt. template <typename T> -[[nodiscard]] -LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> +[[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> first_leading_one(T value) { return first_leading_zero(~value); } // Specialization of first_trailing_zero ('math_extras.h') for BigInt. template <typename T> -[[nodiscard]] -LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> +[[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> first_trailing_zero(T value) { return value == cpp::numeric_limits<T>::max() ? 0 : cpp::countr_zero(~value) + 1; @@ -1265,8 +1273,7 @@ first_trailing_zero(T value) { // Specialization of first_trailing_one ('math_extras.h') for BigInt. template <typename T> -[[nodiscard]] -LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> +[[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<is_big_int_v<T>, int> first_trailing_one(T value) { return value == cpp::numeric_limits<T>::max() ? 0 : cpp::countr_zero(value) + 1; diff --git a/libc/src/__support/float_to_string.h b/libc/src/__support/float_to_string.h index 1287c3e..4c59cfd 100644 --- a/libc/src/__support/float_to_string.h +++ b/libc/src/__support/float_to_string.h @@ -689,7 +689,7 @@ template <> class FloatToString<long double> { wide_int float_as_int = mantissa; - float_as_int.shift_left(exponent); + float_as_int <<= exponent; int_block_index = 0; while (float_as_int > 0) { @@ -708,10 +708,11 @@ template <> class FloatToString<long double> { const int SHIFT_AMOUNT = FLOAT_AS_INT_WIDTH + exponent; static_assert(EXTRA_INT_WIDTH >= sizeof(long double) * 8); - float_as_fixed.shift_left(SHIFT_AMOUNT); + float_as_fixed <<= SHIFT_AMOUNT; // If there are still digits above the decimal point, handle those. - if (float_as_fixed.clz() < static_cast<int>(EXTRA_INT_WIDTH)) { + if (cpp::countl_zero(float_as_fixed) < + static_cast<int>(EXTRA_INT_WIDTH)) { UInt<EXTRA_INT_WIDTH> above_decimal_point = float_as_fixed >> FLOAT_AS_INT_WIDTH; diff --git a/libc/src/__support/integer_literals.h b/libc/src/__support/integer_literals.h index de1f88f..e99799c 100644 --- a/libc/src/__support/integer_literals.h +++ b/libc/src/__support/integer_literals.h @@ -151,12 +151,15 @@ template <size_t N> struct Parser<LIBC_NAMESPACE::UInt<N>> { template <typename T> LIBC_INLINE constexpr T parse_with_prefix(const char *ptr) { using P = Parser<T>; - if (ptr[0] == '0' && ptr[1] == 'x') - return P::template parse<16>(ptr + 2); - else if (ptr[0] == '0' && ptr[1] == 'b') - return P::template parse<2>(ptr + 2); - else - return P::template parse<10>(ptr); + if (ptr == nullptr) + return T(); + if (ptr[0] == '0') { + if (ptr[1] == 'b') + return P::template parse<2>(ptr + 2); + if (ptr[1] == 'x') + return P::template parse<16>(ptr + 2); + } + return P::template parse<10>(ptr); } } // namespace internal @@ -169,6 +172,16 @@ LIBC_INLINE constexpr auto operator""_u256(const char *x) { return internal::parse_with_prefix<UInt<256>>(x); } +template <typename T> LIBC_INLINE constexpr T parse_bigint(const char *ptr) { + if (ptr == nullptr) + return T(); + if (ptr[0] == '-' || ptr[0] == '+') { + auto positive = internal::parse_with_prefix<T>(ptr + 1); + return ptr[0] == '-' ? -positive : positive; + } + return internal::parse_with_prefix<T>(ptr); +} + } // namespace LIBC_NAMESPACE #endif // LLVM_LIBC_SRC___SUPPORT_INTEGER_LITERALS_H diff --git a/libc/src/__support/math_extras.h b/libc/src/__support/math_extras.h index 70a8800..4bd8719 100644 --- a/libc/src/__support/math_extras.h +++ b/libc/src/__support/math_extras.h @@ -10,9 +10,9 @@ #ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXTRAS_H #define LLVM_LIBC_SRC___SUPPORT_MATH_EXTRAS_H -#include "src/__support/CPP/bit.h" // countl_one, countr_zero -#include "src/__support/CPP/limits.h" // CHAR_BIT, numeric_limits -#include "src/__support/CPP/type_traits.h" // is_unsigned_v +#include "src/__support/CPP/bit.h" // countl_one, countr_zero +#include "src/__support/CPP/limits.h" // CHAR_BIT, numeric_limits +#include "src/__support/CPP/type_traits.h" // is_unsigned_v, is_constant_evaluated #include "src/__support/macros/attributes.h" // LIBC_INLINE namespace LIBC_NAMESPACE { @@ -32,199 +32,94 @@ mask_trailing_ones() { template <typename T, size_t count> LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_unsigned_v<T>, T> mask_leading_ones() { - constexpr T MASK(mask_trailing_ones<T, CHAR_BIT * sizeof(T) - count>()); - return T(~MASK); // bitwise NOT performs integer promotion. + return T(~mask_trailing_ones<T, CHAR_BIT * sizeof(T) - count>()); } -// Add with carry -template <typename T> struct SumCarry { - T sum; - T carry; -}; - -// This version is always valid for constexpr. -template <typename T> -LIBC_INLINE constexpr cpp::enable_if_t< - cpp::is_integral_v<T> && cpp::is_unsigned_v<T>, SumCarry<T>> -add_with_carry_const(T a, T b, T carry_in) { - T tmp = a + carry_in; - T sum = b + tmp; - T carry_out = (sum < b) + (tmp < a); - return {sum, carry_out}; -} - -template <typename T> -LIBC_INLINE constexpr cpp::enable_if_t< - cpp::is_integral_v<T> && cpp::is_unsigned_v<T>, SumCarry<T>> -add_with_carry(T a, T b, T carry_in) { - return add_with_carry_const<T>(a, b, carry_in); -} - -#if __has_builtin(__builtin_addc) -// https://clang.llvm.org/docs/LanguageExtensions.html#multiprecision-arithmetic-builtins - -template <> -LIBC_INLINE constexpr SumCarry<unsigned char> -add_with_carry<unsigned char>(unsigned char a, unsigned char b, - unsigned char carry_in) { - if (__builtin_is_constant_evaluated()) { - return add_with_carry_const<unsigned char>(a, b, carry_in); - } else { - SumCarry<unsigned char> result{0, 0}; - result.sum = __builtin_addcb(a, b, carry_in, &result.carry); - return result; - } -} - -template <> -LIBC_INLINE constexpr SumCarry<unsigned short> -add_with_carry<unsigned short>(unsigned short a, unsigned short b, - unsigned short carry_in) { - if (__builtin_is_constant_evaluated()) { - return add_with_carry_const<unsigned short>(a, b, carry_in); - } else { - SumCarry<unsigned short> result{0, 0}; - result.sum = __builtin_addcs(a, b, carry_in, &result.carry); - return result; - } -} - -template <> -LIBC_INLINE constexpr SumCarry<unsigned int> -add_with_carry<unsigned int>(unsigned int a, unsigned int b, - unsigned int carry_in) { - if (__builtin_is_constant_evaluated()) { - return add_with_carry_const<unsigned int>(a, b, carry_in); - } else { - SumCarry<unsigned int> result{0, 0}; - result.sum = __builtin_addc(a, b, carry_in, &result.carry); - return result; - } -} - -template <> -LIBC_INLINE constexpr SumCarry<unsigned long> -add_with_carry<unsigned long>(unsigned long a, unsigned long b, - unsigned long carry_in) { - if (__builtin_is_constant_evaluated()) { - return add_with_carry_const<unsigned long>(a, b, carry_in); - } else { - SumCarry<unsigned long> result{0, 0}; - result.sum = __builtin_addcl(a, b, carry_in, &result.carry); - return result; - } +// Create a bitmask with the count right-most bits set to 0, and all other bits +// set to 1. Only unsigned types are allowed. +template <typename T, size_t count> +LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_unsigned_v<T>, T> +mask_trailing_zeros() { + return mask_leading_ones<T, CHAR_BIT * sizeof(T) - count>(); } -template <> -LIBC_INLINE constexpr SumCarry<unsigned long long> -add_with_carry<unsigned long long>(unsigned long long a, unsigned long long b, - unsigned long long carry_in) { - if (__builtin_is_constant_evaluated()) { - return add_with_carry_const<unsigned long long>(a, b, carry_in); - } else { - SumCarry<unsigned long long> result{0, 0}; - result.sum = __builtin_addcll(a, b, carry_in, &result.carry); - return result; - } +// Create a bitmask with the count left-most bits set to 0, and all other bits +// set to 1. Only unsigned types are allowed. +template <typename T, size_t count> +LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_unsigned_v<T>, T> +mask_leading_zeros() { + return mask_trailing_ones<T, CHAR_BIT * sizeof(T) - count>(); } -#endif // __has_builtin(__builtin_addc) - -// Subtract with borrow -template <typename T> struct DiffBorrow { - T diff; - T borrow; -}; - -// This version is always valid for constexpr. +// Returns whether 'a + b' overflows, the result is stored in 'res'. template <typename T> -LIBC_INLINE constexpr cpp::enable_if_t< - cpp::is_integral_v<T> && cpp::is_unsigned_v<T>, DiffBorrow<T>> -sub_with_borrow_const(T a, T b, T borrow_in) { - T tmp = a - b; - T diff = tmp - borrow_in; - T borrow_out = (diff > tmp) + (tmp > a); - return {diff, borrow_out}; +[[nodiscard]] LIBC_INLINE constexpr bool add_overflow(T a, T b, T &res) { + return __builtin_add_overflow(a, b, &res); } -// This version is not always valid for constepxr because it's overriden below -// if builtins are available. +// Returns whether 'a - b' overflows, the result is stored in 'res'. template <typename T> -LIBC_INLINE constexpr cpp::enable_if_t< - cpp::is_integral_v<T> && cpp::is_unsigned_v<T>, DiffBorrow<T>> -sub_with_borrow(T a, T b, T borrow_in) { - return sub_with_borrow_const<T>(a, b, borrow_in); -} - -#if __has_builtin(__builtin_subc) -// https://clang.llvm.org/docs/LanguageExtensions.html#multiprecision-arithmetic-builtins - -template <> -LIBC_INLINE constexpr DiffBorrow<unsigned char> -sub_with_borrow<unsigned char>(unsigned char a, unsigned char b, - unsigned char borrow_in) { - if (__builtin_is_constant_evaluated()) { - return sub_with_borrow_const<unsigned char>(a, b, borrow_in); - } else { - DiffBorrow<unsigned char> result{0, 0}; - result.diff = __builtin_subcb(a, b, borrow_in, &result.borrow); - return result; - } +[[nodiscard]] LIBC_INLINE constexpr bool sub_overflow(T a, T b, T &res) { + return __builtin_sub_overflow(a, b, &res); } -template <> -LIBC_INLINE constexpr DiffBorrow<unsigned short> -sub_with_borrow<unsigned short>(unsigned short a, unsigned short b, - unsigned short borrow_in) { - if (__builtin_is_constant_evaluated()) { - return sub_with_borrow_const<unsigned short>(a, b, borrow_in); - } else { - DiffBorrow<unsigned short> result{0, 0}; - result.diff = __builtin_subcs(a, b, borrow_in, &result.borrow); - return result; - } -} +#define RETURN_IF(TYPE, BUILTIN) \ + if constexpr (cpp::is_same_v<T, TYPE>) \ + return BUILTIN(a, b, carry_in, carry_out); -template <> -LIBC_INLINE constexpr DiffBorrow<unsigned int> -sub_with_borrow<unsigned int>(unsigned int a, unsigned int b, - unsigned int borrow_in) { - if (__builtin_is_constant_evaluated()) { - return sub_with_borrow_const<unsigned int>(a, b, borrow_in); - } else { - DiffBorrow<unsigned int> result{0, 0}; - result.diff = __builtin_subc(a, b, borrow_in, &result.borrow); - return result; - } -} - -template <> -LIBC_INLINE constexpr DiffBorrow<unsigned long> -sub_with_borrow<unsigned long>(unsigned long a, unsigned long b, - unsigned long borrow_in) { - if (__builtin_is_constant_evaluated()) { - return sub_with_borrow_const<unsigned long>(a, b, borrow_in); - } else { - DiffBorrow<unsigned long> result{0, 0}; - result.diff = __builtin_subcl(a, b, borrow_in, &result.borrow); - return result; +// Returns the result of 'a + b' taking into account 'carry_in'. +// The carry out is stored in 'carry_out' it not 'nullptr', dropped otherwise. +// We keep the pass by pointer interface for consistency with the intrinsic. +template <typename T> +[[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_unsigned_v<T>, T> +add_with_carry(T a, T b, T carry_in, T &carry_out) { + if constexpr (!cpp::is_constant_evaluated()) { +#if __has_builtin(__builtin_addcb) + RETURN_IF(unsigned char, __builtin_addcb) +#elif __has_builtin(__builtin_addcs) + RETURN_IF(unsigned short, __builtin_addcs) +#elif __has_builtin(__builtin_addc) + RETURN_IF(unsigned int, __builtin_addc) +#elif __has_builtin(__builtin_addcl) + RETURN_IF(unsigned long, __builtin_addcl) +#elif __has_builtin(__builtin_addcll) + RETURN_IF(unsigned long long, __builtin_addcll) +#endif } + T sum = {}; + T carry1 = add_overflow(a, b, sum); + T carry2 = add_overflow(sum, carry_in, sum); + carry_out = carry1 | carry2; + return sum; } -template <> -LIBC_INLINE constexpr DiffBorrow<unsigned long long> -sub_with_borrow<unsigned long long>(unsigned long long a, unsigned long long b, - unsigned long long borrow_in) { - if (__builtin_is_constant_evaluated()) { - return sub_with_borrow_const<unsigned long long>(a, b, borrow_in); - } else { - DiffBorrow<unsigned long long> result{0, 0}; - result.diff = __builtin_subcll(a, b, borrow_in, &result.borrow); - return result; +// Returns the result of 'a - b' taking into account 'carry_in'. +// The carry out is stored in 'carry_out' it not 'nullptr', dropped otherwise. +// We keep the pass by pointer interface for consistency with the intrinsic. +template <typename T> +[[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_unsigned_v<T>, T> +sub_with_borrow(T a, T b, T carry_in, T &carry_out) { + if constexpr (!cpp::is_constant_evaluated()) { +#if __has_builtin(__builtin_subcb) + RETURN_IF(unsigned char, __builtin_subcb) +#elif __has_builtin(__builtin_subcs) + RETURN_IF(unsigned short, __builtin_subcs) +#elif __has_builtin(__builtin_subc) + RETURN_IF(unsigned int, __builtin_subc) +#elif __has_builtin(__builtin_subcl) + RETURN_IF(unsigned long, __builtin_subcl) +#elif __has_builtin(__builtin_subcll) + RETURN_IF(unsigned long long, __builtin_subcll) +#endif } + T sub = {}; + T carry1 = sub_overflow(a, b, sub); + T carry2 = sub_overflow(sub, carry_in, sub); + carry_out = carry1 | carry2; + return sub; } -#endif // __has_builtin(__builtin_subc) +#undef RETURN_IF template <typename T> [[nodiscard]] LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_unsigned_v<T>, int> diff --git a/libc/src/__support/number_pair.h b/libc/src/__support/number_pair.h index ee6667b..2f713fc 100644 --- a/libc/src/__support/number_pair.h +++ b/libc/src/__support/number_pair.h @@ -20,17 +20,6 @@ template <typename T> struct NumberPair { T hi = T(0); }; -template <typename T> -cpp::enable_if_t<cpp::is_integral_v<T> && cpp::is_unsigned_v<T>, - NumberPair<T>> constexpr split(T a) { - constexpr size_t HALF_BIT_WIDTH = sizeof(T) * 4; - constexpr T LOWER_HALF_MASK = (T(1) << HALF_BIT_WIDTH) - T(1); - NumberPair<T> result; - result.lo = a & LOWER_HALF_MASK; - result.hi = a >> HALF_BIT_WIDTH; - return result; -} - } // namespace LIBC_NAMESPACE #endif // LLVM_LIBC_SRC___SUPPORT_NUMBER_PAIR_H diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index b67ee3a..c89792b 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -88,6 +88,8 @@ add_math_entrypoint_object(expf) add_math_entrypoint_object(exp2) add_math_entrypoint_object(exp2f) +add_math_entrypoint_object(exp2m1f) + add_math_entrypoint_object(exp10) add_math_entrypoint_object(exp10f) diff --git a/libc/src/math/exp2m1f.h b/libc/src/math/exp2m1f.h new file mode 100644 index 0000000..0eaf6b0 --- /dev/null +++ b/libc/src/math/exp2m1f.h @@ -0,0 +1,18 @@ +//===-- Implementation header for exp2m1f -----------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_EXP2M1F_H +#define LLVM_LIBC_SRC_MATH_EXP2M1F_H + +namespace LIBC_NAMESPACE { + +float exp2m1f(float x); + +} // namespace LIBC_NAMESPACE + +#endif // LLVM_LIBC_SRC_MATH_EXP2M1F_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index b164d33..dc77f8b 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -838,6 +838,27 @@ add_entrypoint_object( ) add_entrypoint_object( + exp2m1f + SRCS + exp2m1f.cpp + HDRS + ../exp2m1f.h + DEPENDS + .explogxf + libc.src.errno.errno + libc.src.__support.common + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.rounding_mode + libc.src.__support.macros.optimization + libc.src.__support.macros.properties.cpu_features + COMPILE_OPTIONS + -O3 +) + +add_entrypoint_object( exp10 SRCS exp10.cpp diff --git a/libc/src/math/generic/exp2m1f.cpp b/libc/src/math/generic/exp2m1f.cpp new file mode 100644 index 0000000..c60930d --- /dev/null +++ b/libc/src/math/generic/exp2m1f.cpp @@ -0,0 +1,183 @@ +//===-- Implementation of exp2m1f function --------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "src/math/exp2m1f.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/except_value_utils.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/rounding_mode.h" +#include "src/__support/common.h" +#include "src/__support/macros/optimization.h" +#include "src/__support/macros/properties/cpu_features.h" +#include "src/errno/libc_errno.h" + +#include "explogxf.h" + +namespace LIBC_NAMESPACE { + +static constexpr size_t N_EXCEPTS_LO = 8; + +static constexpr fputil::ExceptValues<float, N_EXCEPTS_LO> EXP2M1F_EXCEPTS_LO = + {{ + // (input, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.36dc8ep-36, exp2m1f(x) = 0x1.aef212p-37 (RZ) + {0x2d9b'6e47U, 0x2d57'7909U, 1U, 0U, 0U}, + // x = 0x1.224936p-19, exp2m1f(x) = 0x1.926c0ep-20 (RZ) + {0x3611'249bU, 0x35c9'3607U, 1U, 0U, 1U}, + // x = 0x1.d16d2p-20, exp2m1f(x) = 0x1.429becp-20 (RZ) + {0x35e8'b690U, 0x35a1'4df6U, 1U, 0U, 1U}, + // x = 0x1.17949ep-14, exp2m1f(x) = 0x1.8397p-15 (RZ) + {0x388b'ca4fU, 0x3841'cb80U, 1U, 0U, 1U}, + // x = -0x1.9c3e1ep-38, exp2m1f(x) = -0x1.1dbeacp-38 (RZ) + {0xacce'1f0fU, 0xac8e'df56U, 0U, 1U, 0U}, + // x = -0x1.4d89b4p-32, exp2m1f(x) = -0x1.ce61b6p-33 (RZ) + {0xafa6'c4daU, 0xaf67'30dbU, 0U, 1U, 1U}, + // x = -0x1.a6eac4p-10, exp2m1f(x) = -0x1.24fadap-10 (RZ) + {0xbad3'7562U, 0xba92'7d6dU, 0U, 1U, 1U}, + // x = -0x1.e7526ep-6, exp2m1f(x) = -0x1.4e53dep-6 (RZ) + {0xbcf3'a937U, 0xbca7'29efU, 0U, 1U, 1U}, + }}; + +static constexpr size_t N_EXCEPTS_HI = 3; + +static constexpr fputil::ExceptValues<float, N_EXCEPTS_HI> EXP2M1F_EXCEPTS_HI = + {{ + // (input, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.16a972p-1, exp2m1f(x) = 0x1.d545b2p-2 (RZ) + {0x3f0b'54b9U, 0x3eea'a2d9U, 1U, 0U, 0U}, + // x = -0x1.9f12acp-5, exp2m1f(x) = -0x1.1ab68cp-5 (RZ) + {0xbd4f'8956U, 0xbd0d'5b46U, 0U, 1U, 0U}, + // x = -0x1.de7b9cp-5, exp2m1f(x) = -0x1.4508f4p-5 (RZ) + {0xbd6f'3dceU, 0xbd22'847aU, 0U, 1U, 1U}, + }}; + +LLVM_LIBC_FUNCTION(float, exp2m1f, (float x)) { + using FPBits = fputil::FPBits<float>; + FPBits xbits(x); + + uint32_t x_u = xbits.uintval(); + uint32_t x_abs = x_u & 0x7fff'ffffU; + + // When |x| >= 128, or x is nan, or |x| <= 2^-5 + if (LIBC_UNLIKELY(x_abs >= 0x4300'0000U || x_abs <= 0x3d00'0000U)) { + // |x| <= 2^-5 + if (x_abs <= 0x3d00'0000U) { + if (auto r = EXP2M1F_EXCEPTS_LO.lookup(x_u); LIBC_UNLIKELY(r.has_value())) + return r.value(); + + // Minimax polynomial generated by Sollya with: + // > display = hexadecimal; + // > fpminimax((2^x - 1)/x, 5, [|D...|], [-2^-5, 2^-5]); + constexpr double COEFFS[] = { + 0x1.62e42fefa39f3p-1, 0x1.ebfbdff82c57bp-3, 0x1.c6b08d6f2d7aap-5, + 0x1.3b2ab6fc92f5dp-7, 0x1.5d897cfe27125p-10, 0x1.43090e61e6af1p-13}; + double xd = x; + double xsq = xd * xd; + double c0 = fputil::multiply_add(xd, COEFFS[1], COEFFS[0]); + double c1 = fputil::multiply_add(xd, COEFFS[3], COEFFS[2]); + double c2 = fputil::multiply_add(xd, COEFFS[5], COEFFS[4]); + double p = fputil::polyeval(xsq, c0, c1, c2); + return static_cast<float>(p * xd); + } + + // x >= 128, or x is nan + if (xbits.is_pos()) { + if (xbits.is_finite()) { + int rounding = fputil::quick_get_round(); + if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) + return FPBits::max_normal().get_val(); + + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_OVERFLOW); + } + + // x >= 128 and 2^x - 1 rounds to +inf, or x is +inf or nan + return x + FPBits::inf().get_val(); + } + } + + if (LIBC_UNLIKELY(x <= -25.0f)) { + // 2^(-inf) - 1 = -1 + if (xbits.is_inf()) + return -1.0f; + // 2^nan - 1 = nan + if (xbits.is_nan()) + return x; + + int rounding = fputil::quick_get_round(); + if (rounding == FE_UPWARD || rounding == FE_TOWARDZERO) + return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f + + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_UNDERFLOW); + return -1.0f; + } + + if (auto r = EXP2M1F_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value())) + return r.value(); + + // For -25 < x < 128, to compute 2^x, we perform the following range + // reduction: find hi, mid, lo such that: + // x = hi + mid + lo, in which: + // hi is an integer, + // 0 <= mid * 2^5 < 32 is an integer, + // -2^(-6) <= lo <= 2^(-6). + // In particular, + // hi + mid = round(x * 2^5) * 2^(-5). + // Then, + // 2^x = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo. + // 2^mid is stored in the lookup table of 32 elements. + // 2^lo is computed using a degree-4 minimax polynomial generated by Sollya. + // We perform 2^hi * 2^mid by simply add hi to the exponent field of 2^mid. + + // kf = (hi + mid) * 2^5 = round(x * 2^5) + float kf; + int k; +#ifdef LIBC_TARGET_CPU_HAS_NEAREST_INT + kf = fputil::nearest_integer(x * 32.0f); + k = static_cast<int>(kf); +#else + constexpr float HALF[2] = {0.5f, -0.5f}; + k = static_cast<int>(fputil::multiply_add(x, 32.0f, HALF[x < 0.0f])); + kf = static_cast<float>(k); +#endif // LIBC_TARGET_CPU_HAS_NEAREST_INT + + // lo = x - (hi + mid) = x - kf * 2^(-5) + double lo = fputil::multiply_add(-0x1.0p-5f, kf, x); + + // hi = floor(kf * 2^(-4)) + // exp2_hi = shift hi to the exponent field of double precision. + int64_t exp2_hi = + static_cast<int64_t>(static_cast<uint64_t>(k >> ExpBase::MID_BITS) + << fputil::FPBits<double>::FRACTION_LEN); + // mh = 2^hi * 2^mid + // mh_bits = bit field of mh + int64_t mh_bits = ExpBase::EXP_2_MID[k & ExpBase::MID_MASK] + exp2_hi; + double mh = fputil::FPBits<double>(static_cast<uint64_t>(mh_bits)).get_val(); + + // Degree-4 polynomial approximating (2^x - 1)/x generated by Sollya with: + // > display = hexadecimal; + // > fpminimax((2^x - 1)/x, 4, [|D...|], [-2^-6, 2^-6]); + constexpr double COEFFS[5] = {0x1.62e42fefa39efp-1, 0x1.ebfbdff8131c4p-3, + 0x1.c6b08d7061695p-5, 0x1.3b2b1bee74b2ap-7, + 0x1.5d88091198529p-10}; + double lo_sq = lo * lo; + double c1 = fputil::multiply_add(lo, COEFFS[0], 1.0); + double c2 = fputil::multiply_add(lo, COEFFS[2], COEFFS[1]); + double c3 = fputil::multiply_add(lo, COEFFS[4], COEFFS[3]); + double exp2_lo = fputil::polyeval(lo_sq, c1, c2, c3); + // 2^x - 1 = 2^(hi + mid + lo) - 1 + // = 2^(hi + mid) * 2^lo - 1 + // ~ mh * (1 + lo * P(lo)) - 1 + // = mh * exp2_lo - 1 + return static_cast<float>(fputil::multiply_add(exp2_lo, mh, -1.0)); +} + +} // namespace LIBC_NAMESPACE diff --git a/libc/src/stdio/fseeko.h b/libc/src/stdio/fseeko.h index 77fb412..3202ed2 100644 --- a/libc/src/stdio/fseeko.h +++ b/libc/src/stdio/fseeko.h @@ -10,6 +10,7 @@ #define LLVM_LIBC_SRC_STDIO_FSEEKO_H #include <stdio.h> +#include <unistd.h> namespace LIBC_NAMESPACE { diff --git a/libc/src/stdio/ftello.h b/libc/src/stdio/ftello.h index 5ab17f9..0fdf13a 100644 --- a/libc/src/stdio/ftello.h +++ b/libc/src/stdio/ftello.h @@ -10,6 +10,7 @@ #define LLVM_LIBC_SRC_STDIO_FTELLO_H #include <stdio.h> +#include <unistd.h> namespace LIBC_NAMESPACE { diff --git a/libc/test/src/__support/integer_literals_test.cpp b/libc/test/src/__support/integer_literals_test.cpp index 5298cf3..cbc906a 100644 --- a/libc/test/src/__support/integer_literals_test.cpp +++ b/libc/test/src/__support/integer_literals_test.cpp @@ -133,3 +133,24 @@ TEST(LlvmLibcIntegerLiteralTest, u256) { U256_MAX, 0xFFFFFFFF'FFFFFFFF'FFFFFFFF'FFFFFFFF'FFFFFFFF'FFFFFFFF'FFFFFFFF'FFFFFFFF_u256); } + +TEST(LlvmLibcIntegerLiteralTest, parse_bigint) { + using T = LIBC_NAMESPACE::Int<128>; + struct { + const char *str; + T expected; + } constexpr TEST_CASES[] = { + {"0", 0}, {"-1", -1}, {"+1", 1}, {"-0xFF", -255}, {"-0b11", -3}, + }; + for (auto tc : TEST_CASES) { + T actual = LIBC_NAMESPACE::parse_bigint<T>(tc.str); + EXPECT_EQ(actual, tc.expected); + } +} + +TEST(LlvmLibcIntegerLiteralTest, parse_bigint_invalid) { + using T = LIBC_NAMESPACE::Int<128>; + const T expected; // default construction + EXPECT_EQ(LIBC_NAMESPACE::parse_bigint<T>(nullptr), expected); + EXPECT_EQ(LIBC_NAMESPACE::parse_bigint<T>(""), expected); +} diff --git a/libc/test/src/__support/math_extras_test.cpp b/libc/test/src/__support/math_extras_test.cpp index e88b3e1..401e631e 100644 --- a/libc/test/src/__support/math_extras_test.cpp +++ b/libc/test/src/__support/math_extras_test.cpp @@ -101,4 +101,61 @@ TYPED_TEST(LlvmLibcBitTest, CountZeros, UnsignedTypesNoBigInt) { EXPECT_EQ(count_zeros<T>(cpp::numeric_limits<T>::max() >> i), i); } +using UnsignedTypes = testing::TypeList< +#if defined(__SIZEOF_INT128__) + __uint128_t, +#endif + unsigned char, unsigned short, unsigned int, unsigned long, + unsigned long long>; + +TYPED_TEST(LlvmLibcBlockMathExtrasTest, add_overflow, UnsignedTypes) { + constexpr T ZERO = cpp::numeric_limits<T>::min(); + constexpr T ONE(1); + constexpr T MAX = cpp::numeric_limits<T>::max(); + constexpr T BEFORE_MAX = MAX - 1; + + const struct { + T lhs; + T rhs; + T sum; + bool carry; + } TESTS[] = { + {ZERO, ONE, ONE, false}, // 0x00 + 0x01 = 0x01 + {BEFORE_MAX, ONE, MAX, false}, // 0xFE + 0x01 = 0xFF + {MAX, ONE, ZERO, true}, // 0xFF + 0x01 = 0x00 (carry) + {MAX, MAX, BEFORE_MAX, true}, // 0xFF + 0xFF = 0xFE (carry) + }; + for (auto tc : TESTS) { + T sum; + bool carry = add_overflow<T>(tc.lhs, tc.rhs, sum); + EXPECT_EQ(sum, tc.sum); + EXPECT_EQ(carry, tc.carry); + } +} + +TYPED_TEST(LlvmLibcBlockMathExtrasTest, sub_overflow, UnsignedTypes) { + constexpr T ZERO = cpp::numeric_limits<T>::min(); + constexpr T ONE(1); + constexpr T MAX = cpp::numeric_limits<T>::max(); + constexpr T BEFORE_MAX = MAX - 1; + + const struct { + T lhs; + T rhs; + T sub; + bool carry; + } TESTS[] = { + {ONE, ZERO, ONE, false}, // 0x01 - 0x00 = 0x01 + {MAX, MAX, ZERO, false}, // 0xFF - 0xFF = 0x00 + {ZERO, ONE, MAX, true}, // 0x00 - 0x01 = 0xFF (carry) + {BEFORE_MAX, MAX, MAX, true}, // 0xFE - 0xFF = 0xFF (carry) + }; + for (auto tc : TESTS) { + T sub; + bool carry = sub_overflow<T>(tc.lhs, tc.rhs, sub); + EXPECT_EQ(sub, tc.sub); + EXPECT_EQ(carry, tc.carry); + } +} + } // namespace LIBC_NAMESPACE diff --git a/libc/test/src/__support/uint_test.cpp b/libc/test/src/__support/uint_test.cpp index 5764324..5696e54 100644 --- a/libc/test/src/__support/uint_test.cpp +++ b/libc/test/src/__support/uint_test.cpp @@ -8,6 +8,7 @@ #include "src/__support/CPP/optional.h" #include "src/__support/UInt.h" +#include "src/__support/integer_literals.h" // parse_unsigned_bigint #include "src/__support/macros/properties/types.h" // LIBC_TYPES_HAS_INT128 #include "include/llvm-libc-macros/math-macros.h" // HUGE_VALF, HUGE_VALF @@ -15,6 +16,195 @@ namespace LIBC_NAMESPACE { +enum Value { ZERO, ONE, TWO, MIN, MAX }; + +template <typename T> auto create(Value value) { + switch (value) { + case ZERO: + return T(0); + case ONE: + return T(1); + case TWO: + return T(2); + case MIN: + return T::min(); + case MAX: + return T::max(); + } +} + +using Types = testing::TypeList< // +#ifdef LIBC_TYPES_HAS_INT64 + BigInt<64, false, uint64_t>, // 64-bits unsigned (1 x uint64_t) + BigInt<64, true, uint64_t>, // 64-bits signed (1 x uint64_t) +#endif +#ifdef LIBC_TYPES_HAS_INT128 + BigInt<128, false, __uint128_t>, // 128-bits unsigned (1 x __uint128_t) + BigInt<128, true, __uint128_t>, // 128-bits signed (1 x __uint128_t) +#endif + BigInt<16, false, uint16_t>, // 16-bits unsigned (1 x uint16_t) + BigInt<16, true, uint16_t>, // 16-bits signed (1 x uint16_t) + BigInt<64, false, uint16_t>, // 64-bits unsigned (4 x uint16_t) + BigInt<64, true, uint16_t> // 64-bits signed (4 x uint16_t) + >; + +#define ASSERT_SAME(A, B) ASSERT_TRUE((A) == (B)) + +TYPED_TEST(LlvmLibcUIntClassTest, Additions, Types) { + ASSERT_SAME(create<T>(ZERO) + create<T>(ZERO), create<T>(ZERO)); + ASSERT_SAME(create<T>(ONE) + create<T>(ZERO), create<T>(ONE)); + ASSERT_SAME(create<T>(ZERO) + create<T>(ONE), create<T>(ONE)); + ASSERT_SAME(create<T>(ONE) + create<T>(ONE), create<T>(TWO)); + // 2's complement addition works for signed and unsigned types. + // - unsigned : 0xff + 0x01 = 0x00 (255 + 1 = 0) + // - signed : 0xef + 0x01 = 0xf0 (127 + 1 = -128) + ASSERT_SAME(create<T>(MAX) + create<T>(ONE), create<T>(MIN)); +} + +TYPED_TEST(LlvmLibcUIntClassTest, Subtraction, Types) { + ASSERT_SAME(create<T>(ZERO) - create<T>(ZERO), create<T>(ZERO)); + ASSERT_SAME(create<T>(ONE) - create<T>(ONE), create<T>(ZERO)); + ASSERT_SAME(create<T>(ONE) - create<T>(ZERO), create<T>(ONE)); + // 2's complement subtraction works for signed and unsigned types. + // - unsigned : 0x00 - 0x01 = 0xff ( 0 - 1 = 255) + // - signed : 0xf0 - 0x01 = 0xef (-128 - 1 = 127) + ASSERT_SAME(create<T>(MIN) - create<T>(ONE), create<T>(MAX)); +} + +TYPED_TEST(LlvmLibcUIntClassTest, Multiplication, Types) { + ASSERT_SAME(create<T>(ZERO) * create<T>(ZERO), create<T>(ZERO)); + ASSERT_SAME(create<T>(ZERO) * create<T>(ONE), create<T>(ZERO)); + ASSERT_SAME(create<T>(ONE) * create<T>(ZERO), create<T>(ZERO)); + ASSERT_SAME(create<T>(ONE) * create<T>(ONE), create<T>(ONE)); + ASSERT_SAME(create<T>(ONE) * create<T>(TWO), create<T>(TWO)); + ASSERT_SAME(create<T>(TWO) * create<T>(ONE), create<T>(TWO)); + // - unsigned : 0xff x 0xff = 0x01 (mod 0xff) + // - signed : 0xef x 0xef = 0x01 (mod 0xff) + ASSERT_SAME(create<T>(MAX) * create<T>(MAX), create<T>(ONE)); +} + +template <typename T> void print(const char *msg, T value) { + testing::tlog << msg; + IntegerToString<T, radix::Hex> buffer(value); + testing::tlog << buffer.view() << "\n"; +} + +TEST(LlvmLibcUIntClassTest, SignedAddSub) { + // Computations performed by https://www.wolframalpha.com/ + using T = BigInt<128, true, uint32_t>; + const T a = parse_bigint<T>("1927508279017230597"); + const T b = parse_bigint<T>("278789278723478925"); + const T s = parse_bigint<T>("2206297557740709522"); + // Addition + ASSERT_SAME(a + b, s); + ASSERT_SAME(b + a, s); // commutative + // Subtraction + ASSERT_SAME(a - s, -b); + ASSERT_SAME(s - a, b); +} + +TEST(LlvmLibcUIntClassTest, SignedMulDiv) { + // Computations performed by https://www.wolframalpha.com/ + using T = BigInt<128, true, uint16_t>; + struct { + const char *a; + const char *b; + const char *mul; + } const test_cases[] = {{"-4", "3", "-12"}, + {"-3", "-3", "9"}, + {"1927508279017230597", "278789278723478925", + "537368642840747885329125014794668225"}}; + for (auto tc : test_cases) { + const T a = parse_bigint<T>(tc.a); + const T b = parse_bigint<T>(tc.b); + const T mul = parse_bigint<T>(tc.mul); + // Multiplication + ASSERT_SAME(a * b, mul); + ASSERT_SAME(b * a, mul); // commutative + ASSERT_SAME(a * -b, -mul); // sign + ASSERT_SAME(-a * b, -mul); // sign + ASSERT_SAME(-a * -b, mul); // sign + // Division + ASSERT_SAME(mul / a, b); + ASSERT_SAME(mul / b, a); + ASSERT_SAME(-mul / a, -b); // sign + ASSERT_SAME(mul / -a, -b); // sign + ASSERT_SAME(-mul / -a, b); // sign + } +} + +TYPED_TEST(LlvmLibcUIntClassTest, Division, Types) { + ASSERT_SAME(create<T>(ZERO) / create<T>(ONE), create<T>(ZERO)); + ASSERT_SAME(create<T>(MAX) / create<T>(ONE), create<T>(MAX)); + ASSERT_SAME(create<T>(MAX) / create<T>(MAX), create<T>(ONE)); + ASSERT_SAME(create<T>(ONE) / create<T>(ONE), create<T>(ONE)); + if constexpr (T::SIGNED) { + // Special case found by fuzzing. + ASSERT_SAME(create<T>(MIN) / create<T>(MIN), create<T>(ONE)); + } + // - unsigned : 0xff / 0x02 = 0x7f + // - signed : 0xef / 0x02 = 0x77 + ASSERT_SAME(create<T>(MAX) / create<T>(TWO), (create<T>(MAX) >> 1)); + + using word_type = typename T::word_type; + const T zero_one_repeated = T::all_ones() / T(0xff); + const word_type pattern = word_type(~0) / word_type(0xff); + for (const word_type part : zero_one_repeated.val) { + if constexpr (T::SIGNED == false) { + EXPECT_EQ(part, pattern); + } + } +} + +TYPED_TEST(LlvmLibcUIntClassTest, is_neg, Types) { + EXPECT_FALSE(create<T>(ZERO).is_neg()); + EXPECT_FALSE(create<T>(ONE).is_neg()); + EXPECT_FALSE(create<T>(TWO).is_neg()); + EXPECT_EQ(create<T>(MIN).is_neg(), T::SIGNED); + EXPECT_FALSE(create<T>(MAX).is_neg()); +} + +TYPED_TEST(LlvmLibcUIntClassTest, Masks, Types) { + if constexpr (!T::SIGNED) { + constexpr size_t BITS = T::BITS; + // mask_trailing_ones + ASSERT_SAME((mask_trailing_ones<T, 0>()), T::zero()); + ASSERT_SAME((mask_trailing_ones<T, 1>()), T::one()); + ASSERT_SAME((mask_trailing_ones<T, BITS - 1>()), T::all_ones() >> 1); + ASSERT_SAME((mask_trailing_ones<T, BITS>()), T::all_ones()); + // mask_leading_ones + ASSERT_SAME((mask_leading_ones<T, 0>()), T::zero()); + ASSERT_SAME((mask_leading_ones<T, 1>()), T::one() << (BITS - 1)); + ASSERT_SAME((mask_leading_ones<T, BITS - 1>()), T::all_ones() - T::one()); + ASSERT_SAME((mask_leading_ones<T, BITS>()), T::all_ones()); + // mask_trailing_zeros + ASSERT_SAME((mask_trailing_zeros<T, 0>()), T::all_ones()); + ASSERT_SAME((mask_trailing_zeros<T, 1>()), T::all_ones() - T::one()); + ASSERT_SAME((mask_trailing_zeros<T, BITS - 1>()), T::one() << (BITS - 1)); + ASSERT_SAME((mask_trailing_zeros<T, BITS>()), T::zero()); + // mask_trailing_zeros + ASSERT_SAME((mask_leading_zeros<T, 0>()), T::all_ones()); + ASSERT_SAME((mask_leading_zeros<T, 1>()), T::all_ones() >> 1); + ASSERT_SAME((mask_leading_zeros<T, BITS - 1>()), T::one()); + ASSERT_SAME((mask_leading_zeros<T, BITS>()), T::zero()); + } +} + +TYPED_TEST(LlvmLibcUIntClassTest, CountBits, Types) { + if constexpr (!T::SIGNED) { + for (size_t i = 0; i <= T::BITS; ++i) { + const auto l_one = T::all_ones() << i; // 0b111...000 + const auto r_one = T::all_ones() >> i; // 0b000...111 + const int zeros = i; + const int ones = T::BITS - zeros; + ASSERT_EQ(cpp::countr_one(r_one), ones); + ASSERT_EQ(cpp::countl_one(l_one), ones); + ASSERT_EQ(cpp::countr_zero(l_one), zeros); + ASSERT_EQ(cpp::countl_zero(r_one), zeros); + } + } +} + using LL_UInt64 = UInt<64>; // We want to test UInt<128> explicitly. So, for // convenience, we use a sugar which does not conflict with the UInt128 type @@ -561,7 +751,7 @@ TEST(LlvmLibcUIntClassTest, FullMulTests) { LL_UInt##Bits a = ~LL_UInt##Bits(0); \ LL_UInt##Bits hi = a.quick_mul_hi(a); \ LL_UInt##Bits trunc = static_cast<LL_UInt##Bits>(a.ful_mul(a) >> Bits); \ - uint64_t overflow = trunc.sub(hi); \ + uint64_t overflow = trunc.sub_overflow(hi); \ EXPECT_EQ(overflow, uint64_t(0)); \ EXPECT_LE(uint64_t(trunc), uint64_t(Error)); \ } while (0) diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index f8f0f8b..bbf8f07 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -638,6 +638,21 @@ add_fp_unittest( ) add_fp_unittest( + exp2m1f_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + exp2m1f_test.cpp + DEPENDS + libc.include.llvm-libc-macros.math_macros + libc.src.errno.errno + libc.src.math.exp2m1f + libc.src.__support.CPP.array + libc.src.__support.FPUtil.fp_bits +) + +add_fp_unittest( exp10f_test NEED_MPFR SUITE diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt index df32dd4..6b2f3dd 100644 --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -143,6 +143,21 @@ add_fp_unittest( ) add_fp_unittest( + exp2m1f_test + NO_RUN_POSTBUILD + NEED_MPFR + SUITE + libc_math_exhaustive_tests + SRCS + exp2m1f_test.cpp + DEPENDS + .exhaustive_test + libc.src.math.exp2m1f + LINK_LIBRARIES + -lpthread +) + +add_fp_unittest( exp10f_test NO_RUN_POSTBUILD NEED_MPFR diff --git a/libc/test/src/math/exhaustive/exp2m1f_test.cpp b/libc/test/src/math/exhaustive/exp2m1f_test.cpp new file mode 100644 index 0000000..2111024 --- /dev/null +++ b/libc/test/src/math/exhaustive/exp2m1f_test.cpp @@ -0,0 +1,33 @@ +//===-- Exhaustive test for exp2m1f ---------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "exhaustive_test.h" +#include "src/math/exp2m1f.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +namespace mpfr = LIBC_NAMESPACE::testing::mpfr; + +using LlvmLibcExp2m1fExhaustiveTest = + LlvmLibcUnaryOpExhaustiveMathTest<float, mpfr::Operation::Exp2m1, + LIBC_NAMESPACE::exp2m1f>; + +// Range: [0, Inf]; +static constexpr uint32_t POS_START = 0x0000'0000U; +static constexpr uint32_t POS_STOP = 0x7f80'0000U; + +TEST_F(LlvmLibcExp2m1fExhaustiveTest, PostiveRange) { + test_full_range_all_roundings(POS_START, POS_STOP); +} + +// Range: [-Inf, 0]; +static constexpr uint32_t NEG_START = 0x8000'0000U; +static constexpr uint32_t NEG_STOP = 0xff80'0000U; + +TEST_F(LlvmLibcExp2m1fExhaustiveTest, NegativeRange) { + test_full_range_all_roundings(NEG_START, NEG_STOP); +} diff --git a/libc/test/src/math/exp2m1f_test.cpp b/libc/test/src/math/exp2m1f_test.cpp new file mode 100644 index 0000000..a0f0da8 --- /dev/null +++ b/libc/test/src/math/exp2m1f_test.cpp @@ -0,0 +1,66 @@ +//===-- Unittests for exp2m1f ---------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "include/llvm-libc-macros/math-macros.h" +#include "src/__support/CPP/array.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/errno/libc_errno.h" +#include "src/math/exp2m1f.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +#include <stdint.h> + +using LlvmLibcExp2m1fTest = LIBC_NAMESPACE::testing::FPTest<float>; + +namespace mpfr = LIBC_NAMESPACE::testing::mpfr; + +TEST_F(LlvmLibcExp2m1fTest, TrickyInputs) { + constexpr LIBC_NAMESPACE::cpp::array<float, 10> INPUTS = { + // EXP2M1F_EXCEPTS_LO + 0x1.36dc8ep-36, + 0x1.224936p-19, + 0x1.d16d2p-20, + 0x1.17949ep-14, + -0x1.9c3e1ep-38, + -0x1.4d89b4p-32, + -0x1.a6eac4p-10, + -0x1.e7526ep-6, + // EXP2M1F_EXCEPTS_HI + 0x1.16a972p-1, + -0x1.9f12acp-5, + }; + + for (float x : INPUTS) { + LIBC_NAMESPACE::libc_errno = 0; + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2m1, x, + LIBC_NAMESPACE::exp2m1f(x), 0.5); + } +} + +TEST_F(LlvmLibcExp2m1fTest, InFloatRange) { + constexpr uint32_t COUNT = 100'000; + constexpr uint32_t STEP = UINT32_MAX / COUNT; + for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) { + float x = FPBits(v).get_val(); + if (isnan(x) || isinf(x)) + continue; + LIBC_NAMESPACE::libc_errno = 0; + float result = LIBC_NAMESPACE::exp2m1f(x); + + // If the computation resulted in an error or did not produce valid result + // in the single-precision floating point range, then ignore comparing with + // MPFR result as MPFR can still produce valid results because of its + // wider precision. + if (isnan(result) || isinf(result) || LIBC_NAMESPACE::libc_errno != 0) + continue; + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2m1, x, + LIBC_NAMESPACE::exp2m1f(x), 0.5); + } +} diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index 5d269dd..4ac1842 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -774,6 +774,17 @@ add_fp_unittest( ) add_fp_unittest( + exp2m1f_test + SUITE + libc-math-smoke-tests + SRCS + exp2m1f_test.cpp + DEPENDS + libc.src.errno.errno + libc.src.math.exp2m1f +) + +add_fp_unittest( exp10f_test SUITE libc-math-smoke-tests diff --git a/libc/test/src/math/smoke/exp2m1f_test.cpp b/libc/test/src/math/smoke/exp2m1f_test.cpp new file mode 100644 index 0000000..2df4353 --- /dev/null +++ b/libc/test/src/math/smoke/exp2m1f_test.cpp @@ -0,0 +1,63 @@ +//===-- Unittests for exp2m1f ---------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "src/errno/libc_errno.h" +#include "src/math/exp2m1f.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" + +using LlvmLibcExp2m1fTest = LIBC_NAMESPACE::testing::FPTest<float>; +using LIBC_NAMESPACE::fputil::testing::ForceRoundingMode; +using LIBC_NAMESPACE::fputil::testing::RoundingMode; + +TEST_F(LlvmLibcExp2m1fTest, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::exp2m1f(aNaN)); + EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::exp2m1f(inf)); + EXPECT_FP_EQ_ALL_ROUNDING(-1.0f, LIBC_NAMESPACE::exp2m1f(neg_inf)); + EXPECT_FP_EQ_ALL_ROUNDING(0.0f, LIBC_NAMESPACE::exp2m1f(0.0f)); + EXPECT_FP_EQ_ALL_ROUNDING(-0.0f, LIBC_NAMESPACE::exp2m1f(-0.0f)); + + EXPECT_FP_EQ_ALL_ROUNDING(1.0f, LIBC_NAMESPACE::exp2m1f(1.0f)); + EXPECT_FP_EQ_ALL_ROUNDING(-0.5f, LIBC_NAMESPACE::exp2m1f(-1.0f)); + EXPECT_FP_EQ_ALL_ROUNDING(3.0f, LIBC_NAMESPACE::exp2m1f(2.0f)); + EXPECT_FP_EQ_ALL_ROUNDING(-0.75f, LIBC_NAMESPACE::exp2m1f(-2.0f)); +} + +TEST_F(LlvmLibcExp2m1fTest, Overflow) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp2m1f(0x1.fffffep+127), + FE_OVERFLOW); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp2m1f(128.0f), + FE_OVERFLOW); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp2m1f(0x1.000002p+7), + FE_OVERFLOW); + EXPECT_MATH_ERRNO(ERANGE); +} + +TEST_F(LlvmLibcExp2m1fTest, Underflow) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ_WITH_EXCEPTION(-1.0f, LIBC_NAMESPACE::exp2m1f(-0x1.fffffep+127), + FE_UNDERFLOW); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ_WITH_EXCEPTION(-1.0f, LIBC_NAMESPACE::exp2m1f(-25.0f), + FE_UNDERFLOW); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ_WITH_EXCEPTION(-1.0f, LIBC_NAMESPACE::exp2m1f(-0x1.900002p4), + FE_UNDERFLOW); + EXPECT_MATH_ERRNO(ERANGE); +} diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index a83f7a7..eaa47da 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -89,7 +89,8 @@ public: // precision. template <typename XType, cpp::enable_if_t<cpp::is_same_v<float, XType>, int> = 0> - explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE, + explicit MPFRNumber(XType x, + unsigned int precision = ExtraPrecision<XType>::VALUE, RoundingMode rounding = RoundingMode::Nearest) : mpfr_precision(precision), mpfr_rounding(get_mpfr_rounding_mode(rounding)) { @@ -99,7 +100,8 @@ public: template <typename XType, cpp::enable_if_t<cpp::is_same_v<double, XType>, int> = 0> - explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE, + explicit MPFRNumber(XType x, + unsigned int precision = ExtraPrecision<XType>::VALUE, RoundingMode rounding = RoundingMode::Nearest) : mpfr_precision(precision), mpfr_rounding(get_mpfr_rounding_mode(rounding)) { @@ -109,7 +111,8 @@ public: template <typename XType, cpp::enable_if_t<cpp::is_same_v<long double, XType>, int> = 0> - explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE, + explicit MPFRNumber(XType x, + unsigned int precision = ExtraPrecision<XType>::VALUE, RoundingMode rounding = RoundingMode::Nearest) : mpfr_precision(precision), mpfr_rounding(get_mpfr_rounding_mode(rounding)) { @@ -119,7 +122,8 @@ public: template <typename XType, cpp::enable_if_t<cpp::is_integral_v<XType>, int> = 0> - explicit MPFRNumber(XType x, int precision = ExtraPrecision<float>::VALUE, + explicit MPFRNumber(XType x, + unsigned int precision = ExtraPrecision<float>::VALUE, RoundingMode rounding = RoundingMode::Nearest) : mpfr_precision(precision), mpfr_rounding(get_mpfr_rounding_mode(rounding)) { @@ -134,6 +138,12 @@ public: mpfr_set(value, other.value, mpfr_rounding); } + MPFRNumber(const MPFRNumber &other, unsigned int precision) + : mpfr_precision(precision), mpfr_rounding(other.mpfr_rounding) { + mpfr_init2(value, mpfr_precision); + mpfr_set(value, other.value, mpfr_rounding); + } + ~MPFRNumber() { mpfr_clear(value); } MPFRNumber &operator=(const MPFRNumber &rhs) { @@ -229,6 +239,36 @@ public: return result; } + MPFRNumber exp2m1() const { + // TODO: Only use mpfr_exp2m1 once CI and buildbots get MPFR >= 4.2.0. +#if MPFR_VERSION_MAJOR > 4 || \ + (MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2) + MPFRNumber result(*this); + mpfr_exp2m1(result.value, value, mpfr_rounding); + return result; +#else + unsigned int prec = mpfr_precision * 3; + MPFRNumber result(*this, prec); + + float f = mpfr_get_flt(abs().value, mpfr_rounding); + if (f > 0.5f && f < 0x1.0p30f) { + mpfr_exp2(result.value, value, mpfr_rounding); + mpfr_sub_ui(result.value, result.value, 1, mpfr_rounding); + return result; + } + + MPFRNumber ln2(2.0f, prec); + // log(2) + mpfr_log(ln2.value, ln2.value, mpfr_rounding); + // x * log(2) + mpfr_mul(result.value, value, ln2.value, mpfr_rounding); + // e^(x * log(2)) - 1 + int ex = mpfr_expm1(result.value, result.value, mpfr_rounding); + mpfr_subnormalize(result.value, ex, mpfr_rounding); + return result; +#endif + } + MPFRNumber exp10() const { MPFRNumber result(*this); mpfr_exp10(result.value, value, mpfr_rounding); @@ -570,6 +610,8 @@ unary_operation(Operation op, InputType input, unsigned int precision, return mpfrInput.exp(); case Operation::Exp2: return mpfrInput.exp2(); + case Operation::Exp2m1: + return mpfrInput.exp2m1(); case Operation::Exp10: return mpfrInput.exp10(); case Operation::Expm1: diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h index d5ff590..0a41ac6 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -37,6 +37,7 @@ enum class Operation : int { Erf, Exp, Exp2, + Exp2m1, Exp10, Expm1, Floor, |