diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt --- a/libc/config/darwin/arm/entrypoints.txt +++ b/libc/config/darwin/arm/entrypoints.txt @@ -131,6 +131,8 @@ libc.src.math.fmin libc.src.math.fminf libc.src.math.fminl + libc.src.math.fmod + libc.src.math.fmodf libc.src.math.frexp libc.src.math.frexpf libc.src.math.frexpl diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -150,6 +150,8 @@ libc.src.math.fmin libc.src.math.fminf libc.src.math.fminl + libc.src.math.fmod + libc.src.math.fmodf libc.src.math.frexp libc.src.math.frexpf libc.src.math.frexpl diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -156,6 +156,8 @@ libc.src.math.fmax libc.src.math.fmaxf libc.src.math.fmaxl + libc.src.math.fmod + libc.src.math.fmodf libc.src.math.frexp libc.src.math.frexpf libc.src.math.frexpl diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -133,6 +133,8 @@ libc.src.math.fmax libc.src.math.fmaxf libc.src.math.fmaxl + libc.src.math.fmod + libc.src.math.fmodf libc.src.math.frexp libc.src.math.frexpf libc.src.math.frexpl diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -378,6 +378,10 @@ FunctionSpec<"fma", RetValSpec, [ArgSpec, ArgSpec, ArgSpec]>, FunctionSpec<"fmaf", RetValSpec, [ArgSpec, ArgSpec, ArgSpec]>, + FunctionSpec<"fmod", RetValSpec, [ArgSpec, ArgSpec]>, + + FunctionSpec<"fmodf", RetValSpec, [ArgSpec, ArgSpec]>, + FunctionSpec<"frexp", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"frexpf", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"frexpl", RetValSpec, [ArgSpec, ArgSpec]>, diff --git a/libc/src/__support/FPUtil/CMakeLists.txt b/libc/src/__support/FPUtil/CMakeLists.txt --- a/libc/src/__support/FPUtil/CMakeLists.txt +++ b/libc/src/__support/FPUtil/CMakeLists.txt @@ -13,6 +13,7 @@ NormalFloat.h PlatformDefs.h UInt.h + builtin_wrappers.h DEPENDS libc.include.math libc.include.errno diff --git a/libc/src/__support/FPUtil/FPBits.h b/libc/src/__support/FPUtil/FPBits.h --- a/libc/src/__support/FPUtil/FPBits.h +++ b/libc/src/__support/FPUtil/FPBits.h @@ -13,6 +13,8 @@ #include "src/__support/CPP/Bit.h" #include "src/__support/CPP/TypeTraits.h" +#include "src/__support/FPUtil/builtin_wrappers.h" +#include "src/__support/common.h" #include "FloatProperties.h" #include @@ -57,6 +59,11 @@ UIntType get_mantissa() const { return bits & FloatProp::MANTISSA_MASK; } + // The function return mantissa with implicit bit set for normal values. + constexpr UIntType get_explicit_mantissa() { + return (FloatProp::MANTISSA_MASK + 1) | (FloatProp::MANTISSA_MASK & bits); + } + void set_unbiased_exponent(UIntType expVal) { expVal = (expVal << (FloatProp::MANTISSA_WIDTH)) & FloatProp::EXPONENT_MASK; bits &= ~(FloatProp::EXPONENT_MASK); @@ -69,14 +76,12 @@ } void set_sign(bool signVal) { - bits &= ~(FloatProp::SIGN_MASK); - UIntType sign = UIntType(signVal) << (FloatProp::BIT_WIDTH - 1); - bits |= sign; + bits |= FloatProp::SIGN_MASK; + if (!signVal) + bits -= FloatProp::SIGN_MASK; } - bool get_sign() const { - return ((bits & FloatProp::SIGN_MASK) >> (FloatProp::BIT_WIDTH - 1)); - } + bool get_sign() const { return (bits & FloatProp::SIGN_MASK) != 0; } static_assert(sizeof(T) == sizeof(UIntType), "Data type and integral representation have different sizes."); @@ -92,7 +97,7 @@ static constexpr UIntType MAX_NORMAL = ((UIntType(MAX_EXPONENT) - 1) << MantissaWidth::VALUE) | MAX_SUBNORMAL; - // We don't want accidental type promotions/conversions so we require exact + // We don't want accidental type promotions/conversions, so we require exact // type match. template ::Value, int> = 0> @@ -118,42 +123,67 @@ } bool is_zero() const { - return get_mantissa() == 0 && get_unbiased_exponent() == 0; + // Remove sign bit by shift + return (bits << 1) == 0; } bool is_inf() const { - return get_mantissa() == 0 && get_unbiased_exponent() == MAX_EXPONENT; + // return get_mantissa() == 0 && get_unbiased_exponent() == MAX_EXPONENT; + return (bits & (~(FloatProp::SIGN_MASK))) == FloatProp::EXPONENT_MASK; } bool is_nan() const { - return get_unbiased_exponent() == MAX_EXPONENT && get_mantissa() != 0; + return (bits & (~(FloatProp::SIGN_MASK))) > FloatProp::EXPONENT_MASK; } - bool is_inf_or_nan() const { return get_unbiased_exponent() == MAX_EXPONENT; } - - static FPBits zero() { return FPBits(); } + bool is_inf_or_nan() const { + return (bits & FloatProp::EXPONENT_MASK) == FloatProp::EXPONENT_MASK; + } - static FPBits neg_zero() { - return FPBits(UIntType(1) << (sizeof(UIntType) * 8 - 1)); + static constexpr FPBits zero(bool sign = false) { + return FPBits(sign ? FloatProp::SIGN_MASK : UIntType(0)); } - static FPBits inf() { + static constexpr FPBits neg_zero() { return zero(true); } + + static constexpr FPBits inf() { FPBits bits; bits.set_unbiased_exponent(MAX_EXPONENT); return bits; } - static FPBits neg_inf() { + static constexpr FPBits neg_inf() { FPBits bits = inf(); bits.set_sign(1); return bits; } - static T build_nan(UIntType v) { + static constexpr T build_nan(UIntType v) { FPBits bits = inf(); bits.set_mantissa(v); return T(bits); } + + // The function convert integer number and unbiased exponent to proper Float. + // Result = number * 2^(ep+1) + // Be careful! The function add to +1 to ep for seamless normalized to + // denormalized transition. + inline static constexpr FPBits make_value(UIntType number, int ep) { + FPBits result; + // offset: +1 for sign, but -1 for implicit first bit + int lz = fputil::clz(number) - FloatProp::EXPONENT_WIDTH; + number <<= lz; + ep -= lz; + + if (likely(ep >= 0)) { + // Implicit number bit will be removed by mask + result.set_mantissa(number); + result.set_unbiased_exponent(ep + 1); + } else { + result.set_mantissa(number >> -ep); + } + return result; + } }; } // namespace fputil diff --git a/libc/src/__support/FPUtil/FloatProperties.h b/libc/src/__support/FPUtil/FloatProperties.h --- a/libc/src/__support/FPUtil/FloatProperties.h +++ b/libc/src/__support/FPUtil/FloatProperties.h @@ -22,7 +22,7 @@ static_assert(sizeof(BitsType) == sizeof(float), "Unexpected size of 'float' type."); - static constexpr uint32_t BIT_WIDTH = sizeof(BitsType) << 3; + static constexpr uint32_t BIT_WIDTH = sizeof(BitsType) * 8; static constexpr uint32_t MANTISSA_WIDTH = 23; static constexpr uint32_t EXPONENT_WIDTH = 8; @@ -43,7 +43,7 @@ static_assert(sizeof(BitsType) == sizeof(double), "Unexpected size of 'double' type."); - static constexpr uint32_t BIT_WIDTH = sizeof(BitsType) << 3; + static constexpr uint32_t BIT_WIDTH = sizeof(BitsType) * 8; static constexpr uint32_t MANTISSA_WIDTH = 52; static constexpr uint32_t EXPONENT_WIDTH = 11; @@ -95,7 +95,7 @@ static_assert(sizeof(BitsType) == sizeof(long double), "Unexpected size of 'long double' type."); - static constexpr uint32_t BIT_WIDTH = (sizeof(BitsType) << 3) - 48; + static constexpr uint32_t BIT_WIDTH = (sizeof(BitsType) * 8) - 48; static constexpr uint32_t MANTISSA_WIDTH = 63; static constexpr uint32_t EXPONENT_WIDTH = 15; diff --git a/libc/src/__support/FPUtil/builtin_wrappers.h b/libc/src/__support/FPUtil/builtin_wrappers.h new file mode 100644 --- /dev/null +++ b/libc/src/__support/FPUtil/builtin_wrappers.h @@ -0,0 +1,61 @@ +//===--Convenient template for builtins -------------------------*- C++ -*-===// +// (Count Lead Zeroes) and (Count Trailing Zeros) +// +// 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_FPUTIL_CLZ_H +#define LLVM_LIBC_SRC_SUPPORT_FPUTIL_CLZ_H + +namespace __llvm_libc { +namespace fputil { + +// The following overloads are matched based on what is accepted by +// __builtin_clz/ctz* rather than using the exactly-sized aliases from stdint.h. +// This way, we can avoid making any assumptions about integer sizes and let the +// compiler match for us. +template static inline int clz(T val); +template <> inline int clz(unsigned int val) { + return __builtin_clz(val); +} +template <> inline int clz(unsigned long int val) { + return __builtin_clzl(val); +} +template <> inline int clz(unsigned long long int val) { + return __builtin_clzll(val); +} + +template static inline int ctz(T val); +template <> inline int ctz(unsigned int val) { + return __builtin_ctz(val); +} +template <> inline int ctz(unsigned long int val) { + return __builtin_ctzl(val); +} +template <> inline int ctz(unsigned long long int val) { + return __builtin_ctzll(val); +} + +template static inline bool isnan(T val) { + return __builtin_isnan(val); +} + +template static inline bool isinf(T val) { + return __builtin_isinf(val); +} + +template static inline bool isfinite(T val) { + return __builtin_isfinite(val); +} + +inline float quiet_NaN(float) { return __builtin_nanf(""); } +inline double quiet_NaN(double) { return __builtin_nan(""); } +inline long double quiet_NaN(long double) { return __builtin_nanl(""); } + +} // namespace fputil +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_CLZ_H diff --git a/libc/src/__support/FPUtil/generic/CMakeLists.txt b/libc/src/__support/FPUtil/generic/CMakeLists.txt --- a/libc/src/__support/FPUtil/generic/CMakeLists.txt +++ b/libc/src/__support/FPUtil/generic/CMakeLists.txt @@ -10,3 +10,9 @@ HDRS FMA.h ) + +add_header_library( + fmod + HDRS + FMod.h +) diff --git a/libc/src/__support/FPUtil/generic/FMod.h b/libc/src/__support/FPUtil/generic/FMod.h new file mode 100644 --- /dev/null +++ b/libc/src/__support/FPUtil/generic/FMod.h @@ -0,0 +1,250 @@ +//===-- Common header for fmod implementations ------------------*- 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_FPUTIL_GENERIC_FMOD_H +#define LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_FMOD_H + +#include "src/__support/CPP/TypeTraits.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/builtin_wrappers.h" +#include "src/__support/common.h" +#include "src/math/generic/math_utils.h" + +#include + +namespace __llvm_libc { +namespace fputil { +namespace generic { + +// Objective: +// The algorithm uses integer arithmetic (max uint64_t) for general case. +// Some common cases, like abs(x) < abs(y) or abs(x) < 1000 * abs(y) are +// treated specially to increase performance. The part of checking special +// cases, numbers NaN, INF etc. is inspired by GCC libm implementation. +// +// Mathematics: +// Input x,y in the algorithm is represented (mathematically) like hx * 2^ix +// and hy * 2^iy. This is an unambiguous number representation. For example: +// h * 2^i = (2 * h) * 2^(i-1) +// The algorithm uses the fact that +// r = a % b = (a % (N * b)) % b, +// where N is positive integer number. a and b - positive. Let's adopt the +// formula for representation above. +// a = hx * 2^ix, b = hy * 2^iy, N = 2^k +// r(k) = a % b = (hx * 2^ix) % (2 ^k * hy * 2^iy) +// = 2^(iy + k) * (hx * 2^(ix - iy - k) % hy) +// r(k) = hr * 2^ir = (hx % hy) * 2^(iy + k) +// = (2^p * (hx % hy) * 2^(iy + k - p)) +// hr = 2^p * (hx % hy), ir = iy + k - p +// +// Algorithm description: +// The examples below use byte for simplicity. +// 1) Shift hy maximum to right without losing bits and increase iy value +// hy = 0b00101100 iy = 20 after shift hy = 0b00001011 iy = 22. +// 2) hx = hx % hy. +// 3) Move hx maximum to left. Note that after (hx = hx % hy) CLZ in hx is +// not +// lower than CLZ in hy. hx=0b00001001 ix = 100, hx=0b10010000, ix = 100-4 +// = 96. +// 4) Repeat (2) until ix == iy. +// +// Complexity analysis (double): +// Converting x,y to (hx,ix),(hy, iy): CTZ/shift/AND/OR/if. Loop count: +// (ix - iy) / (64 - "length of hy"). +// max("length of hy") = 53, +// max(ix - iy) = 2048 +// Maximum operation is 186. For rare "unrealistic" cases. +// +// Special cases (double): +// Supposing that case where |y| > 1e-292 and |x/y|<2000 is very common +// special processing is implemented. No hy alignment, no loop: +// result = (hx * 2^(ix - iy)) % hy. +// When x and y are both subnormal (rare case but...) the +// result = hx % hy. +// Simplified conversion back to double. + +template class FModGCCInspiredWrapper { + + static_assert(cpp::IsFloatingPointType::Value, + "FModGCCInspiredWrapper instantiated with invalid type."); + +public: + static bool PreCheck(T x, T y, T &out) { + if (unlikely(y == 0 || fputil::isnan(y) || !fputil::isfinite(x))) { + force_eval((x * y) / (y * x)); + // y=0,or x not finite or y is NaN + if (fputil::isnan(x) || fputil::isnan(y)) { + out = fputil::quiet_NaN(T(0)); + } else { + out = with_errno(quiet_NaN(T(0)), EDOM); + } + return true; + } + return false; + } +}; + +template class FModFastMathWrapper { + + static_assert(cpp::IsFloatingPointType::Value, + "FModFastMathWrapper instantiated with invalid type."); + +public: + static bool PreCheck(T, T, T &) { return false; } +}; + +template , + bool InverseMultiplication = false> +class FMod { + static_assert(cpp::IsFloatingPointType::Value, + "FMod instantiated with invalid type."); + +private: + using FPB = FPBits; + using intU_t = typename FPB::UIntType; + + template + typename std::enable_if< + !B, intU_t>::type inline constexpr static division_loop(int n, + int hyltzeroes, + intU_t hx, + intU_t hy) { + while (n > hyltzeroes) { + n -= hyltzeroes; + hx <<= hyltzeroes; + hx %= hy; + } + hx <<= n; + hx %= hy; + return hx; + } + + template + typename std::enable_if< + B, intU_t>::type inline constexpr static division_loop(int n, + int hyltzeroes, + intU_t hx, + intU_t hy) { + if (n > hyltzeroes) { + intU_t inv_hy = (std::numeric_limits::max() / hy); + while (n > hyltzeroes) { + n -= hyltzeroes; + intU_t hd = (hx * inv_hy) >> (FPB::FloatProp::BIT_WIDTH - hyltzeroes); + hx <<= hyltzeroes; + hx -= hd * hy; + while (unlikely(hx > hy)) + hx -= hy; + } + intU_t hd = (hx * inv_hy) >> (FPB::FloatProp::BIT_WIDTH - n); + hx <<= n; + hx -= hd * hy; + while (unlikely(hx > hy)) + hx -= hy; + } else { + hx <<= n; + hx %= hy; + } + return hx; + } + + inline static constexpr FPB make_internal(FPB sx, FPB sy) { + + if (likely(sx.uintval() <= sy.uintval())) { + if (sx.uintval() < sy.uintval()) + return sx; // |x|<|y| return x + return FPB::zero(); // |x|=|y| return 0.0 + } + + int ix = sx.get_unbiased_exponent(); + int iy = sy.get_unbiased_exponent(); + + // Most common case where |y| is "very normal" and |x/y| < 2^EXPONENT_WIDTH + if (likely(iy > int(FPB::FloatProp::MANTISSA_WIDTH) && + ix - iy <= int(FPB::FloatProp::EXPONENT_WIDTH))) { + intU_t hx = sx.get_explicit_mantissa(); + intU_t hy = sy.get_explicit_mantissa(); + intU_t d = (ix == iy) ? (hx - hy) : (hx << (ix - iy)) % hy; + if (d == 0) + return FPB::zero(); + // iy - 1 because of "zero power" for number with power 1 + return FPB::make_value(d, iy - 1); + } + /* Both subnormal special case. */ + if (unlikely(ix == 0 && iy == 0)) { + FPB d; + d.set_mantissa(sx.uintval() % sy.uintval()); + return d; + } + + // Note that hx is not subnormal by conditions above. + intU_t hx = sx.get_explicit_mantissa(); + ix--; + + intU_t hy = sy.get_explicit_mantissa(); + int lzhy = FPB::FloatProp::EXPONENT_WIDTH; + if (likely(iy > 0)) { + iy--; + } else { + hy = sy.get_mantissa(); + lzhy = clz(hy); + } + + // Assume hy != 0 + int tzhy = ctz(hy); + int hyltzeroes = lzhy + tzhy; + // n > 0 byt conditions above + int n = ix - iy; + + { + // Shift hy right until the end or n = 0 + int right_shift = n < tzhy ? n : tzhy; + hy >>= right_shift; + n -= right_shift; + iy += right_shift; + } + + { + // Shift hx left until the end or n = 0 + int left_shift = n < int(FPB::FloatProp::EXPONENT_WIDTH) + ? n + : FPB::FloatProp::EXPONENT_WIDTH; + hx <<= left_shift; + n -= left_shift; + } + + hx %= hy; + if (hx == 0) + return FPB::zero(); + + if (n == 0) + return FPB::make_value(hx, iy); + + /* hx next can't be 0, because hx < hy, hy % 2 == 1 hx * 2^i % hy != 0 */ + hx = division_loop(n, hyltzeroes, hx, hy); + return FPB::make_value(hx, iy); + } + +public: + static inline T make(T x, T y) { + if (T out; Wrapper::PreCheck(x, y, out)) + return out; + FPB sx(x), sy(y); + bool sign = sx.get_sign(); + sx.set_sign(false); + sy.set_sign(false); + FPB result = make_internal(sx, sy); + result.set_sign(sign); + return result.get_val(); + } +}; + +} // namespace generic +} // namespace fputil +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_FMOD_H diff --git a/libc/src/__support/FPUtil/generic/sqrt.h b/libc/src/__support/FPUtil/generic/sqrt.h --- a/libc/src/__support/FPUtil/generic/sqrt.h +++ b/libc/src/__support/FPUtil/generic/sqrt.h @@ -15,6 +15,7 @@ #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PlatformDefs.h" +#include "src/__support/FPUtil/builtin_wrappers.h" namespace __llvm_libc { namespace fputil { @@ -31,20 +32,7 @@ }; #endif // SPECIAL_X86_LONG_DOUBLE -// The following overloads are matched based on what is accepted by -// __builtin_clz* rather than using the exactly-sized aliases from stdint.h. -// This way, we can avoid making any assumptions about integer sizes and let the -// compiler match for us. -template static inline int clz(T val); -template <> inline int clz(unsigned int val) { - return __builtin_clz(val); -} -template <> inline int clz(unsigned long int val) { - return __builtin_clzl(val); -} -template <> inline int clz(unsigned long long int val) { - return __builtin_clzll(val); -} +using fputil::ctz; template static inline void normalize(int &exponent, diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -105,6 +105,9 @@ add_math_entrypoint_object(fminf) add_math_entrypoint_object(fminl) +add_math_entrypoint_object(fmod) +add_math_entrypoint_object(fmodf) + add_math_entrypoint_object(frexp) add_math_entrypoint_object(frexpf) add_math_entrypoint_object(frexpl) diff --git a/libc/src/math/fmod.h b/libc/src/math/fmod.h new file mode 100644 --- /dev/null +++ b/libc/src/math/fmod.h @@ -0,0 +1,18 @@ +//===-- Implementation header for fmod --------------------------*- 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_FMOD_H +#define LLVM_LIBC_SRC_MATH_FMOD_H + +namespace __llvm_libc { + +double fmod(double x, double y); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_FMOD_H diff --git a/libc/src/math/fmodf.h b/libc/src/math/fmodf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/fmodf.h @@ -0,0 +1,18 @@ +//===-- Implementation header for fmodf -------------------------*- 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_FMODF_H +#define LLVM_LIBC_SRC_MATH_FMODF_H + +namespace __llvm_libc { + +float fmodf(float x, float y); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_FMODF_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -1090,3 +1090,29 @@ COMPILE_OPTIONS -O3 ) + +add_entrypoint_object( + fmod + SRCS + fmod.cpp + HDRS + ../fmod.h + DEPENDS + libc.include.math + libc.src.__support.FPUtil.generic.fmod + COMPILE_OPTIONS + -O3 +) + +add_entrypoint_object( + fmodf + SRCS + fmodf.cpp + HDRS + ../fmodf.h + DEPENDS + libc.include.math + libc.src.__support.FPUtil.generic.fmod + COMPILE_OPTIONS + -O3 +) \ No newline at end of file diff --git a/libc/src/math/generic/fmod.cpp b/libc/src/math/generic/fmod.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/generic/fmod.cpp @@ -0,0 +1,19 @@ +//===-- Double-precision fmod 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/fmod.h" +#include "src/__support/FPUtil/generic/FMod.h" +#include "src/__support/common.h" + +namespace __llvm_libc { + +LLVM_LIBC_FUNCTION(double, fmod, (double x, double y)) { + return fputil::generic::FMod::make(x, y); +} + +} // namespace __llvm_libc diff --git a/libc/src/math/generic/fmodf.cpp b/libc/src/math/generic/fmodf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/generic/fmodf.cpp @@ -0,0 +1,19 @@ +//===-- Single-precision fmodf 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/fmodf.h" +#include "src/__support/FPUtil/generic/FMod.h" +#include "src/__support/common.h" + +namespace __llvm_libc { + +LLVM_LIBC_FUNCTION(float, fmodf, (float x, float y)) { + return fputil::generic::FMod::make(x, y); +} + +} // namespace __llvm_libc diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -1267,6 +1267,34 @@ libc.src.__support.FPUtil.fputil ) +add_fp_unittest( + fmodf_test + SUITE + libc_math_unittests + SRCS + fmodf_test.cpp + HDRS + FModTest.h + DEPENDS + libc.include.math + libc.src.math.fmodf + libc.src.__support.FPUtil.fputil +) + +add_fp_unittest( + fmod_test + SUITE + libc_math_unittests + SRCS + fmod_test.cpp + HDRS + FModTest.h + DEPENDS + libc.include.math + libc.src.math.fmod + libc.src.__support.FPUtil.fputil +) + add_subdirectory(generic) add_subdirectory(exhaustive) add_subdirectory(differential_testing) diff --git a/libc/test/src/math/FModTest.h b/libc/test/src/math/FModTest.h new file mode 100644 --- /dev/null +++ b/libc/test/src/math/FModTest.h @@ -0,0 +1,270 @@ +//===-- Utility class to test fmod special numbers ------------------------===// +// +// 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_TEST_SRC_MATH_FMODTEST_H +#define LLVM_LIBC_TEST_SRC_MATH_FMODTEST_H + +#include "src/__support/FPUtil/BasicOperations.h" +#include "src/__support/FPUtil/NearestIntegerOperations.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" +#include "utils/UnitTest/Test.h" + +#include + +#define TEST_SPECIAL(x, y, expected, dom_err, expected_exception) \ + EXPECT_EQ(FPBits(T(expected)).uintval(), FPBits(f(x, y)).uintval()); \ + EXPECT_MATH_ERRNO((dom_err) ? EDOM : 0); \ + EXPECT_FP_EXCEPTION(expected_exception); \ + __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT) + +#define TEST_REGULAR(x, y, expected) TEST_SPECIAL(x, y, expected, false, 0) + +template class FmodTest : public __llvm_libc::testing::Test { + + DECLARE_SPECIAL_CONSTANTS(T) + +public: + typedef T (*FModFunc)(T, T); + + void testSpecialNumbers(FModFunc f) { + using nl = std::numeric_limits; + + // fmod (+0, y) == +0 for y != 0. + TEST_SPECIAL(0, 3, 0, false, 0); + TEST_SPECIAL(0, nl::denorm_min(), 0, false, 0); + TEST_SPECIAL(0, -nl::denorm_min(), 0, false, 0); + TEST_SPECIAL(0, nl::min(), 0, false, 0); + TEST_SPECIAL(0, -nl::min(), 0, false, 0); + TEST_SPECIAL(0, nl::max(), 0, false, 0); + TEST_SPECIAL(0, -nl::max(), 0, false, 0); + + // fmod (-0, y) == -0 for y != 0. + TEST_SPECIAL(neg_zero, 3, neg_zero, false, 0); + TEST_SPECIAL(neg_zero, nl::denorm_min(), neg_zero, false, 0); + TEST_SPECIAL(neg_zero, -nl::denorm_min(), neg_zero, false, 0); + TEST_SPECIAL(neg_zero, nl::min(), neg_zero, false, 0); + TEST_SPECIAL(neg_zero, -nl::min(), neg_zero, false, 0); + TEST_SPECIAL(neg_zero, nl::max(), neg_zero, false, 0); + TEST_SPECIAL(neg_zero, -nl::max(), neg_zero, false, 0); + + // fmod (+inf, y) == nl::quiet_NaN() plus invalid exception. + TEST_SPECIAL(inf, 3, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(inf, -1.1L, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(inf, 0, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(inf, neg_zero, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(inf, nl::denorm_min(), nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(inf, nl::min(), nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(inf, nl::max(), nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(inf, inf, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(inf, neg_inf, nl::quiet_NaN(), true, FE_INVALID); + + // fmod (-inf, y) == nl::quiet_NaN() plus invalid exception. + TEST_SPECIAL(neg_inf, 3, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_inf, -1.1L, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_inf, 0, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_inf, neg_zero, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_inf, nl::denorm_min(), nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_inf, nl::min(), nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_inf, nl::max(), nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_inf, inf, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_inf, neg_inf, nl::quiet_NaN(), true, FE_INVALID); + + // fmod (x, +0) == nl::quiet_NaN() plus invalid exception. + TEST_SPECIAL(3, 0, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(-1.1L, 0, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(0, 0, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_zero, 0, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(nl::denorm_min(), 0, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(nl::min(), 0, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(nl::max(), 0, nl::quiet_NaN(), true, FE_INVALID); + + // fmod (x, -0) == nl::quiet_NaN() plus invalid exception. + TEST_SPECIAL(3, neg_zero, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(-1.1L, neg_zero, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(0, neg_zero, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(neg_zero, neg_zero, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(nl::denorm_min(), neg_zero, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(nl::min(), neg_zero, nl::quiet_NaN(), true, FE_INVALID); + TEST_SPECIAL(nl::max(), neg_zero, nl::quiet_NaN(), true, FE_INVALID); + + // fmod (x, +inf) == x for x not infinite. + TEST_SPECIAL(0, inf, 0, false, 0); + TEST_SPECIAL(neg_zero, inf, neg_zero, false, 0); + TEST_SPECIAL(nl::denorm_min(), inf, nl::denorm_min(), false, 0); + TEST_SPECIAL(nl::min(), inf, nl::min(), false, 0); + TEST_SPECIAL(nl::max(), inf, nl::max(), false, 0); + TEST_SPECIAL(3.0, inf, 3.0, false, 0); + // fmod (x, -inf) == x for x not infinite. + TEST_SPECIAL(0, neg_inf, 0, false, 0); + TEST_SPECIAL(neg_zero, neg_inf, neg_zero, false, 0); + TEST_SPECIAL(nl::denorm_min(), neg_inf, nl::denorm_min(), false, 0); + TEST_SPECIAL(nl::min(), neg_inf, nl::min(), false, 0); + TEST_SPECIAL(nl::max(), neg_inf, nl::max(), false, 0); + TEST_SPECIAL(3.0, neg_inf, 3.0, false, 0); + + TEST_SPECIAL(0, nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(0, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(neg_zero, nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(neg_zero, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(1, nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(1, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(inf, nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(inf, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(neg_inf, nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(neg_inf, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(0, nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(0, -nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(neg_zero, nl::signaling_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(neg_zero, -nl::signaling_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(1, nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(1, -nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(inf, nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(inf, -nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(neg_inf, nl::signaling_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(neg_inf, -nl::signaling_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(nl::quiet_NaN(), 0, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(-nl::quiet_NaN(), 0, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(nl::quiet_NaN(), neg_zero, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(-nl::quiet_NaN(), neg_zero, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(nl::quiet_NaN(), 1, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(-nl::quiet_NaN(), 1, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(nl::quiet_NaN(), inf, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(-nl::quiet_NaN(), inf, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(nl::quiet_NaN(), neg_inf, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(-nl::quiet_NaN(), neg_inf, nl::quiet_NaN(), false, 0); + TEST_SPECIAL(nl::signaling_NaN(), 0, nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), 0, nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(nl::signaling_NaN(), neg_zero, nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), neg_zero, nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(nl::signaling_NaN(), 1, nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), 1, nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(nl::signaling_NaN(), inf, nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), inf, nl::quiet_NaN(), false, FE_INVALID); + TEST_SPECIAL(nl::signaling_NaN(), neg_inf, nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), neg_inf, nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(nl::quiet_NaN(), nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(nl::quiet_NaN(), -nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(-nl::quiet_NaN(), nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(-nl::quiet_NaN(), -nl::quiet_NaN(), nl::quiet_NaN(), false, 0); + TEST_SPECIAL(nl::quiet_NaN(), nl::signaling_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(nl::quiet_NaN(), -nl::signaling_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(-nl::quiet_NaN(), nl::signaling_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(-nl::quiet_NaN(), -nl::signaling_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(nl::signaling_NaN(), nl::quiet_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(nl::signaling_NaN(), -nl::quiet_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), nl::quiet_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), -nl::quiet_NaN(), nl::quiet_NaN(), false, + FE_INVALID); + TEST_SPECIAL(nl::signaling_NaN(), nl::signaling_NaN(), nl::quiet_NaN(), + false, FE_INVALID); + TEST_SPECIAL(nl::signaling_NaN(), -nl::signaling_NaN(), nl::quiet_NaN(), + false, FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), nl::signaling_NaN(), nl::quiet_NaN(), + false, FE_INVALID); + TEST_SPECIAL(-nl::signaling_NaN(), -nl::signaling_NaN(), nl::quiet_NaN(), + false, FE_INVALID); + + TEST_SPECIAL(6.5, 2.25L, 2.0L, false, 0); + TEST_SPECIAL(-6.5, 2.25L, -2.0L, false, 0); + TEST_SPECIAL(6.5, -2.25L, 2.0L, false, 0); + TEST_SPECIAL(-6.5, -2.25L, -2.0L, false, 0); + + TEST_SPECIAL(nl::max(), nl::max(), 0, false, 0); + TEST_SPECIAL(nl::max(), -nl::max(), 0, false, 0); + TEST_SPECIAL(nl::max(), nl::min(), 0, false, 0); + TEST_SPECIAL(nl::max(), -nl::min(), 0, false, 0); + TEST_SPECIAL(nl::max(), nl::denorm_min(), 0, false, 0); + TEST_SPECIAL(nl::max(), -nl::denorm_min(), 0, false, 0); + TEST_SPECIAL(-nl::max(), nl::max(), neg_zero, false, 0); + TEST_SPECIAL(-nl::max(), -nl::max(), neg_zero, false, 0); + TEST_SPECIAL(-nl::max(), nl::min(), neg_zero, false, 0); + TEST_SPECIAL(-nl::max(), -nl::min(), neg_zero, false, 0); + TEST_SPECIAL(-nl::max(), nl::denorm_min(), neg_zero, false, 0); + TEST_SPECIAL(-nl::max(), -nl::denorm_min(), neg_zero, false, 0); + + TEST_SPECIAL(nl::min(), nl::max(), nl::min(), false, 0); + TEST_SPECIAL(nl::min(), -nl::max(), nl::min(), false, 0); + TEST_SPECIAL(nl::min(), nl::min(), 0, false, 0); + TEST_SPECIAL(nl::min(), -nl::min(), 0, false, 0); + TEST_SPECIAL(nl::min(), nl::denorm_min(), 0, false, 0); + TEST_SPECIAL(nl::min(), -nl::denorm_min(), 0, false, 0); + TEST_SPECIAL(-nl::min(), nl::max(), -nl::min(), false, 0); + TEST_SPECIAL(-nl::min(), -nl::max(), -nl::min(), false, 0); + TEST_SPECIAL(-nl::min(), nl::min(), neg_zero, false, 0); + TEST_SPECIAL(-nl::min(), -nl::min(), neg_zero, false, 0); + TEST_SPECIAL(-nl::min(), nl::denorm_min(), neg_zero, false, 0); + TEST_SPECIAL(-nl::min(), -nl::denorm_min(), neg_zero, false, 0); + + TEST_SPECIAL(nl::denorm_min(), nl::max(), nl::denorm_min(), false, 0); + TEST_SPECIAL(nl::denorm_min(), -nl::max(), nl::denorm_min(), false, 0); + TEST_SPECIAL(nl::denorm_min(), nl::min(), nl::denorm_min(), false, 0); + TEST_SPECIAL(nl::denorm_min(), -nl::min(), nl::denorm_min(), false, 0); + TEST_SPECIAL(nl::denorm_min(), nl::denorm_min(), 0, false, 0); + TEST_SPECIAL(nl::denorm_min(), -nl::denorm_min(), 0, false, 0); + TEST_SPECIAL(-nl::denorm_min(), nl::max(), -nl::denorm_min(), false, 0); + TEST_SPECIAL(-nl::denorm_min(), -nl::max(), -nl::denorm_min(), false, 0); + TEST_SPECIAL(-nl::denorm_min(), nl::min(), -nl::denorm_min(), false, 0); + TEST_SPECIAL(-nl::denorm_min(), -nl::min(), -nl::denorm_min(), false, 0); + TEST_SPECIAL(-nl::denorm_min(), nl::denorm_min(), neg_zero, false, 0); + TEST_SPECIAL(-nl::denorm_min(), -nl::denorm_min(), neg_zero, false, 0); + } + + void testRegularExtreme(FModFunc f) { + + TEST_REGULAR(0x1p127L, 0x3p-149L, 0x1p-149L); + TEST_REGULAR(0x1p127L, -0x3p-149L, 0x1p-149L); + TEST_REGULAR(0x1p127L, 0x3p-148L, 0x1p-147L); + TEST_REGULAR(0x1p127L, -0x3p-148L, 0x1p-147L); + TEST_REGULAR(0x1p127L, 0x3p-126L, 0x1p-125L); + TEST_REGULAR(0x1p127L, -0x3p-126L, 0x1p-125L); + TEST_REGULAR(-0x1p127L, 0x3p-149L, -0x1p-149L); + TEST_REGULAR(-0x1p127L, -0x3p-149L, -0x1p-149L); + TEST_REGULAR(-0x1p127L, 0x3p-148L, -0x1p-147L); + TEST_REGULAR(-0x1p127L, -0x3p-148L, -0x1p-147L); + TEST_REGULAR(-0x1p127L, 0x3p-126L, -0x1p-125L); + TEST_REGULAR(-0x1p127L, -0x3p-126L, -0x1p-125L); + + if constexpr (sizeof(T) >= sizeof(double)) { + TEST_REGULAR(0x1p1023L, 0x3p-1074L, 0x1p-1073L); + TEST_REGULAR(0x1p1023L, -0x3p-1074L, 0x1p-1073L); + TEST_REGULAR(0x1p1023L, 0x3p-1073L, 0x1p-1073L); + TEST_REGULAR(0x1p1023L, -0x3p-1073L, 0x1p-1073L); + TEST_REGULAR(0x1p1023L, 0x3p-1022L, 0x1p-1021L); + TEST_REGULAR(0x1p1023L, -0x3p-1022L, 0x1p-1021L); + TEST_REGULAR(-0x1p1023L, 0x3p-1074L, -0x1p-1073L); + TEST_REGULAR(-0x1p1023L, -0x3p-1074L, -0x1p-1073L); + TEST_REGULAR(-0x1p1023L, 0x3p-1073L, -0x1p-1073L); + TEST_REGULAR(-0x1p1023L, -0x3p-1073L, -0x1p-1073L); + TEST_REGULAR(-0x1p1023L, 0x3p-1022L, -0x1p-1021L); + TEST_REGULAR(-0x1p1023L, -0x3p-1022L, -0x1p-1021L); + } + } +}; + +#define LIST_FMOD_TESTS(T, func) \ + using LlvmLibcFmodTest = FmodTest; \ + TEST_F(LlvmLibcFmodTest, SpecialNumbers) { testSpecialNumbers(&func); } \ + TEST_F(LlvmLibcFmodTest, RegularExtreme) { testRegularExtreme(&func); } + +#endif // LLVM_LIBC_TEST_SRC_MATH_FMODTEST_H diff --git a/libc/test/src/math/differential_testing/BinaryOpSingleOutputDiff.h b/libc/test/src/math/differential_testing/BinaryOpSingleOutputDiff.h --- a/libc/test/src/math/differential_testing/BinaryOpSingleOutputDiff.h +++ b/libc/test/src/math/differential_testing/BinaryOpSingleOutputDiff.h @@ -22,6 +22,40 @@ public: typedef T Func(T, T); + static uint64_t run_diff_in_range(Func myFunc, Func otherFunc, + UIntType startingBit, UIntType endingBit, + UIntType N, + testutils::OutputFileStream &log) { + uint64_t result = 0; + if (endingBit < startingBit) { + return result; + } + + UIntType step = (endingBit - startingBit) / N; + for (UIntType bitsX = startingBit, bitsY = endingBit;; + bitsX += step, bitsY -= step) { + T x = T(FPBits(bitsX)); + T y = T(FPBits(bitsY)); + FPBits myBits = FPBits(myFunc(x, y)); + FPBits otherBits = FPBits(otherFunc(x, y)); + if (myBits.uintval() != otherBits.uintval()) { + result++; + log << " Input: " << bitsX << ", " << bitsY << " (" << x << ", " + << y << ")\n" + << " My result: " << myBits.uintval() << " (" << myBits.get_val() + << ")\n" + << "Other result: " << otherBits.uintval() << " (" + << otherBits.get_val() << ")\n" + << '\n'; + } + + if (endingBit - bitsX < step) { + break; + } + } + return result; + } + static void run_perf_in_range(Func myFunc, Func otherFunc, UIntType startingBit, UIntType endingBit, UIntType N, testutils::OutputFileStream &log) { @@ -84,6 +118,26 @@ myFunc, otherFunc, /* startingBit= */ FPBits(T(0x1.0p-10)).uintval(), /* endingBit= */ FPBits(T(0x1.0p+10)).uintval(), 10'000'001, log); } + + static void run_diff(Func myFunc, Func otherFunc, const char *logFile) { + uint64_t diffCount = 0; + testutils::OutputFileStream log(logFile); + log << " Diff tests with inputs in denormal range:\n"; + diffCount += run_diff_in_range( + myFunc, otherFunc, /* startingBit= */ UIntType(0), + /* endingBit= */ FPBits::MAX_SUBNORMAL, 1'000'001, log); + log << "\n Diff tests with inputs in normal range:\n"; + diffCount += run_diff_in_range( + myFunc, otherFunc, /* startingBit= */ FPBits::MIN_NORMAL, + /* endingBit= */ FPBits::MAX_NORMAL, 100'000'001, log); + log << "\n Diff tests with inputs in normal range with exponents " + "close to each other:\n"; + diffCount += run_diff_in_range( + myFunc, otherFunc, /* startingBit= */ FPBits(T(0x1.0p-10)).uintval(), + /* endingBit= */ FPBits(T(0x1.0p+10)).uintval(), 10'000'001, log); + + log << "Total number of differing results: " << diffCount << '\n'; + } }; } // namespace testing diff --git a/libc/test/src/math/differential_testing/CMakeLists.txt b/libc/test/src/math/differential_testing/CMakeLists.txt --- a/libc/test/src/math/differential_testing/CMakeLists.txt +++ b/libc/test/src/math/differential_testing/CMakeLists.txt @@ -448,3 +448,43 @@ COMPILE_OPTIONS -fno-builtin ) + +add_diff_binary( + fmodf_diff + SRCS + fmodf_diff.cpp + DEPENDS + .single_input_single_output_diff + libc.src.math.fmodf +) + +add_diff_binary( + fmodf_perf + SRCS + fmodf_perf.cpp + DEPENDS + .single_input_single_output_diff + libc.src.math.fmodf + COMPILE_OPTIONS + -fno-builtin +) + +add_diff_binary( + fmod_diff + SRCS + fmod_diff.cpp + DEPENDS + .single_input_single_output_diff + libc.src.math.fmod +) + +add_diff_binary( + fmod_perf + SRCS + fmod_perf.cpp + DEPENDS + .single_input_single_output_diff + libc.src.math.fmod + COMPILE_OPTIONS + -fno-builtin +) diff --git a/libc/test/src/math/differential_testing/fmod_diff.cpp b/libc/test/src/math/differential_testing/fmod_diff.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/differential_testing/fmod_diff.cpp @@ -0,0 +1,15 @@ +//===-- Differential test for fmod ----------------------------------------===// +// +// 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 "BinaryOpSingleOutputDiff.h" + +#include "src/math/fmod.h" + +#include + +BINARY_OP_SINGLE_OUTPUT_DIFF(double, __llvm_libc::fmod, ::fmod, "fmod_diff.log") diff --git a/libc/test/src/math/differential_testing/fmod_perf.cpp b/libc/test/src/math/differential_testing/fmod_perf.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/differential_testing/fmod_perf.cpp @@ -0,0 +1,15 @@ +//===-- Differential test for fmod ----------------------------------------===// +// +// 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 "BinaryOpSingleOutputDiff.h" + +#include "src/math/fmod.h" + +#include + +BINARY_OP_SINGLE_OUTPUT_PERF(double, __llvm_libc::fmod, ::fmod, "fmod_perf.log") diff --git a/libc/test/src/math/differential_testing/fmodf_diff.cpp b/libc/test/src/math/differential_testing/fmodf_diff.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/differential_testing/fmodf_diff.cpp @@ -0,0 +1,16 @@ +//===-- Differential test for fmodf ---------------------------------------===// +// +// 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 "BinaryOpSingleOutputDiff.h" + +#include "src/math/fmodf.h" + +#include + +BINARY_OP_SINGLE_OUTPUT_DIFF(float, __llvm_libc::fmodf, ::fmodf, + "fmodf_diff.log") diff --git a/libc/test/src/math/differential_testing/fmodf_perf.cpp b/libc/test/src/math/differential_testing/fmodf_perf.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/differential_testing/fmodf_perf.cpp @@ -0,0 +1,16 @@ +//===-- Differential test for fmodf ---------------------------------------===// +// +// 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 "BinaryOpSingleOutputDiff.h" + +#include "src/math/fmodf.h" + +#include + +BINARY_OP_SINGLE_OUTPUT_PERF(float, __llvm_libc::fmodf, ::fmodf, + "fmodf_perf.log") diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt --- a/libc/test/src/math/exhaustive/CMakeLists.txt +++ b/libc/test/src/math/exhaustive/CMakeLists.txt @@ -184,3 +184,16 @@ LINK_LIBRARIES -lpthread ) + +add_fp_unittest( + FMod_test + NO_RUN_POSTBUILD + NEED_MPFR + SUITE + libc_math_exhaustive_tests + SRCS + FMod_test.cpp + DEPENDS + libc.src.__support.FPUtil.fputil + libc.src.__support.FPUtil.generic.fmod +) diff --git a/libc/test/src/math/exhaustive/FMod_test.cpp b/libc/test/src/math/exhaustive/FMod_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/exhaustive/FMod_test.cpp @@ -0,0 +1,71 @@ +//===-- Utility class to test FMod generic implementation -------*- 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/FPUtil/generic/FMod.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/FPMatcher.h" +#include "utils/UnitTest/Test.h" + +#include +#include + +namespace mpfr = __llvm_libc::testing::mpfr; + +template +class LlvmLibcFModTest : public __llvm_libc::testing::Test { + static constexpr std::array test_bases = { + T(0.0), + T(1.0), + T(3.0), + T(27.0), + T(11.0 / 8.0), + T(2.764443), + T(1.0) - std::numeric_limits::epsilon(), + T(1.0) + std::numeric_limits::epsilon(), + T(M_PI), + T(M_SQRT2), + T(M_E)}; + +public: + void testExtensive() { + using FMod = __llvm_libc::fputil::generic::FMod< + T, __llvm_libc::fputil::generic::FModFastMathWrapper, + InverseMultiplication>; + using nl = std::numeric_limits; + int min2 = nl::min_exponent - nl::digits - 5; + int max2 = nl::max_exponent + 3; + for (T by : test_bases) { + for (int iy = min2; iy < max2; iy++) { + T y = by * std::ldexp(2, iy); + if (y == 0 || !std::isfinite(y)) + continue; + for (T bx : test_bases) { + for (int ix = min2; ix < max2; ix++) { + T x = bx * std::ldexp(2, ix); + if (!std::isfinite(x)) + continue; + T result = FMod::make(x, y); + mpfr::BinaryInput input{x, y}; + EXPECT_MPFR_MATCH(mpfr::Operation::Fmod, input, result, 0.0); + } + } + } + } + } +}; + +using LlvmLibcFModFloatTest = LlvmLibcFModTest; +TEST_F(LlvmLibcFModFloatTest, ExtensiveTest) { testExtensive(); } + +using LlvmLibcFModFloatInvTest = LlvmLibcFModTest; +TEST_F(LlvmLibcFModFloatInvTest, ExtensiveTest) { testExtensive(); } + +using LlvmLibcFModDoubleTest = LlvmLibcFModTest; +TEST_F(LlvmLibcFModDoubleTest, ExtensiveTest) { testExtensive(); } + +using LlvmLibcFModDoubleInvTest = LlvmLibcFModTest; +TEST_F(LlvmLibcFModDoubleInvTest, ExtensiveTest) { testExtensive(); } diff --git a/libc/test/src/math/fmod_test.cpp b/libc/test/src/math/fmod_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/fmod_test.cpp @@ -0,0 +1,13 @@ +//===-- Unittests for fmod ------------------------------------------------===// +// +// 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 "FModTest.h" + +#include "src/math/fmod.h" + +LIST_FMOD_TESTS(double, __llvm_libc::fmod) diff --git a/libc/test/src/math/fmodf_test.cpp b/libc/test/src/math/fmodf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/fmodf_test.cpp @@ -0,0 +1,13 @@ +//===-- Unittests for fmodf -----------------------------------------------===// +// +// 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 "FModTest.h" + +#include "src/math/fmodf.h" + +LIST_FMOD_TESTS(float, __llvm_libc::fmodf) diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -56,6 +56,7 @@ // input and produce a single floating point number of the same type as // output. BeginBinaryOperationsSingleOutput, + Fmod, Hypot, EndBinaryOperationsSingleOutput, diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -247,6 +247,12 @@ return result; } + MPFRNumber fmod(const MPFRNumber &b) { + MPFRNumber result(*this); + mpfr_fmod(result.value, value, b.value, mpfr_rounding); + return result; + } + MPFRNumber frexp(int &exp) { MPFRNumber result(*this); mpfr_exp_t resultExp; @@ -561,6 +567,8 @@ MPFRNumber inputX(x, precision, rounding); MPFRNumber inputY(y, precision, rounding); switch (op) { + case Operation::Fmod: + return inputX.fmod(inputY); case Operation::Hypot: return inputX.hypot(inputY); default: