diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt index 1a03683d72e61..b783dd4b04c74 100644 --- a/libc/config/linux/riscv/entrypoints.txt +++ b/libc/config/linux/riscv/entrypoints.txt @@ -986,6 +986,7 @@ if(LIBC_COMPILER_HAS_FIXED_POINT) libc.src.stdfix.idivulr libc.src.stdfix.idivuk libc.src.stdfix.idivulk + libc.src.stdfix.rdivi ) endif() diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 5590f1a15ac57..f6beccc7229dd 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -1019,6 +1019,7 @@ if(LIBC_COMPILER_HAS_FIXED_POINT) libc.src.stdfix.idivulr libc.src.stdfix.idivuk libc.src.stdfix.idivulk + libc.src.stdfix.rdivi ) endif() diff --git a/libc/docs/headers/stdfix.rst b/libc/docs/headers/stdfix.rst index a2f5e94774798..e75b9d8aac440 100644 --- a/libc/docs/headers/stdfix.rst +++ b/libc/docs/headers/stdfix.rst @@ -81,7 +81,7 @@ The following functions are included in the ISO/IEC TR 18037:2008 standard. +---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+ | muli | | | | | | | | | | | | | +---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+ -| \*divi | | | | | | | | | | | | | +| \*divi | | | | |check| | | | | | | | | | +---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+ | round | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | +---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+ diff --git a/libc/include/stdfix.yaml b/libc/include/stdfix.yaml index 5b385124eb63d..451330c3478d2 100644 --- a/libc/include/stdfix.yaml +++ b/libc/include/stdfix.yaml @@ -544,3 +544,11 @@ functions: arguments: - type: unsigned long accum guard: LIBC_COMPILER_HAS_FIXED_POINT + - name: rdivi + standards: + - stdc_ext + return_type: fract + arguments: + - type: int + - type: int + guard: LIBC_COMPILER_HAS_FIXED_POINT diff --git a/libc/src/__support/fixed_point/fx_bits.h b/libc/src/__support/fixed_point/fx_bits.h index 00c6119b4f353..25221d9a33161 100644 --- a/libc/src/__support/fixed_point/fx_bits.h +++ b/libc/src/__support/fixed_point/fx_bits.h @@ -10,9 +10,11 @@ #define LLVM_LIBC_SRC___SUPPORT_FIXED_POINT_FX_BITS_H #include "include/llvm-libc-macros/stdfix-macros.h" +#include "src/__support/CPP/algorithm.h" #include "src/__support/CPP/bit.h" #include "src/__support/CPP/limits.h" // numeric_limits #include "src/__support/CPP/type_traits.h" +#include "src/__support/libc_assert.h" #include "src/__support/macros/attributes.h" // LIBC_INLINE #include "src/__support/macros/config.h" // LIBC_NAMESPACE_DECL #include "src/__support/macros/null_check.h" // LIBC_CRASH_ON_VALUE @@ -21,6 +23,8 @@ #include "fx_rep.h" +#include + #ifdef LIBC_COMPILER_HAS_FIXED_POINT namespace LIBC_NAMESPACE_DECL { @@ -224,6 +228,113 @@ idiv(T x, T y) { return static_cast(result); } +LIBC_INLINE long accum nrstep(long accum d, long accum x0) { + auto v = x0 * (2.lk - (d * x0)); + return v; +} + +// Divide the two integers and return a fixed_point value +// +// For reference, see: +// https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division +// https://stackoverflow.com/a/9231996 + +template LIBC_INLINE constexpr XType divi(int n, int d) { + // If the value of the second operand of the / operator is zero, the + // behavior is undefined. Ref: ISO/IEC TR 18037:2008(E) p.g. 16 + LIBC_CRASH_ON_VALUE(d, 0); + + if (LIBC_UNLIKELY(n == 0)) { + return FXRep::ZERO(); + } + auto is_power_of_two = [](int n) { return (n > 0) && ((n & (n - 1)) == 0); }; + long accum max_val = static_cast(FXRep::MAX()); + long accum min_val = static_cast(FXRep::MIN()); + + if (is_power_of_two(cpp::abs(d))) { + int k = cpp::countr_zero(static_cast(cpp::abs(d))); + constexpr int F = FXRep::FRACTION_LEN; + int64_t scaled_n = static_cast(n) << F; + int64_t res64 = scaled_n >> k; + constexpr int TOTAL_BITS = sizeof(XType) * 8; + const int64_t max_limit = (1LL << (TOTAL_BITS - 1)) - 1; + const int64_t min_limit = -(1LL << (TOTAL_BITS - 1)); + if (res64 > max_limit) { + return FXRep::MAX(); + } else if (res64 < min_limit) { + return FXRep::MIN(); + } + long accum res_accum = + static_cast(res64) / static_cast(1 << F); + res_accum = (d < 0) ? static_cast(-1) * res_accum : res_accum; + if (res_accum > max_val) { + return FXRep::MAX(); + } else if (res_accum < min_val) { + return FXRep::MIN(); + } + return static_cast(res_accum); + } + + bool result_is_negative = ((n < 0) != (d < 0)); + int64_t n64 = static_cast(n); + int64_t d64 = static_cast(d); + + uint64_t nv = static_cast(n64 < 0 ? -n64 : n64); + uint64_t dv = static_cast(d64 < 0 ? -d64 : d64); + + if (d == INT_MIN) { + nv <<= 1; + dv >>= 1; + } + + uint32_t clz = cpp::countl_zero(static_cast(dv)) - 1; + uint64_t scaled_val = dv << clz; + // Scale denominator to be in the range of [0.5,1] + FXBits d_scaled{scaled_val}; + uint64_t scaled_val_n = nv << clz; + // Scale the numerator as much as the denominator to maintain correctness of + // the original equation + FXBits n_scaled{scaled_val_n}; + long accum n_scaled_val = n_scaled.get_val(); + long accum d_scaled_val = d_scaled.get_val(); + // x0 = (48/17) - (32/17) * d_n + long accum a = 0x2.d89d89d8p0lk; // 48/17 = 2.8235294... + long accum b = 0x1.e1e1e1e1p0lk; // 32/17 = 1.8823529... + // Error of the initial approximation, as derived + // from the wikipedia article is + // E0 = 1/17 = 0.059 (5.9%) + long accum initial_approx = a - (b * d_scaled_val); + // Since, 0.5 <= d_scaled_val <= 1.0, 0.9412 <= initial_approx <= 1.88235 + LIBC_ASSERT((initial_approx >= 0x0.78793dd9p0lk) && + (initial_approx <= 0x1.f0f0d845p0lk)); + // Each newton-raphson iteration will square the error, due + // to quadratic convergence. So, + // E1 = (0.059)^2 = 0.0034 + long accum val = nrstep(d_scaled_val, initial_approx); + if constexpr (FXRep::FRACTION_LEN > 8) { + // E2 = 0.0000121 + val = nrstep(d_scaled_val, val); + if constexpr (FXRep::FRACTION_LEN > 16) { + // E3 = 1.468e−10 + val = nrstep(d_scaled_val, val); + } + } + long accum res = n_scaled_val * val; + + if (result_is_negative) { + res *= static_cast(-1); + } + + // Per clause 7.18a.6.1, saturate values on overflow + if (res > max_val) { + return FXRep::MAX(); + } else if (res < min_val) { + return FXRep::MIN(); + } else { + return static_cast(res); + } +} + } // namespace fixed_point } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/stdfix/CMakeLists.txt b/libc/src/stdfix/CMakeLists.txt index 843111e3f80b1..3cbabd17ad34c 100644 --- a/libc/src/stdfix/CMakeLists.txt +++ b/libc/src/stdfix/CMakeLists.txt @@ -89,6 +89,20 @@ foreach(suffix IN ITEMS r lr k lk ur ulr uk ulk) ) endforeach() +foreach(suffix IN ITEMS r) + add_entrypoint_object( + ${suffix}divi + HDRS + ${suffix}divi.h + SRCS + ${suffix}divi.cpp + COMPILE_OPTIONS + ${libc_opt_high_flag} + DEPENDS + libc.src.__support.fixed_point.fx_bits + ) +endforeach() + add_entrypoint_object( uhksqrtus HDRS diff --git a/libc/src/stdfix/rdivi.cpp b/libc/src/stdfix/rdivi.cpp new file mode 100644 index 0000000000000..7021a8b1c8b2a --- /dev/null +++ b/libc/src/stdfix/rdivi.cpp @@ -0,0 +1,21 @@ +//===-- Implementation of rdivi 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 "rdivi.h" +#include "include/llvm-libc-macros/stdfix-macros.h" // fract +#include "src/__support/common.h" // LLVM_LIBC_FUNCTION +#include "src/__support/fixed_point/fx_bits.h" // fixed_point +#include "src/__support/macros/config.h" // LIBC_NAMESPACE_DECL + +namespace LIBC_NAMESPACE_DECL { + +LLVM_LIBC_FUNCTION(fract, rdivi, (int a, int b)) { + return fixed_point::divi(a, b); +} + +} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/stdfix/rdivi.h b/libc/src/stdfix/rdivi.h new file mode 100644 index 0000000000000..aeda1ee9d40f0 --- /dev/null +++ b/libc/src/stdfix/rdivi.h @@ -0,0 +1,21 @@ +//===-- Implementation header for rdivi ------------------------*- 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_STDFIX_RDIVI_H +#define LLVM_LIBC_SRC_STDFIX_RDIVI_H + +#include "include/llvm-libc-macros/stdfix-macros.h" +#include "src/__support/macros/config.h" + +namespace LIBC_NAMESPACE_DECL { + +fract rdivi(int a, int b); + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC_STDFIX_RDIVI_H diff --git a/libc/test/src/stdfix/CMakeLists.txt b/libc/test/src/stdfix/CMakeLists.txt index e2b4bc1805f7c..741522276feaa 100644 --- a/libc/test/src/stdfix/CMakeLists.txt +++ b/libc/test/src/stdfix/CMakeLists.txt @@ -120,6 +120,22 @@ foreach(suffix IN ITEMS r lr k lk ur ulr uk ulk) ) endforeach() +foreach(suffix IN ITEMS r) + add_libc_test( + ${suffix}divi_test + SUITE + libc-stdfix-tests + HDRS + DivITest.h + SRCS + ${suffix}divi_test.cpp + DEPENDS + libc.src.stdfix.${suffix}divi + libc.src.__support.fixed_point.fx_bits + libc.hdr.signal_macros + ) +endforeach() + add_libc_test( uhksqrtus_test SUITE diff --git a/libc/test/src/stdfix/DivITest.h b/libc/test/src/stdfix/DivITest.h new file mode 100644 index 0000000000000..ad15637a3455f --- /dev/null +++ b/libc/test/src/stdfix/DivITest.h @@ -0,0 +1,74 @@ +//===-- Utility class to test fxdivi functions ------------------*- 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 +// +//===----------------------------------------------------------------------===// + +#include "src/__support/CPP/type_traits.h" +#include "src/__support/fixed_point/fx_bits.h" +#include "src/__support/fixed_point/fx_rep.h" +#include "test/UnitTest/Test.h" + +template XType get_epsilon() = delete; +template <> fract get_epsilon() { return FRACT_EPSILON; } +template <> unsigned fract get_epsilon() { return UFRACT_EPSILON; } +template <> long fract get_epsilon() { return LFRACT_EPSILON; } + +template +class DivITest : public LIBC_NAMESPACE::testing::Test { + using FXRep = LIBC_NAMESPACE::fixed_point::FXRep; + using FXBits = LIBC_NAMESPACE::fixed_point::FXBits; + +public: + typedef XType (*DivIFunc)(int, int); + + void testBasic(DivIFunc func) { + XType epsilon = get_epsilon(); + EXPECT_LT((func(2, 3) - 0.666656494140625r), epsilon); + EXPECT_LT((func(3, 4) - 0.75r), epsilon); + EXPECT_LT((func(1043, 2764) - 0.3773516643r), epsilon); + EXPECT_LT((func(60000, 720293) - 0.08329943509r), epsilon); + + EXPECT_EQ(func(128, 256), 0.5r); + EXPECT_EQ(func(1, 2), 0.5r); + EXPECT_EQ(func(1, 4), 0.25r); + EXPECT_EQ(func(1, 8), 0.125r); + EXPECT_EQ(func(1, 16), 0.0625r); + + EXPECT_EQ(func(-1, 2), -0.5r); + EXPECT_EQ(func(1, -4), -0.25r); + EXPECT_EQ(func(-1, 8), -0.125r); + EXPECT_EQ(func(1, -16), -0.0625r); + } + + void testSpecial(DivIFunc func) { + XType epsilon = get_epsilon(); + EXPECT_EQ(func(0, 10), 0.r); + EXPECT_EQ(func(0, -10), 0.r); + EXPECT_EQ(func(-(1 << FRACT_FBIT), 1 << FRACT_FBIT), FRACT_MIN); + EXPECT_EQ(func((1 << FRACT_FBIT) - 1, 1 << FRACT_FBIT), FRACT_MAX); + // From Section 7.18a.6.1, functions returning a fixed-point value, the + // return value is saturated on overflow. + EXPECT_EQ(func(INT_MAX, INT_MAX), FRACT_MAX); + EXPECT_LT(func(INT_MAX - 1, INT_MAX) - 0.99999999r, epsilon); + EXPECT_EQ(func(INT_MIN, INT_MAX), FRACT_MIN); + // Expecting 0 here as fract is not precise enough to + // handle 1/INT_MAX + EXPECT_LT(func(1, INT_MAX) - 0.r, epsilon); + // This results in 1.1739, which should be saturated to FRACT_MAX + EXPECT_EQ(func(27, 23), FRACT_MAX); + + EXPECT_EQ(func(INT_MIN, 1), FRACT_MIN); + EXPECT_LT(func(1, INT_MIN) - 0.r, epsilon); + + EXPECT_EQ(func(INT_MIN, INT_MIN), 1.r); + } +}; + +#define LIST_DIVI_TESTS(Name, XType, func) \ + using LlvmLibc##Name##diviTest = DivITest; \ + TEST_F(LlvmLibc##Name##diviTest, Basic) { testBasic(&func); } \ + TEST_F(LlvmLibc##Name##diviTest, Special) { testSpecial(&func); } \ + static_assert(true, "Require semicolon.") diff --git a/libc/test/src/stdfix/rdivi_test.cpp b/libc/test/src/stdfix/rdivi_test.cpp new file mode 100644 index 0000000000000..10ab366679a36 --- /dev/null +++ b/libc/test/src/stdfix/rdivi_test.cpp @@ -0,0 +1,14 @@ +//===-- Unittests for rdivi -----------------------------------------------===// +// +// 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 "DivITest.h" + +#include "llvm-libc-macros/stdfix-macros.h" // fract +#include "src/stdfix/rdivi.h" + +LIST_DIVI_TESTS(r, fract, LIBC_NAMESPACE::rdivi);