diff options
Diffstat (limited to 'libc/src')
36 files changed, 450 insertions, 365 deletions
diff --git a/libc/src/__support/GPU/allocator.cpp b/libc/src/__support/GPU/allocator.cpp index 866aea7..2b78c4d 100644 --- a/libc/src/__support/GPU/allocator.cpp +++ b/libc/src/__support/GPU/allocator.cpp @@ -16,6 +16,7 @@ #include "allocator.h" +#include "src/__support/CPP/algorithm.h" #include "src/__support/CPP/atomic.h" #include "src/__support/CPP/bit.h" #include "src/__support/CPP/new.h" @@ -31,14 +32,12 @@ constexpr static uint64_t SLAB_SIZE = /* 2 MiB */ 2ull * 1024 * 1024; constexpr static uint64_t ARRAY_SIZE = MAX_SIZE / SLAB_SIZE; constexpr static uint64_t SLAB_ALIGNMENT = SLAB_SIZE - 1; constexpr static uint32_t BITS_IN_WORD = sizeof(uint32_t) * 8; +constexpr static uint32_t BITS_IN_DWORD = sizeof(uint64_t) * 8; constexpr static uint32_t MIN_SIZE = 16; constexpr static uint32_t MIN_ALIGNMENT = MIN_SIZE - 1; // The number of times to attempt claiming an in-progress slab allocation. -constexpr static uint32_t MAX_TRIES = 128; - -// A sentinel used to indicate an invalid but non-null pointer value. -constexpr static uint64_t SENTINEL = cpp::numeric_limits<uint64_t>::max(); +constexpr static uint32_t MAX_TRIES = 1024; static_assert(!(ARRAY_SIZE & (ARRAY_SIZE - 1)), "Must be a power of two"); @@ -70,8 +69,8 @@ static void rpc_free(void *ptr) { // Convert a potentially disjoint bitmask into an increasing integer per-lane // for use with indexing between gpu lanes. -static inline uint32_t lane_count(uint64_t lane_mask) { - return cpp::popcount(lane_mask & ((uint64_t(1) << gpu::get_lane_id()) - 1)); +static inline uint32_t lane_count(uint64_t lane_mask, uint32_t id) { + return cpp::popcount(lane_mask & ((uint64_t(1) << id) - 1)); } // Obtain an initial value to seed a random number generator. We use the rounded @@ -133,7 +132,8 @@ static inline constexpr T round_up(const T x) { void uniform_memset(uint32_t *s, uint32_t c, uint32_t n, uint64_t uniform) { uint64_t mask = gpu::get_lane_mask(); uint32_t workers = cpp::popcount(uniform); - for (uint32_t i = impl::lane_count(mask & uniform); i < n; i += workers) + for (uint32_t i = impl::lane_count(mask & uniform, gpu::get_lane_id()); i < n; + i += workers) s[i] = c; } @@ -142,10 +142,27 @@ static inline constexpr bool is_pow2(uint64_t x) { return x && (x & (x - 1)) == 0; } -// Where this chunk size should start looking in the global array. -static inline constexpr uint32_t start_index(uint32_t chunk_index) { - return (ARRAY_SIZE * impl::get_chunk_id(chunk_index)) / - impl::get_chunk_id(SLAB_SIZE / 2); +// Where this chunk size should start looking in the global array. Small +// allocations are much more likely than large ones, so we give them the most +// space. We use a cubic easing function normalized on the possible chunks. +static inline constexpr uint32_t get_start_index(uint32_t chunk_size) { + constexpr uint32_t max_chunk = impl::get_chunk_id(SLAB_SIZE / 2); + uint64_t norm = + (1 << 16) - (impl::get_chunk_id(chunk_size) << 16) / max_chunk; + uint64_t bias = (norm * norm * norm) >> 32; + uint64_t inv = (1 << 16) - bias; + return static_cast<uint32_t>(((ARRAY_SIZE - 1) * inv) >> 16); +} + +// Returns the id of the lane below this one that acts as its leader. +static inline uint32_t get_leader_id(uint64_t ballot, uint32_t id) { + uint64_t mask = id < BITS_IN_DWORD ? ~0ull << (id + 1) : 0; + return BITS_IN_DWORD - cpp::countl_zero(ballot & ~mask) - 1; +} + +// We use a sentinal value to indicate a failed or in-progress allocation. +template <typename T> bool is_sentinel(const T &x) { + return x == cpp::numeric_limits<T>::max(); } } // namespace impl @@ -264,28 +281,33 @@ struct Slab { continue; // We try using any known empty bits from the previous attempt first. - uint32_t start = gpu::shuffle(mask, cpp::countr_zero(uniform & mask), - ~after ? (old_index & ~(BITS_IN_WORD - 1)) + - cpp::countr_zero(~after) - : impl::xorshift32(state)); + uint32_t start = gpu::shuffle( + mask, cpp::countr_zero(uniform & mask), + ~after ? (old_index & ~(BITS_IN_WORD - 1)) + cpp::countr_zero(~after) + : __builtin_align_down(impl::xorshift32(state), BITS_IN_WORD)); - uint32_t id = impl::lane_count(uniform & mask); + // Each lane tries to claim one bit in a single contiguous mask. + uint32_t id = impl::lane_count(uniform & mask, gpu::get_lane_id()); uint32_t index = (start + id) % usable_bits(chunk_size); uint32_t slot = index / BITS_IN_WORD; uint32_t bit = index % BITS_IN_WORD; // Get the mask of bits destined for the same slot and coalesce it. - uint64_t match = uniform & gpu::match_any(mask, slot); - uint32_t length = cpp::popcount(match); - uint32_t bitmask = gpu::shuffle( - mask, cpp::countr_zero(match), - static_cast<uint32_t>((uint64_t(1) << length) - 1) << bit); + uint32_t leader = impl::get_leader_id( + uniform & gpu::ballot(mask, !id || index % BITS_IN_WORD == 0), + gpu::get_lane_id()); + uint32_t length = cpp::popcount(uniform & mask) - + impl::lane_count(uniform & mask, leader); + uint32_t bitmask = + static_cast<uint32_t>( + (uint64_t(1) << cpp::min(length, BITS_IN_WORD)) - 1) + << bit; uint32_t before = 0; - if (gpu::get_lane_id() == static_cast<uint32_t>(cpp::countr_zero(match))) + if (gpu::get_lane_id() == leader) before = cpp::AtomicRef(get_bitfield()[slot]) .fetch_or(bitmask, cpp::MemoryOrder::RELAXED); - before = gpu::shuffle(mask, cpp::countr_zero(match), before); + before = gpu::shuffle(mask, leader, before); if (~before & (1 << bit)) result = ptr_from_index(index, chunk_size); else @@ -323,20 +345,20 @@ struct GuardPtr { private: struct RefCounter { // Indicates that the object is in its deallocation phase and thus invalid. - static constexpr uint64_t INVALID = uint64_t(1) << 63; + static constexpr uint32_t INVALID = uint32_t(1) << 31; // If a read preempts an unlock call we indicate this so the following // unlock call can swap out the helped bit and maintain exclusive ownership. - static constexpr uint64_t HELPED = uint64_t(1) << 62; + static constexpr uint32_t HELPED = uint32_t(1) << 30; // Resets the reference counter, cannot be reset to zero safely. - void reset(uint32_t n, uint64_t &count) { + void reset(uint32_t n, uint32_t &count) { counter.store(n, cpp::MemoryOrder::RELAXED); count = n; } // Acquire a slot in the reference counter if it is not invalid. - bool acquire(uint32_t n, uint64_t &count) { + bool acquire(uint32_t n, uint32_t &count) { count = counter.fetch_add(n, cpp::MemoryOrder::RELAXED) + n; return (count & INVALID) == 0; } @@ -349,7 +371,7 @@ private: // another thread resurrected the counter and we quit, or a parallel read // helped us invalidating it. For the latter, claim that flag and return. if (counter.fetch_sub(n, cpp::MemoryOrder::RELAXED) == n) { - uint64_t expected = 0; + uint32_t expected = 0; if (counter.compare_exchange_strong(expected, INVALID, cpp::MemoryOrder::RELAXED, cpp::MemoryOrder::RELAXED)) @@ -372,28 +394,29 @@ private: return (val & INVALID) ? 0 : val; } - cpp::Atomic<uint64_t> counter{0}; + cpp::Atomic<uint32_t> counter{0}; }; - cpp::Atomic<Slab *> ptr{nullptr}; - RefCounter ref{}; + cpp::Atomic<Slab *> ptr; + RefCounter ref; // Should be called be a single lane for each different pointer. template <typename... Args> - Slab *try_lock_impl(uint32_t n, uint64_t &count, Args &&...args) { + Slab *try_lock_impl(uint32_t n, uint32_t &count, Args &&...args) { Slab *expected = ptr.load(cpp::MemoryOrder::RELAXED); if (!expected && ptr.compare_exchange_strong( - expected, reinterpret_cast<Slab *>(SENTINEL), + expected, + reinterpret_cast<Slab *>(cpp::numeric_limits<uintptr_t>::max()), cpp::MemoryOrder::RELAXED, cpp::MemoryOrder::RELAXED)) { - count = cpp::numeric_limits<uint64_t>::max(); + count = cpp::numeric_limits<uint32_t>::max(); void *raw = impl::rpc_allocate(sizeof(Slab)); if (!raw) return nullptr; return new (raw) Slab(cpp::forward<Args>(args)...); } - if (!expected || expected == reinterpret_cast<Slab *>(SENTINEL)) + if (!expected || impl::is_sentinel(reinterpret_cast<uintptr_t>(expected))) return nullptr; if (!ref.acquire(n, count)) @@ -405,7 +428,7 @@ private: // Finalize the associated memory and signal that it is ready to use by // resetting the counter. - void finalize(Slab *mem, uint32_t n, uint64_t &count) { + void finalize(Slab *mem, uint32_t n, uint32_t &count) { cpp::atomic_thread_fence(cpp::MemoryOrder::RELEASE); ptr.store(mem, cpp::MemoryOrder::RELAXED); cpp::atomic_thread_fence(cpp::MemoryOrder::ACQUIRE); @@ -418,7 +441,7 @@ public: // The uniform mask represents which lanes share the same pointer. For each // uniform value we elect a leader to handle it on behalf of the other lanes. template <typename... Args> - Slab *try_lock(uint64_t lane_mask, uint64_t uniform, uint64_t &count, + Slab *try_lock(uint64_t lane_mask, uint64_t uniform, uint32_t &count, Args &&...args) { count = 0; Slab *result = nullptr; @@ -433,14 +456,15 @@ public: // We defer storing the newly allocated slab until now so that we can use // multiple lanes to initialize it and release it for use. - if (count == cpp::numeric_limits<uint64_t>::max()) { + if (impl::is_sentinel(count)) { result->initialize(uniform); if (gpu::get_lane_id() == uint32_t(cpp::countr_zero(uniform))) finalize(result, cpp::popcount(uniform), count); } - if (count != cpp::numeric_limits<uint64_t>::max()) - count = count - cpp::popcount(uniform) + impl::lane_count(uniform) + 1; + if (!impl::is_sentinel(count)) + count = count - cpp::popcount(uniform) + + impl::lane_count(uniform, gpu::get_lane_id()) + 1; return result; } @@ -469,7 +493,7 @@ static GuardPtr slots[ARRAY_SIZE] = {}; // Keep a cache of the last successful slot for each chunk size. Initialize it // to an even spread of the total size. Must be updated if the chunking scheme // changes. -#define S(X) (impl::start_index(X)) +#define S(X) (impl::get_start_index(X)) static cpp::Atomic<uint32_t> indices[] = { S(16), S(32), S(48), S(64), S(96), S(112), S(128), S(192), S(224), S(256), S(384), S(448), S(512), S(768), @@ -481,26 +505,28 @@ static cpp::Atomic<uint32_t> indices[] = { #undef S // Tries to find a slab in the table that can support the given chunk size. -static Slab *find_slab(uint32_t chunk_size) { +static Slab *find_slab(uint32_t chunk_size, uint64_t &uniform) { // We start at the index of the last successful allocation for this kind. uint32_t chunk_id = impl::get_chunk_id(chunk_size); uint32_t start = indices[chunk_id].load(cpp::MemoryOrder::RELAXED); - uint64_t uniform = gpu::match_any(gpu::get_lane_mask(), chunk_size); - for (uint32_t offset = 0; offset < ARRAY_SIZE; ++offset) { + for (uint32_t offset = 0; offset <= ARRAY_SIZE; ++offset) { uint32_t index = - !offset ? start : (impl::start_index(chunk_size) + offset) % ARRAY_SIZE; + !offset ? start + : (impl::get_start_index(chunk_size) + offset - 1) % ARRAY_SIZE; - if (slots[index].use_count() < Slab::available_chunks(chunk_size)) { + if (!offset || + slots[index].use_count() < Slab::available_chunks(chunk_size)) { uint64_t lane_mask = gpu::get_lane_mask(); - uint64_t reserved = 0; + uint32_t reserved = 0; Slab *slab = slots[index].try_lock(lane_mask, uniform & lane_mask, reserved, chunk_size, index); // If there is a slab allocation in progress we retry a few times. for (uint32_t retries = 0; - retries < MAX_TRIES && !slab && reserved != SENTINEL; retries++) { + !slab && !impl::is_sentinel(reserved) && retries < MAX_TRIES; + retries++) { uint64_t lane_mask = gpu::get_lane_mask(); slab = slots[index].try_lock(lane_mask, uniform & lane_mask, reserved, chunk_size, index); @@ -514,13 +540,17 @@ static Slab *find_slab(uint32_t chunk_size) { slab->get_chunk_size() == chunk_size) { if (index != start) indices[chunk_id].store(index, cpp::MemoryOrder::RELAXED); + uniform = uniform & gpu::get_lane_mask(); return slab; } else if (slab && (reserved > Slab::available_chunks(chunk_size) || slab->get_chunk_size() != chunk_size)) { slots[index].unlock(gpu::get_lane_mask(), gpu::get_lane_mask() & uniform); - } else if (!slab && reserved == SENTINEL) { + } else if (!slab && impl::is_sentinel(reserved)) { + uniform = uniform & gpu::get_lane_mask(); return nullptr; + } else { + sleep_briefly(); } } } @@ -547,12 +577,12 @@ void *allocate(uint64_t size) { // Try to find a slab for the rounded up chunk size and allocate from it. uint32_t chunk_size = impl::get_chunk_size(static_cast<uint32_t>(size)); - Slab *slab = find_slab(chunk_size); - if (!slab || slab == reinterpret_cast<Slab *>(SENTINEL)) + uint64_t uniform = gpu::match_any(gpu::get_lane_mask(), chunk_size); + Slab *slab = find_slab(chunk_size, uniform); + if (!slab || impl::is_sentinel(reinterpret_cast<uintptr_t>(slab))) return nullptr; uint64_t lane_mask = gpu::get_lane_mask(); - uint64_t uniform = gpu::match_any(lane_mask, slab->get_global_index()); void *ptr = slab->allocate(lane_mask, uniform); return ptr; } diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index b096c61..13f46a1 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -141,6 +141,20 @@ add_header_library( ) add_header_library( + asinhf + HDRS + asinhf.h + DEPENDS + .acoshf_utils + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.sqrt + libc.src.__support.macros.config + libc.src.__support.macros.optimization +) + +add_header_library( asinf HDRS asinf.h @@ -157,6 +171,20 @@ add_header_library( ) add_header_library( + asinf16 + HDRS + asinf16.h + DEPENDS + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.cast + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.sqrt + libc.src.__support.macros.optimization +) + +add_header_library( erff HDRS erff.h diff --git a/libc/src/__support/math/acos.h b/libc/src/__support/math/acos.h index a52ead7..0e1e413 100644 --- a/libc/src/__support/math/acos.h +++ b/libc/src/__support/math/acos.h @@ -24,7 +24,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr double acos(double x) { +LIBC_INLINE static constexpr double acos(double x) { using DoubleDouble = fputil::DoubleDouble; using namespace asin_internal; using FPBits = fputil::FPBits<double>; diff --git a/libc/src/__support/math/acosf.h b/libc/src/__support/math/acosf.h index 153087e..7a0c0e5 100644 --- a/libc/src/__support/math/acosf.h +++ b/libc/src/__support/math/acosf.h @@ -45,7 +45,7 @@ static constexpr fputil::ExceptValues<float, N_EXCEPTS> ACOSF_EXCEPTS = {{ } // namespace acosf_internal -static constexpr float acosf(float x) { +LIBC_INLINE static constexpr float acosf(float x) { using namespace acosf_internal; using namespace inv_trigf_utils_internal; using FPBits = typename fputil::FPBits<float>; diff --git a/libc/src/__support/math/acosf16.h b/libc/src/__support/math/acosf16.h index 58d3761..3f0e002 100644 --- a/libc/src/__support/math/acosf16.h +++ b/libc/src/__support/math/acosf16.h @@ -26,7 +26,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float16 acosf16(float16 x) { +LIBC_INLINE static constexpr float16 acosf16(float16 x) { // Generated by Sollya using the following command: // > round(pi/2, SG, RN); diff --git a/libc/src/__support/math/acoshf.h b/libc/src/__support/math/acoshf.h index f18f169..4e00311 100644 --- a/libc/src/__support/math/acoshf.h +++ b/libc/src/__support/math/acoshf.h @@ -21,7 +21,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float acoshf(float x) { +LIBC_INLINE static constexpr float acoshf(float x) { using namespace acoshf_internal; using FPBits_t = typename fputil::FPBits<float>; FPBits_t xbits(x); diff --git a/libc/src/__support/math/acoshf16.h b/libc/src/__support/math/acoshf16.h index a02b7b6..e5be2a8 100644 --- a/libc/src/__support/math/acoshf16.h +++ b/libc/src/__support/math/acoshf16.h @@ -28,7 +28,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float16 acoshf16(float16 x) { +LIBC_INLINE static constexpr float16 acoshf16(float16 x) { using namespace acoshf_internal; constexpr size_t N_EXCEPTS = 2; diff --git a/libc/src/__support/math/acospif16.h b/libc/src/__support/math/acospif16.h index 5829aed..cf29c76 100644 --- a/libc/src/__support/math/acospif16.h +++ b/libc/src/__support/math/acospif16.h @@ -25,7 +25,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float16 acospif16(float16 x) { +LIBC_INLINE static constexpr float16 acospif16(float16 x) { using FPBits = fputil::FPBits<float16>; FPBits xbits(x); diff --git a/libc/src/__support/math/asin.h b/libc/src/__support/math/asin.h index 84fc1cf..5e06d04 100644 --- a/libc/src/__support/math/asin.h +++ b/libc/src/__support/math/asin.h @@ -25,7 +25,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr double asin(double x) { +LIBC_INLINE static constexpr double asin(double x) { using namespace asin_internal; using FPBits = fputil::FPBits<double>; diff --git a/libc/src/__support/math/asinf16.h b/libc/src/__support/math/asinf16.h new file mode 100644 index 0000000..3d032a4 --- /dev/null +++ b/libc/src/__support/math/asinf16.h @@ -0,0 +1,146 @@ +//===-- Implementation header for asinf16 -----------------------*- 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___SUPPORT_MATH_ASINF16_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINF16_H + +#include "include/llvm-libc-macros/float16-macros.h" + +#ifdef LIBC_TYPES_HAS_FLOAT16 + +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/cast.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/sqrt.h" +#include "src/__support/macros/optimization.h" + +namespace LIBC_NAMESPACE_DECL { + +namespace math { + +LIBC_INLINE static constexpr float16 asinf16(float16 x) { + + // Generated by Sollya using the following command: + // > round(pi/2, D, RN); + constexpr float PI_2 = 0x1.921fb54442d18p0f; + + using FPBits = fputil::FPBits<float16>; + FPBits xbits(x); + + uint16_t x_u = xbits.uintval(); + uint16_t x_abs = x_u & 0x7fff; + float xf = x; + + // |x| > 0x1p0, |x| > 1, or x is NaN. + if (LIBC_UNLIKELY(x_abs > 0x3c00)) { + // asinf16(NaN) = NaN + if (xbits.is_nan()) { + if (xbits.is_signaling_nan()) { + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + + return x; + } + + // 1 < |x| <= +/-inf + fputil::raise_except_if_required(FE_INVALID); + fputil::set_errno_if_required(EDOM); + + return FPBits::quiet_nan().get_val(); + } + + float xsq = xf * xf; + + // |x| <= 0x1p-1, |x| <= 0.5 + if (x_abs <= 0x3800) { + // asinf16(+/-0) = +/-0 + if (LIBC_UNLIKELY(x_abs == 0)) + return x; + + // Exhaustive tests show that, + // for |x| <= 0x1.878p-9, when: + // x > 0, and rounding upward, or + // x < 0, and rounding downward, then, + // asin(x) = x * 2^-11 + x + // else, in other rounding modes, + // asin(x) = x + if (LIBC_UNLIKELY(x_abs <= 0x1a1e)) { + int rounding = fputil::quick_get_round(); + + if ((xbits.is_pos() && rounding == FE_UPWARD) || + (xbits.is_neg() && rounding == FE_DOWNWARD)) + return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf)); + return x; + } + + // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with: + // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]); + float result = + fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f, + 0x1.43b2d6p-5f, 0x1.a0d73ep-5f); + return fputil::cast<float16>(xf * result); + } + + // When |x| > 0.5, assume that 0.5 < |x| <= 1, + // + // Step-by-step range-reduction proof: + // 1: Let y = asin(x), such that, x = sin(y) + // 2: From complimentary angle identity: + // x = sin(y) = cos(pi/2 - y) + // 3: Let z = pi/2 - y, such that x = cos(z) + // 4: From double angle formula; cos(2A) = 1 - sin^2(A): + // z = 2A, z/2 = A + // cos(z) = 1 - 2 * sin^2(z/2) + // 5: Make sin(z/2) subject of the formula: + // sin(z/2) = sqrt((1 - cos(z))/2) + // 6: Recall [3]; x = cos(z). Therefore: + // sin(z/2) = sqrt((1 - x)/2) + // 7: Let u = (1 - x)/2 + // 8: Therefore: + // asin(sqrt(u)) = z/2 + // 2 * asin(sqrt(u)) = z + // 9: Recall [3], z = pi/2 - y. Therefore: + // y = pi/2 - z + // y = pi/2 - 2 * asin(sqrt(u)) + // 10: Recall [1], y = asin(x). Therefore: + // asin(x) = pi/2 - 2 * asin(sqrt(u)) + // + // WHY? + // 11: Recall [7], u = (1 - x)/2 + // 12: Since 0.5 < x <= 1, therefore: + // 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5 + // + // Hence, we can reuse the same [0, 0.5] domain polynomial approximation for + // Step [10] as `sqrt(u)` is in range. + + // 0x1p-1 < |x| <= 0x1p0, 0.5 < |x| <= 1.0 + float xf_abs = (xf < 0 ? -xf : xf); + float sign = (xbits.uintval() >> 15 == 1 ? -1.0 : 1.0); + float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f); + float u_sqrt = fputil::sqrt<float>(u); + + // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with: + // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]); + float asin_sqrt_u = + u_sqrt * fputil::polyeval(u, 0x1.000002p0f, 0x1.554c2ap-3f, + 0x1.3541ccp-4f, 0x1.43b2d6p-5f, 0x1.a0d73ep-5f); + + return fputil::cast<float16>(sign * + fputil::multiply_add(-2.0f, asin_sqrt_u, PI_2)); +} + +} // namespace math + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LIBC_TYPES_HAS_FLOAT16 + +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINF16_H diff --git a/libc/src/__support/math/asinhf.h b/libc/src/__support/math/asinhf.h new file mode 100644 index 0000000..1c08a6e --- /dev/null +++ b/libc/src/__support/math/asinhf.h @@ -0,0 +1,125 @@ +//===-- Implementation header for asinf -------------------------*- 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___SUPPORT_MATH_ASINHF_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINHF_H + +#include "acoshf_utils.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/sqrt.h" +#include "src/__support/macros/config.h" +#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY + +namespace LIBC_NAMESPACE_DECL { + +namespace math { + +LIBC_INLINE static constexpr float asinhf(float x) { + using namespace acoshf_internal; + using FPBits_t = typename fputil::FPBits<float>; + FPBits_t xbits(x); + uint32_t x_u = xbits.uintval(); + uint32_t x_abs = xbits.abs().uintval(); + + // |x| <= 2^-3 + if (LIBC_UNLIKELY(x_abs <= 0x3e80'0000U)) { + // |x| <= 2^-26 + if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) { + return static_cast<float>(LIBC_UNLIKELY(x_abs == 0) + ? x + : (x - 0x1.5555555555555p-3 * x * x * x)); + } + + double x_d = x; + double x_sq = x_d * x_d; + // Generated by Sollya with: + // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16|], [|D...|], + // [0, 2^-2]); + double p = fputil::polyeval( + x_sq, 0.0, -0x1.555555555551ep-3, 0x1.3333333325495p-4, + -0x1.6db6db5a7622bp-5, 0x1.f1c70f82928c6p-6, -0x1.6e893934266b7p-6, + 0x1.1c0b41d3fbe78p-6, -0x1.c0f47810b3c4fp-7, 0x1.2c8602690143dp-7); + return static_cast<float>(fputil::multiply_add(x_d, p, x_d)); + } + + const double SIGN[2] = {1.0, -1.0}; + double x_sign = SIGN[x_u >> 31]; + double x_d = x; + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + // Helper functions to set results for exceptional cases. + auto round_result_slightly_down = [x_sign](float r) -> float { + return fputil::multiply_add(static_cast<float>(x_sign), r, + static_cast<float>(x_sign) * (-0x1.0p-24f)); + }; + auto round_result_slightly_up = [x_sign](float r) -> float { + return fputil::multiply_add(static_cast<float>(x_sign), r, + static_cast<float>(x_sign) * 0x1.0p-24f); + }; + + if (LIBC_UNLIKELY(x_abs >= 0x4bdd'65a5U)) { + if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) { + if (xbits.is_signaling_nan()) { + fputil::raise_except_if_required(FE_INVALID); + return FPBits_t::quiet_nan().get_val(); + } + + return x; + } + + // Exceptional cases when x > 2^24. + switch (x_abs) { + case 0x4bdd65a5: // |x| = 0x1.bacb4ap24f + return round_result_slightly_down(0x1.1e0696p4f); + case 0x4c803f2c: // |x| = 0x1.007e58p26f + return round_result_slightly_down(0x1.2b786cp4f); + case 0x4f8ffb03: // |x| = 0x1.1ff606p32f + return round_result_slightly_up(0x1.6fdd34p4f); + case 0x5c569e88: // |x| = 0x1.ad3d1p57f + return round_result_slightly_up(0x1.45c146p5f); + case 0x5e68984e: // |x| = 0x1.d1309cp61f + return round_result_slightly_up(0x1.5c9442p5f); + case 0x655890d3: // |x| = 0x1.b121a6p75f + return round_result_slightly_down(0x1.a9a3f2p5f); + case 0x65de7ca6: // |x| = 0x1.bcf94cp76f + return round_result_slightly_up(0x1.af66cp5f); + case 0x6eb1a8ec: // |x| = 0x1.6351d8p94f + return round_result_slightly_down(0x1.08b512p6f); + case 0x7997f30a: // |x| = 0x1.2fe614p116f + return round_result_slightly_up(0x1.451436p6f); + } + } else { + // Exceptional cases when x < 2^24. + if (LIBC_UNLIKELY(x_abs == 0x45abaf26)) { + // |x| = 0x1.575e4cp12f + return round_result_slightly_down(0x1.29becap3f); + } + if (LIBC_UNLIKELY(x_abs == 0x49d29048)) { + // |x| = 0x1.a5209p20f + return round_result_slightly_down(0x1.e1b92p3f); + } + } +#else + if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) + return x; +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + // asinh(x) = log(x + sqrt(x^2 + 1)) + return static_cast<float>( + x_sign * log_eval(fputil::multiply_add( + x_d, x_sign, + fputil::sqrt<double>(fputil::multiply_add(x_d, x_d, 1.0))))); +} + +} // namespace math + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINHF_H diff --git a/libc/src/__support/math/erff.h b/libc/src/__support/math/erff.h index e54ec77..b81be30 100644 --- a/libc/src/__support/math/erff.h +++ b/libc/src/__support/math/erff.h @@ -19,7 +19,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float erff(float x) { +LIBC_INLINE static constexpr float erff(float x) { // Polynomials approximating erf(x)/x on ( k/8, (k + 1)/8 ) generated by // Sollya with: > P = fpminimax(erf(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14|], diff --git a/libc/src/__support/math/exp.h b/libc/src/__support/math/exp.h index 14f0592..83638e8 100644 --- a/libc/src/__support/math/exp.h +++ b/libc/src/__support/math/exp.h @@ -67,7 +67,7 @@ namespace { // Return expm1(dx) / x ~ 1 + dx / 2 + dx^2 / 6 + dx^3 / 24. // For |dx| < 2^-13 + 2^-30: // | output - expm1(dx) / dx | < 2^-51. -static double poly_approx_d(double dx) { +LIBC_INLINE static double poly_approx_d(double dx) { // dx^2 double dx2 = dx * dx; // c0 = 1 + dx / 2 @@ -85,7 +85,7 @@ static double poly_approx_d(double dx) { // Return exp(dx) ~ 1 + dx + dx^2 / 2 + ... + dx^6 / 720 // For |dx| < 2^-13 + 2^-30: // | output - exp(dx) | < 2^-101 -static DoubleDouble poly_approx_dd(const DoubleDouble &dx) { +LIBC_INLINE static DoubleDouble poly_approx_dd(const DoubleDouble &dx) { // Taylor polynomial. constexpr DoubleDouble COEFFS[] = { {0, 0x1p0}, // 1 @@ -106,7 +106,7 @@ static DoubleDouble poly_approx_dd(const DoubleDouble &dx) { // Return exp(dx) ~ 1 + dx + dx^2 / 2 + ... + dx^7 / 5040 // For |dx| < 2^-13 + 2^-30: // | output - exp(dx) | < 2^-126. -static Float128 poly_approx_f128(const Float128 &dx) { +LIBC_INLINE static Float128 poly_approx_f128(const Float128 &dx) { constexpr Float128 COEFFS_128[]{ {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0 {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0 @@ -127,7 +127,7 @@ static Float128 poly_approx_f128(const Float128 &dx) { // Compute exp(x) using 128-bit precision. // TODO(lntue): investigate triple-double precision implementation for this // step. -static Float128 exp_f128(double x, double kd, int idx1, int idx2) { +LIBC_INLINE static Float128 exp_f128(double x, double kd, int idx1, int idx2) { // Recalculate dx: double t1 = fputil::multiply_add(kd, MLOG_2_EXP2_M12_HI, x); // exact @@ -160,8 +160,8 @@ static Float128 exp_f128(double x, double kd, int idx1, int idx2) { } // Compute exp(x) with double-double precision. -static DoubleDouble exp_double_double(double x, double kd, - const DoubleDouble &exp_mid) { +LIBC_INLINE static DoubleDouble exp_double_double(double x, double kd, + const DoubleDouble &exp_mid) { // Recalculate dx: // dx = x - k * 2^-12 * log(2) double t1 = fputil::multiply_add(kd, MLOG_2_EXP2_M12_HI, x); // exact @@ -184,7 +184,7 @@ static DoubleDouble exp_double_double(double x, double kd, // Check for exceptional cases when // |x| <= 2^-53 or x < log(2^-1075) or x >= 0x1.6232bdd7abcd3p+9 -static double set_exceptional(double x) { +LIBC_INLINE static double set_exceptional(double x) { using FPBits = typename fputil::FPBits<double>; FPBits xbits(x); diff --git a/libc/src/__support/math/exp10.h b/libc/src/__support/math/exp10.h index fa60e40c..12a09d7 100644 --- a/libc/src/__support/math/exp10.h +++ b/libc/src/__support/math/exp10.h @@ -83,7 +83,8 @@ LIBC_INLINE static double exp10_poly_approx_d(double dx) { // > P = fpminimax((10^x - 1)/x, 5, [|DD...|], [-2^-14, 2^-14]); // Error bounds: // | output - 10^(dx) | < 2^-101 -static constexpr DoubleDouble exp10_poly_approx_dd(const DoubleDouble &dx) { +LIBC_INLINE static constexpr DoubleDouble +exp10_poly_approx_dd(const DoubleDouble &dx) { // Taylor polynomial. constexpr DoubleDouble COEFFS[] = { {0, 0x1p0}, @@ -105,7 +106,8 @@ static constexpr DoubleDouble exp10_poly_approx_dd(const DoubleDouble &dx) { // Return exp(dx) ~ 1 + a0 * dx + a1 * dx^2 + ... + a6 * dx^7 // For |dx| < 2^-14: // | output - 10^dx | < 1.5 * 2^-124. -static constexpr Float128 exp10_poly_approx_f128(const Float128 &dx) { +LIBC_INLINE static constexpr Float128 +exp10_poly_approx_f128(const Float128 &dx) { constexpr Float128 COEFFS_128[]{ {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0 {Sign::POS, -126, 0x935d8ddd'aaa8ac16'ea56d62b'82d30a2d_u128}, @@ -126,7 +128,8 @@ static constexpr Float128 exp10_poly_approx_f128(const Float128 &dx) { // Compute 10^(x) using 128-bit precision. // TODO(lntue): investigate triple-double precision implementation for this // step. -static Float128 exp10_f128(double x, double kd, int idx1, int idx2) { +LIBC_INLINE static Float128 exp10_f128(double x, double kd, int idx1, + int idx2) { double t1 = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact double t2 = kd * MLOG10_2_EXP2_M12_MID_32; // exact double t3 = kd * MLOG10_2_EXP2_M12_LO; // Error < 2^-144 @@ -157,8 +160,8 @@ static Float128 exp10_f128(double x, double kd, int idx1, int idx2) { } // Compute 10^x with double-double precision. -static DoubleDouble exp10_double_double(double x, double kd, - const DoubleDouble &exp_mid) { +LIBC_INLINE static DoubleDouble +exp10_double_double(double x, double kd, const DoubleDouble &exp_mid) { // Recalculate dx: // dx = x - k * 2^-12 * log10(2) double t1 = fputil::multiply_add(kd, MLOG10_2_EXP2_M12_HI, x); // exact @@ -180,7 +183,7 @@ static DoubleDouble exp10_double_double(double x, double kd, #endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS // When output is denormal. -static double exp10_denorm(double x) { +LIBC_INLINE static double exp10_denorm(double x) { // Range reduction. double tmp = fputil::multiply_add(x, LOG2_10, 0x1.8000'0000'4p21); int k = static_cast<int>(cpp::bit_cast<uint64_t>(tmp) >> 19); @@ -234,7 +237,7 @@ static double exp10_denorm(double x) { // * x >= log10(2^1024) // * x <= log10(2^-1022) // * x is inf or nan -static constexpr double exp10_set_exceptional(double x) { +LIBC_INLINE static constexpr double exp10_set_exceptional(double x) { using FPBits = typename fputil::FPBits<double>; FPBits xbits(x); @@ -285,7 +288,7 @@ static constexpr double exp10_set_exceptional(double x) { namespace math { -static constexpr double exp10(double x) { +LIBC_INLINE static constexpr double exp10(double x) { using FPBits = typename fputil::FPBits<double>; FPBits xbits(x); diff --git a/libc/src/__support/math/exp10f.h b/libc/src/__support/math/exp10f.h index 807b4f0..76ae197 100644 --- a/libc/src/__support/math/exp10f.h +++ b/libc/src/__support/math/exp10f.h @@ -20,7 +20,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float exp10f(float x) { +LIBC_INLINE static constexpr float exp10f(float x) { using FPBits = typename fputil::FPBits<float>; FPBits xbits(x); diff --git a/libc/src/__support/math/exp10f16.h b/libc/src/__support/math/exp10f16.h index 0d8b125..3eca867 100644 --- a/libc/src/__support/math/exp10f16.h +++ b/libc/src/__support/math/exp10f16.h @@ -57,7 +57,7 @@ static constexpr fputil::ExceptValues<float16, N_EXP10F16_EXCEPTS> }}; #endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS -static constexpr float16 exp10f16(float16 x) { +LIBC_INLINE static constexpr float16 exp10f16(float16 x) { using FPBits = fputil::FPBits<float16>; FPBits x_bits(x); diff --git a/libc/src/__support/math/exp10f_utils.h b/libc/src/__support/math/exp10f_utils.h index c30def9..010a2f1 100644 --- a/libc/src/__support/math/exp10f_utils.h +++ b/libc/src/__support/math/exp10f_utils.h @@ -89,7 +89,7 @@ struct Exp10Base : public ExpBase { 0x1.0470591dff149p1, 0x1.2bd7c0a9fbc4dp0, 0x1.1429e74a98f43p-1}; - static double powb_lo(double dx) { + LIBC_INLINE static double powb_lo(double dx) { using fputil::multiply_add; double dx2 = dx * dx; // c0 = 1 + COEFFS[0] * dx diff --git a/libc/src/__support/math/exp_utils.h b/libc/src/__support/math/exp_utils.h index fc9ab10..ef408ed 100644 --- a/libc/src/__support/math/exp_utils.h +++ b/libc/src/__support/math/exp_utils.h @@ -22,8 +22,8 @@ namespace LIBC_NAMESPACE_DECL { // So if we scale x up by 2^1022, we can use // double(1.0 + 2^1022 * x) - 1.0 to test how x is rounded in denormal range. template <bool SKIP_ZIV_TEST = false> -static constexpr cpp::optional<double> ziv_test_denorm(int hi, double mid, - double lo, double err) { +LIBC_INLINE static constexpr cpp::optional<double> +ziv_test_denorm(int hi, double mid, double lo, double err) { using FPBits = typename fputil::FPBits<double>; // Scaling factor = 1/(min normal number) = 2^1022 diff --git a/libc/src/__support/math/expf.h b/libc/src/__support/math/expf.h index 88c1514..f7e11be 100644 --- a/libc/src/__support/math/expf.h +++ b/libc/src/__support/math/expf.h @@ -24,7 +24,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float expf(float x) { +LIBC_INLINE static constexpr float expf(float x) { using FPBits = typename fputil::FPBits<float>; FPBits xbits(x); diff --git a/libc/src/__support/math/expf16.h b/libc/src/__support/math/expf16.h index ded28c7..14302a7 100644 --- a/libc/src/__support/math/expf16.h +++ b/libc/src/__support/math/expf16.h @@ -31,7 +31,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float16 expf16(float16 x) { +LIBC_INLINE static constexpr float16 expf16(float16 x) { #ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS constexpr fputil::ExceptValues<float16, 2> EXPF16_EXCEPTS_LO = {{ // (input, RZ output, RU offset, RD offset, RN offset) diff --git a/libc/src/__support/math/expf16_utils.h b/libc/src/__support/math/expf16_utils.h index 8a2fc94..4204dab7 100644 --- a/libc/src/__support/math/expf16_utils.h +++ b/libc/src/__support/math/expf16_utils.h @@ -47,7 +47,8 @@ struct ExpRangeReduction { float exp_lo; }; -[[maybe_unused]] static ExpRangeReduction exp_range_reduction(float16 x) { +[[maybe_unused]] LIBC_INLINE static ExpRangeReduction +exp_range_reduction(float16 x) { // For -18 < x < 12, to compute exp(x), we perform the following range // reduction: find hi, mid, lo, such that: // x = hi + mid + lo, in which diff --git a/libc/src/__support/math/frexpf.h b/libc/src/__support/math/frexpf.h index 4d2f494..7834a12 100644 --- a/libc/src/__support/math/frexpf.h +++ b/libc/src/__support/math/frexpf.h @@ -17,7 +17,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float frexpf(float x, int *exp) { +LIBC_INLINE static constexpr float frexpf(float x, int *exp) { return fputil::frexp(x, *exp); } diff --git a/libc/src/__support/math/frexpf128.h b/libc/src/__support/math/frexpf128.h index 2fd5bc4..5218b26 100644 --- a/libc/src/__support/math/frexpf128.h +++ b/libc/src/__support/math/frexpf128.h @@ -21,7 +21,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float128 frexpf128(float128 x, int *exp) { +LIBC_INLINE static constexpr float128 frexpf128(float128 x, int *exp) { return fputil::frexp(x, *exp); } diff --git a/libc/src/__support/math/frexpf16.h b/libc/src/__support/math/frexpf16.h index 8deeba0..530b61a 100644 --- a/libc/src/__support/math/frexpf16.h +++ b/libc/src/__support/math/frexpf16.h @@ -21,7 +21,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float16 frexpf16(float16 x, int *exp) { +LIBC_INLINE static constexpr float16 frexpf16(float16 x, int *exp) { return fputil::frexp(x, *exp); } diff --git a/libc/src/__support/math/ldexpf.h b/libc/src/__support/math/ldexpf.h index 3a5ec1d..9ef5d96 100644 --- a/libc/src/__support/math/ldexpf.h +++ b/libc/src/__support/math/ldexpf.h @@ -17,7 +17,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float ldexpf(float x, int exp) { +LIBC_INLINE static constexpr float ldexpf(float x, int exp) { return fputil::ldexp(x, exp); } diff --git a/libc/src/__support/math/ldexpf128.h b/libc/src/__support/math/ldexpf128.h index 3625830..4fba20c 100644 --- a/libc/src/__support/math/ldexpf128.h +++ b/libc/src/__support/math/ldexpf128.h @@ -21,7 +21,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float128 ldexpf128(float128 x, int exp) { +LIBC_INLINE static constexpr float128 ldexpf128(float128 x, int exp) { return fputil::ldexp(x, exp); } diff --git a/libc/src/__support/math/ldexpf16.h b/libc/src/__support/math/ldexpf16.h index fbead87..d978d22 100644 --- a/libc/src/__support/math/ldexpf16.h +++ b/libc/src/__support/math/ldexpf16.h @@ -21,7 +21,7 @@ namespace LIBC_NAMESPACE_DECL { namespace math { -static constexpr float16 ldexpf16(float16 x, int exp) { +LIBC_INLINE static constexpr float16 ldexpf16(float16 x, int exp) { return fputil::ldexp(x, exp); } diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index ecf0967..f91feacb 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3889,12 +3889,7 @@ add_entrypoint_object( HDRS ../asinhf.h DEPENDS - .explogxf - libc.src.__support.FPUtil.fp_bits - libc.src.__support.FPUtil.multiply_add - libc.src.__support.FPUtil.polyeval - libc.src.__support.FPUtil.sqrt - libc.src.__support.macros.optimization + libc.src.__support.math.asinhf ) add_entrypoint_object( @@ -3968,16 +3963,7 @@ add_entrypoint_object( HDRS ../asinf16.h DEPENDS - libc.hdr.errno_macros - libc.hdr.fenv_macros - libc.src.__support.FPUtil.cast - 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.sqrt - libc.src.__support.macros.optimization - libc.src.__support.macros.properties.types + libc.src.__support.math.asinf16 ) add_entrypoint_object( diff --git a/libc/src/math/generic/asinf16.cpp b/libc/src/math/generic/asinf16.cpp index 518c384..af8dbfe 100644 --- a/libc/src/math/generic/asinf16.cpp +++ b/libc/src/math/generic/asinf16.cpp @@ -7,127 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/asinf16.h" -#include "hdr/errno_macros.h" -#include "hdr/fenv_macros.h" -#include "src/__support/FPUtil/FEnvImpl.h" -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/FPUtil/cast.h" -#include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/sqrt.h" -#include "src/__support/macros/optimization.h" +#include "src/__support/math/asinf16.h" namespace LIBC_NAMESPACE_DECL { -// Generated by Sollya using the following command: -// > round(pi/2, D, RN); -static constexpr float PI_2 = 0x1.921fb54442d18p0f; - -LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) { - using FPBits = fputil::FPBits<float16>; - FPBits xbits(x); - - uint16_t x_u = xbits.uintval(); - uint16_t x_abs = x_u & 0x7fff; - float xf = x; - - // |x| > 0x1p0, |x| > 1, or x is NaN. - if (LIBC_UNLIKELY(x_abs > 0x3c00)) { - // asinf16(NaN) = NaN - if (xbits.is_nan()) { - if (xbits.is_signaling_nan()) { - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - - return x; - } - - // 1 < |x| <= +/-inf - fputil::raise_except_if_required(FE_INVALID); - fputil::set_errno_if_required(EDOM); - - return FPBits::quiet_nan().get_val(); - } - - float xsq = xf * xf; - - // |x| <= 0x1p-1, |x| <= 0.5 - if (x_abs <= 0x3800) { - // asinf16(+/-0) = +/-0 - if (LIBC_UNLIKELY(x_abs == 0)) - return x; - - // Exhaustive tests show that, - // for |x| <= 0x1.878p-9, when: - // x > 0, and rounding upward, or - // x < 0, and rounding downward, then, - // asin(x) = x * 2^-11 + x - // else, in other rounding modes, - // asin(x) = x - if (LIBC_UNLIKELY(x_abs <= 0x1a1e)) { - int rounding = fputil::quick_get_round(); - - if ((xbits.is_pos() && rounding == FE_UPWARD) || - (xbits.is_neg() && rounding == FE_DOWNWARD)) - return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf)); - return x; - } - - // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with: - // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]); - float result = - fputil::polyeval(xsq, 0x1.000002p0f, 0x1.554c2ap-3f, 0x1.3541ccp-4f, - 0x1.43b2d6p-5f, 0x1.a0d73ep-5f); - return fputil::cast<float16>(xf * result); - } - - // When |x| > 0.5, assume that 0.5 < |x| <= 1, - // - // Step-by-step range-reduction proof: - // 1: Let y = asin(x), such that, x = sin(y) - // 2: From complimentary angle identity: - // x = sin(y) = cos(pi/2 - y) - // 3: Let z = pi/2 - y, such that x = cos(z) - // 4: From double angle formula; cos(2A) = 1 - sin^2(A): - // z = 2A, z/2 = A - // cos(z) = 1 - 2 * sin^2(z/2) - // 5: Make sin(z/2) subject of the formula: - // sin(z/2) = sqrt((1 - cos(z))/2) - // 6: Recall [3]; x = cos(z). Therefore: - // sin(z/2) = sqrt((1 - x)/2) - // 7: Let u = (1 - x)/2 - // 8: Therefore: - // asin(sqrt(u)) = z/2 - // 2 * asin(sqrt(u)) = z - // 9: Recall [3], z = pi/2 - y. Therefore: - // y = pi/2 - z - // y = pi/2 - 2 * asin(sqrt(u)) - // 10: Recall [1], y = asin(x). Therefore: - // asin(x) = pi/2 - 2 * asin(sqrt(u)) - // - // WHY? - // 11: Recall [7], u = (1 - x)/2 - // 12: Since 0.5 < x <= 1, therefore: - // 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5 - // - // Hence, we can reuse the same [0, 0.5] domain polynomial approximation for - // Step [10] as `sqrt(u)` is in range. - - // 0x1p-1 < |x| <= 0x1p0, 0.5 < |x| <= 1.0 - float xf_abs = (xf < 0 ? -xf : xf); - float sign = (xbits.uintval() >> 15 == 1 ? -1.0 : 1.0); - float u = fputil::multiply_add(-0.5f, xf_abs, 0.5f); - float u_sqrt = fputil::sqrt<float>(u); - - // Degree-6 minimax odd polynomial of asin(x) generated by Sollya with: - // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]); - float asin_sqrt_u = - u_sqrt * fputil::polyeval(u, 0x1.000002p0f, 0x1.554c2ap-3f, - 0x1.3541ccp-4f, 0x1.43b2d6p-5f, 0x1.a0d73ep-5f); - - return fputil::cast<float16>(sign * - fputil::multiply_add(-2.0f, asin_sqrt_u, PI_2)); -} +LLVM_LIBC_FUNCTION(float16, asinf16, (float16 x)) { return math::asinf16(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/asinhf.cpp b/libc/src/math/generic/asinhf.cpp index 3aed3bc..45023c8 100644 --- a/libc/src/math/generic/asinhf.cpp +++ b/libc/src/math/generic/asinhf.cpp @@ -7,112 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/asinhf.h" -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" -#include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/sqrt.h" -#include "src/__support/macros/config.h" -#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY -#include "src/math/generic/common_constants.h" -#include "src/math/generic/explogxf.h" +#include "src/__support/math/asinhf.h" namespace LIBC_NAMESPACE_DECL { -LLVM_LIBC_FUNCTION(float, asinhf, (float x)) { - using namespace acoshf_internal; - using FPBits_t = typename fputil::FPBits<float>; - FPBits_t xbits(x); - uint32_t x_u = xbits.uintval(); - uint32_t x_abs = xbits.abs().uintval(); - - // |x| <= 2^-3 - if (LIBC_UNLIKELY(x_abs <= 0x3e80'0000U)) { - // |x| <= 2^-26 - if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) { - return static_cast<float>(LIBC_UNLIKELY(x_abs == 0) - ? x - : (x - 0x1.5555555555555p-3 * x * x * x)); - } - - double x_d = x; - double x_sq = x_d * x_d; - // Generated by Sollya with: - // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16|], [|D...|], - // [0, 2^-2]); - double p = fputil::polyeval( - x_sq, 0.0, -0x1.555555555551ep-3, 0x1.3333333325495p-4, - -0x1.6db6db5a7622bp-5, 0x1.f1c70f82928c6p-6, -0x1.6e893934266b7p-6, - 0x1.1c0b41d3fbe78p-6, -0x1.c0f47810b3c4fp-7, 0x1.2c8602690143dp-7); - return static_cast<float>(fputil::multiply_add(x_d, p, x_d)); - } - - const double SIGN[2] = {1.0, -1.0}; - double x_sign = SIGN[x_u >> 31]; - double x_d = x; - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - // Helper functions to set results for exceptional cases. - auto round_result_slightly_down = [x_sign](float r) -> float { - return fputil::multiply_add(static_cast<float>(x_sign), r, - static_cast<float>(x_sign) * (-0x1.0p-24f)); - }; - auto round_result_slightly_up = [x_sign](float r) -> float { - return fputil::multiply_add(static_cast<float>(x_sign), r, - static_cast<float>(x_sign) * 0x1.0p-24f); - }; - - if (LIBC_UNLIKELY(x_abs >= 0x4bdd'65a5U)) { - if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) { - if (xbits.is_signaling_nan()) { - fputil::raise_except_if_required(FE_INVALID); - return FPBits_t::quiet_nan().get_val(); - } - - return x; - } - - // Exceptional cases when x > 2^24. - switch (x_abs) { - case 0x4bdd65a5: // |x| = 0x1.bacb4ap24f - return round_result_slightly_down(0x1.1e0696p4f); - case 0x4c803f2c: // |x| = 0x1.007e58p26f - return round_result_slightly_down(0x1.2b786cp4f); - case 0x4f8ffb03: // |x| = 0x1.1ff606p32f - return round_result_slightly_up(0x1.6fdd34p4f); - case 0x5c569e88: // |x| = 0x1.ad3d1p57f - return round_result_slightly_up(0x1.45c146p5f); - case 0x5e68984e: // |x| = 0x1.d1309cp61f - return round_result_slightly_up(0x1.5c9442p5f); - case 0x655890d3: // |x| = 0x1.b121a6p75f - return round_result_slightly_down(0x1.a9a3f2p5f); - case 0x65de7ca6: // |x| = 0x1.bcf94cp76f - return round_result_slightly_up(0x1.af66cp5f); - case 0x6eb1a8ec: // |x| = 0x1.6351d8p94f - return round_result_slightly_down(0x1.08b512p6f); - case 0x7997f30a: // |x| = 0x1.2fe614p116f - return round_result_slightly_up(0x1.451436p6f); - } - } else { - // Exceptional cases when x < 2^24. - if (LIBC_UNLIKELY(x_abs == 0x45abaf26)) { - // |x| = 0x1.575e4cp12f - return round_result_slightly_down(0x1.29becap3f); - } - if (LIBC_UNLIKELY(x_abs == 0x49d29048)) { - // |x| = 0x1.a5209p20f - return round_result_slightly_down(0x1.e1b92p3f); - } - } -#else - if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) - return x; -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - - // asinh(x) = log(x + sqrt(x^2 + 1)) - return static_cast<float>( - x_sign * log_eval(fputil::multiply_add( - x_d, x_sign, - fputil::sqrt<double>(fputil::multiply_add(x_d, x_d, 1.0))))); -} +LLVM_LIBC_FUNCTION(float, asinhf, (float x)) { return math::asinhf(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/stdio/baremetal/CMakeLists.txt b/libc/src/stdio/baremetal/CMakeLists.txt index e879230..548938f 100644 --- a/libc/src/stdio/baremetal/CMakeLists.txt +++ b/libc/src/stdio/baremetal/CMakeLists.txt @@ -72,6 +72,7 @@ add_entrypoint_object( ../scanf.h DEPENDS .scanf_internal + libc.include.inttypes libc.src.stdio.scanf_core.scanf_main libc.src.__support.arg_list libc.src.__support.OSUtil.osutil diff --git a/libc/src/stdio/scanf_core/CMakeLists.txt b/libc/src/stdio/scanf_core/CMakeLists.txt index dee125c..561180c 100644 --- a/libc/src/stdio/scanf_core/CMakeLists.txt +++ b/libc/src/stdio/scanf_core/CMakeLists.txt @@ -35,6 +35,7 @@ add_header_library( core_structs.h DEPENDS .scanf_config + libc.include.inttypes libc.src.__support.CPP.string_view libc.src.__support.CPP.bitset libc.src.__support.FPUtil.fp_bits @@ -97,6 +98,7 @@ add_header_library( DEPENDS .reader .core_structs + libc.include.inttypes libc.src.__support.common libc.src.__support.ctype_utils libc.src.__support.CPP.bitset diff --git a/libc/src/wchar/wchar_utils.h b/libc/src/wchar/wchar_utils.h index e0218c7..55a3cee 100644 --- a/libc/src/wchar/wchar_utils.h +++ b/libc/src/wchar/wchar_utils.h @@ -17,13 +17,10 @@ namespace LIBC_NAMESPACE_DECL { namespace internal { -// returns true if the character exists in the string -LIBC_INLINE static bool wcschr(wchar_t c, const wchar_t *str) { - for (int n = 0; str[n]; ++n) { - if (str[n] == c) - return true; - } - return false; +LIBC_INLINE static const wchar_t *wcschr(const wchar_t *s, wchar_t c) { + for (; *s && *s != c; ++s) + ; + return (*s == c) ? s : nullptr; } // bool should be true for wcscspn for complimentary span @@ -32,7 +29,7 @@ LIBC_INLINE static size_t wcsspn(const wchar_t *s1, const wchar_t *s2, bool not_match_set) { size_t i = 0; for (; s1[i]; ++i) { - bool in_set = wcschr(s1[i], s2); + bool in_set = internal::wcschr(s2, s1[i]); if (in_set == not_match_set) return i; } diff --git a/libc/src/wchar/wcschr.cpp b/libc/src/wchar/wcschr.cpp index defc2ce..8ac4916 100644 --- a/libc/src/wchar/wcschr.cpp +++ b/libc/src/wchar/wcschr.cpp @@ -11,15 +11,14 @@ #include "hdr/types/wchar_t.h" #include "src/__support/common.h" #include "src/__support/macros/config.h" +#include "src/__support/macros/null_check.h" +#include "wchar_utils.h" namespace LIBC_NAMESPACE_DECL { LLVM_LIBC_FUNCTION(const wchar_t *, wcschr, (const wchar_t *s, wchar_t c)) { - for (; *s && *s != c; ++s) - ; - if (*s == c) - return s; - return nullptr; + LIBC_CRASH_ON_NULLPTR(s); + return internal::wcschr(s, c); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/wchar/wcspbrk.cpp b/libc/src/wchar/wcspbrk.cpp index a00ba99..5d86a49 100644 --- a/libc/src/wchar/wcspbrk.cpp +++ b/libc/src/wchar/wcspbrk.cpp @@ -11,17 +11,10 @@ #include "hdr/types/wchar_t.h" #include "src/__support/common.h" #include "src/__support/macros/null_check.h" +#include "wchar_utils.h" namespace LIBC_NAMESPACE_DECL { -bool contains_char(const wchar_t *str, wchar_t target) { - for (; *str != L'\0'; str++) - if (*str == target) - return true; - - return false; -} - LLVM_LIBC_FUNCTION(const wchar_t *, wcspbrk, (const wchar_t *src, const wchar_t *breakset)) { LIBC_CRASH_ON_NULLPTR(src); @@ -29,7 +22,7 @@ LLVM_LIBC_FUNCTION(const wchar_t *, wcspbrk, // currently O(n * m), can be further optimized to O(n + m) with a hash set for (int src_idx = 0; src[src_idx] != 0; src_idx++) - if (contains_char(breakset, src[src_idx])) + if (internal::wcschr(breakset, src[src_idx])) return src + src_idx; return nullptr; diff --git a/libc/src/wchar/wcstok.cpp b/libc/src/wchar/wcstok.cpp index 291efc1..ed4f0aa 100644 --- a/libc/src/wchar/wcstok.cpp +++ b/libc/src/wchar/wcstok.cpp @@ -10,18 +10,12 @@ #include "hdr/types/wchar_t.h" #include "src/__support/common.h" +#include "wchar_utils.h" namespace LIBC_NAMESPACE_DECL { -bool isADelimeter(wchar_t wc, const wchar_t *delimiters) { - for (const wchar_t *delim_ptr = delimiters; *delim_ptr != L'\0'; ++delim_ptr) - if (wc == *delim_ptr) - return true; - return false; -} - LLVM_LIBC_FUNCTION(wchar_t *, wcstok, - (wchar_t *__restrict str, const wchar_t *__restrict delim, + (wchar_t *__restrict str, const wchar_t *__restrict delims, wchar_t **__restrict context)) { if (str == nullptr) { if (*context == nullptr) @@ -30,14 +24,13 @@ LLVM_LIBC_FUNCTION(wchar_t *, wcstok, str = *context; } - wchar_t *tok_start, *tok_end; - for (tok_start = str; *tok_start != L'\0' && isADelimeter(*tok_start, delim); - ++tok_start) - ; + wchar_t *tok_start = str; + while (*tok_start != L'\0' && internal::wcschr(delims, *tok_start)) + ++tok_start; - for (tok_end = tok_start; *tok_end != L'\0' && !isADelimeter(*tok_end, delim); - ++tok_end) - ; + wchar_t *tok_end = tok_start; + while (*tok_end != L'\0' && !internal::wcschr(delims, *tok_end)) + ++tok_end; if (*tok_end != L'\0') { *tok_end = L'\0'; |