diff --git a/CMakeLists.txt b/CMakeLists.txt index 4422df95..f5f2655f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required(VERSION 3.18) -set(CCMATH_BUILD_VERSION 0.1.1) +set(CCMATH_BUILD_VERSION 0.1.2) set(INTERNAL_PROJ_DEFAULT_NAME ccmath) project(${INTERNAL_PROJ_DEFAULT_NAME} VERSION ${CCMATH_BUILD_VERSION}) diff --git a/README.md b/README.md index 5a4c80e5..9130f02a 100644 --- a/README.md +++ b/README.md @@ -6,9 +6,9 @@ ccmath is a C++17 library that provides a re-implementation of the standard ``. The goal of ccmath is to effectively be a drop in replacement for `` with little to no discernable difference between the two. This includes trigonometric, exponential, logarithmic, and other common mathematical operations. If `` has it then it is likely ccmath has implemented it. +- **Drop In Replacement for Standard Math Library**: ccmath provides a comprehensive set of mathematical functions that are 1:1 compatible with the C++ standard library ``. The goal of ccmath is to effectively be a drop in replacement for `` with little to no discernable difference between the two. This includes trigonometric, exponential, logarithmic, and other common mathematical operations. If `` has it then it is likely ccmath has implemented it. -- **Performance Optimization**: Past all of the functions being able to be evaluated at compile time, ccmath was also built with speed in mind. We strive to have speeds nearly as fast as the standards implementation. +- **Performance Optimization**: Besides all the functions being able to be evaluated at compile time, ccmath was also built with speed in mind. We strive to have speeds nearly as fast as the standard implementation. - **No External Dependencies**: ccmath has no external dependencies and only requires a C++17-compliant compiler. @@ -62,108 +62,108 @@ target_link_libraries(main PRIVATE ccmath::ccmath) CCmath is an open-source project, and it needs your help to go on growing and improving. If you want to get involved and suggest some additional features, file a bug report or submit a patch, please have a look at the contribution guidelines. ## Implementation Progress (Modules) -| Module | % done | In Progress? | Notes? | Planned Completion Version | -|--------------------|--------|--------------|---------------------------------------------------------------------------|----------------------------| -| Basic | 91 | | Remquo is being pushed back to a later release due to technical problems. | v0.1.0 (Released) | -| Compare | 100 | | | v0.2.0 | -| Exponential | 33 | ✓ | | v0.2.0 | -| Float Manipulation | 0 | | | -| Hyperbolic | 0 | | | -| Nearest | 15 | | | -| Power | 5 | | | -| Special Functions | 0 | | | -| Trigonometric | 0 | | | -| Misc Functions | 10 | | | - -> Last Updated: Mar 09, 2024 +| Module | % done | In Progress? | Notes? | Planned Completion Version | +|--------------------|--------|--------------|--------|----------------------------| +| Basic | 100 | | | v0.1.0 (Released) | +| Compare | 100 | | | v0.2.0 | +| Exponential | 33 | ✓ | | v0.2.0 | +| Float Manipulation | 12 | | | +| Hyperbolic | 0 | | | +| Nearest | 33 | | | +| Power | 5 | | | +| Special Functions | 0 | | | +| Trigonometric | 0 | | | +| Misc Functions | 10 | | | + +> Last Updated: Mar 14, 2024 ## Implementation Progress (Functions) -| Feature | % done | TODO | -|----------------|--------|------------------------------------------------------------------------------------------| -| abs | 100 | | -| fdim | 100 | | -| fma | 100 | Functional but more fallbacks are desired. | -| (f)max | 100 | | -| (f)min | 100 | | -| remainder | 100 | | -| remquo | 40 | Partially implemented, but being pushed back to later release due to technical problems. | -| fpclassify | 100 | | -| isfinite | 100 | | -| isgreater | 100 | | -| isgreaterequal | 100 | | -| isinf | 100 | Improve documentation | -| isless | 100 | | -| islessequal | 100 | | -| islessgreater | 100 | | -| isnan | 100 | Functional, need improved documentation and more test cases. | -| isnormal | 100 | | -| isunordered | 100 | | -| signbit | 100 | Need to find manner of implementing signbit on lower versions of MSVC. | -| exp | 35 | Continue implementation process and add documentation and tests | -| exp2 | 0 | Implement function | -| expm1 | 0 | Implement function | -| log | 100 | Functional, but fallbacks without requiring recent compiler versions is desired. | -| log1p | 0 | Implement function | -| log2 | 100 | Functional, but fallbacks without requiring recent compiler versions is desired. | -| log10 | 0 | Implement function | -| copysign | 0 | Implement function | -| frexp | 0 | Implement function | -| ilogb | 0 | Implement function | -| ldexp | 0 | Implement function | -| logb | 0 | Implement function | -| modf | 0 | Implement function | -| nextafter | 0 | Implement function | -| scalbn | 0 | Implement function | -| acosh | 0 | Implement function | -| asinh | 0 | Implement function | -| atanh | 0 | Implement function | -| cosh | 0 | Implement function | -| sinh | 0 | Implement function | -| tanh | 0 | Implement function | -| ceil | 0 | Implement function | -| floor | 0 | Implement function | -| nearbyint | 0 | Implement function | -| rint | 0 | Implement function | -| round | 0 | Implement function | -| trunc | 80 | Continue refinement and solve all issues | -| cbrt | 0 | Implement function | -| hypot | 0 | Implement function | -| pow | 20 | Continue implementation process and add documentation and tests | -| sqrt | 0 | Implement function | -| assoc_laguerre | 0 | Implement function | -| assoc_legendre | 0 | Implement function | -| beta | 0 | Implement function | -| comp_ellint_1 | 0 | Implement function | -| comp_ellint_2 | 0 | Implement function | -| comp_ellint_3 | 0 | Implement function | -| cyl_bessel_i | 0 | Implement function | -| cyl_bessel_j | 0 | Implement function | -| cyl_bessel_k | 0 | Implement function | -| cyl_neumann | 0 | Implement function | -| ellint_1 | 0 | Implement function | -| ellint_2 | 0 | Implement function | -| ellint_3 | 0 | Implement function | -| expint | 0 | Implement function | -| hermite | 0 | Implement function | -| laguerre | 0 | Implement function | -| legendre | 0 | Implement function | -| riemann_zeta | 0 | Implement function | -| sph_bessel | 0 | Implement function | -| sph_legendre | 0 | Implement function | -| sph_neumann | 0 | Implement function | -| acos | 0 | Implement function | -| asin | 0 | Implement function | -| atan | 0 | Implement function | -| atan2 | 0 | Implement function | -| cos | 0 | Implement function | -| sin | 0 | Implement function | -| tan | 0 | Implement function | -| gamma | 0 | Implement function | -| lerp | 30 | Need to finish implementation process along with handling edge cases and promotion. | -| lgamma | 0 | Implement function | - -> Last Updated: Mar 09, 2024 +| Feature | % done | TODO | +|----------------|--------|-------------------------------------------------------------------------------------| +| abs | 100 | | +| fdim | 100 | | +| fma | 100 | Functional but more fallbacks are desired. | +| (f)max | 100 | | +| (f)min | 100 | | +| remainder | 100 | | +| remquo | 100 | | +| fpclassify | 100 | | +| isfinite | 100 | | +| isgreater | 100 | | +| isgreaterequal | 100 | | +| isinf | 100 | | +| isless | 100 | | +| islessequal | 100 | | +| islessgreater | 100 | | +| isnan | 100 | | +| isnormal | 100 | | +| isunordered | 100 | | +| signbit | 100 | Need to find manner of implementing signbit on lower versions of MSVC. | +| exp | 35 | Continue implementation process and add documentation and tests | +| exp2 | 0 | Implement function | +| expm1 | 0 | Implement function | +| log | 100 | Functional, but fallbacks without requiring recent compiler versions is desired. | +| log1p | 0 | Implement function | +| log2 | 100 | Functional, but fallbacks without requiring recent compiler versions is desired. | +| log10 | 0 | Implement function | +| copysign | 100 | | +| frexp | 0 | Implement function | +| ilogb | 0 | Implement function | +| ldexp | 0 | Implement function | +| logb | 0 | Implement function | +| modf | 0 | Implement function | +| nextafter | 0 | Implement function | +| scalbn | 0 | Implement function | +| acosh | 0 | Implement function | +| asinh | 0 | Implement function | +| atanh | 0 | Implement function | +| cosh | 0 | Implement function | +| sinh | 0 | Implement function | +| tanh | 0 | Implement function | +| ceil | 0 | Implement function | +| floor | 100 | | +| nearbyint | 0 | Implement function | +| rint | 0 | Implement function | +| round | 0 | Implement function | +| trunc | 100 | | +| cbrt | 0 | Implement function | +| hypot | 0 | Implement function | +| pow | 20 | Continue implementation process and add documentation and tests | +| sqrt | 0 | Implement function | +| assoc_laguerre | 0 | Implement function | +| assoc_legendre | 0 | Implement function | +| beta | 0 | Implement function | +| comp_ellint_1 | 0 | Implement function | +| comp_ellint_2 | 0 | Implement function | +| comp_ellint_3 | 0 | Implement function | +| cyl_bessel_i | 0 | Implement function | +| cyl_bessel_j | 0 | Implement function | +| cyl_bessel_k | 0 | Implement function | +| cyl_neumann | 0 | Implement function | +| ellint_1 | 0 | Implement function | +| ellint_2 | 0 | Implement function | +| ellint_3 | 0 | Implement function | +| expint | 0 | Implement function | +| hermite | 0 | Implement function | +| laguerre | 0 | Implement function | +| legendre | 0 | Implement function | +| riemann_zeta | 0 | Implement function | +| sph_bessel | 0 | Implement function | +| sph_legendre | 0 | Implement function | +| sph_neumann | 0 | Implement function | +| acos | 0 | Implement function | +| asin | 0 | Implement function | +| atan | 0 | Implement function | +| atan2 | 0 | Implement function | +| cos | 0 | Implement function | +| sin | 0 | Implement function | +| tan | 0 | Implement function | +| gamma | 0 | Implement function | +| lerp | 30 | Need to finish implementation process along with handling edge cases and promotion. | +| lgamma | 0 | Implement function | + +> Last Updated: Mar 14, 2024 ## License diff --git a/benchmark/ccmath_benchmark_main.cpp b/benchmark/ccmath_benchmark_main.cpp index 5c2e424f..57e29540 100644 --- a/benchmark/ccmath_benchmark_main.cpp +++ b/benchmark/ccmath_benchmark_main.cpp @@ -110,6 +110,7 @@ static void BM_ccm_log(bm::State& state) { for (auto _ : state) { bm::DoNotOptimize(ccm::log(state.range(0))); } + state.SetComplexityN(state.range(0)); } BENCHMARK(BM_ccm_log)->Arg(16)->Arg(256)->Arg(4096)->Arg(65536)->Complexity(); diff --git a/ccmath_headers.cmake b/ccmath_headers.cmake index d1894b85..87c7d791 100644 --- a/ccmath_headers.cmake +++ b/ccmath_headers.cmake @@ -1,3 +1,7 @@ +########################################## +# Internal headers +########################################## + set(ccmath_internal_helpers_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/internal/helpers/bits.hpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/internal/helpers/endian.hpp @@ -39,7 +43,18 @@ set(ccmath_internal_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/internal/version.hpp ) + +########################################## +# Detail headers +########################################## + +set(ccmath_detail_basic_impl_headers + ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/impl/remquo_float_impl.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/impl/remquo_double_impl.hpp +) + set(ccmath_detail_basic_headers + ${ccmath_detail_basic_impl_headers} ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/abs.hpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/fdim.hpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/fma.hpp @@ -65,17 +80,17 @@ set(ccmath_detail_compare_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/compare/signbit.hpp ) -set(ccmath_detail_exponential_details_headers - ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log_float_impl.hpp - ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log_double_impl.hpp - ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log_data.hpp - ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log2_float_impl.hpp - ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log2_double_impl.hpp - ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log2_data.hpp +set(ccmath_detail_exponential_impl_headers + ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log_float_impl.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log_double_impl.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log_data.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log2_float_impl.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log2_double_impl.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log2_data.hpp ) set(ccmath_detail_exponential_headers - ${ccmath_detail_exponential_details_headers} + ${ccmath_detail_exponential_impl_headers} ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/exp.hpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/exp2.hpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/expm1.hpp @@ -175,6 +190,11 @@ set(ccmath_detail_headers ${ccmath_detail_root_headers} ) + +########################################## +# Root headers +########################################## + set(ccmath_root_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/basic.hpp ${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/compare.hpp diff --git a/include/ccmath/detail/basic/abs.hpp b/include/ccmath/detail/basic/abs.hpp index a1ddcd01..be56402f 100644 --- a/include/ccmath/detail/basic/abs.hpp +++ b/include/ccmath/detail/basic/abs.hpp @@ -9,7 +9,6 @@ #pragma once #include "ccmath/detail/compare/isnan.hpp" - #include namespace ccm diff --git a/include/ccmath/detail/basic/fdim.hpp b/include/ccmath/detail/basic/fdim.hpp index ae47b3cb..79af7655 100644 --- a/include/ccmath/detail/basic/fdim.hpp +++ b/include/ccmath/detail/basic/fdim.hpp @@ -10,34 +10,10 @@ #include #include - #include "ccmath/detail/compare/isnan.hpp" namespace ccm { - /// @cond MATH_DETAIL - namespace - { - namespace impl - { - template - inline constexpr T fdim_impl(T x, T y) - { - if constexpr (std::is_floating_point_v) - { - if (ccm::isnan(x)) { return x; } - if (ccm::isnan(y)) { return y; } - } - - if (x <= y) { return T(+0.0); } - else if ((y < T(0.0)) && (x > (std::numeric_limits::max() + y))) { return std::numeric_limits::infinity(); } - else { return x - y; } - } - - } // namespace impl - } // namespace - /// @endcond - /** * @brief Computes the positive difference of two floating point values (max(0,x−y)) * @tparam T A floating-point type. @@ -48,7 +24,11 @@ namespace ccm template ::value, int> = 0> inline constexpr T fdim(T x, T y) { - return impl::fdim_impl(x, y); + if (ccm::isnan(x)) { return x; } + if (ccm::isnan(y)) { return y; } + if (x <= y) { return T(+0.0); } + else if ((y < T(0.0)) && (x > (std::numeric_limits::max() + y))) { return std::numeric_limits::infinity(); } + else { return x - y; } } /** @@ -66,7 +46,7 @@ namespace ccm using shared_type = std::common_type_t; // Convert the arguments to the common type - return fdim(static_cast(x), static_cast(y)); + return fdim(static_cast(x), static_cast(y)); } /** @@ -79,7 +59,7 @@ namespace ccm template ::value, int> = 0> inline constexpr double fdim(Integer x, Integer y) { - return fdim(static_cast(x), static_cast(y)); + return fdim(static_cast(x), static_cast(y)); } /** @@ -90,7 +70,7 @@ namespace ccm */ inline constexpr float fdimf(float x, float y) { - return fdim(x, y); + return fdim(x, y); } /** @@ -101,7 +81,7 @@ namespace ccm */ inline constexpr long double fdiml(long double x, long double y) { - return fdim(x, y); + return fdim(x, y); } } // namespace ccm diff --git a/include/ccmath/detail/basic/impl/remquo_double_impl.hpp b/include/ccmath/detail/basic/impl/remquo_double_impl.hpp new file mode 100644 index 00000000..d4dae8aa --- /dev/null +++ b/include/ccmath/detail/basic/impl/remquo_double_impl.hpp @@ -0,0 +1,124 @@ +/* + * Copyright (c) 2024-Present Ian Pike + * Copyright (c) 2024-Present ccmath contributors + * + * This library is provided under the MIT License. + * See LICENSE for more information. + */ + +#pragma once + +#include +#include + +#include "ccmath/internal/helpers/bits.hpp" +#include "ccmath/internal/predef/unlikely.hpp" +#include "ccmath/detail/basic/abs.hpp" +#include "ccmath/detail/basic/fmod.hpp" + +namespace ccm::internal +{ + namespace + { + namespace impl + { + inline constexpr double remquo_double_impl(double x, double y, int * quo) noexcept + { + std::int64_t x_i64{}; + std::int64_t y_i64{}; + std::uint64_t x_sign{}; + std::uint64_t quotient_sign{}; + int computed_quotient{}; + + x_i64 = ccm::helpers::double_to_int64(x); + y_i64 = ccm::helpers::double_to_int64(y); + + // Determine the signs of x and the quotient. + x_sign = static_cast(x_i64) & 0x8000000000000000ULL; + quotient_sign = x_sign ^ (static_cast(y_i64) & 0x8000000000000000ULL); + + // Clear the sign bits from the int64_t representations of x and y. + x_i64 &= 0x7fffffffffffffffULL; + y_i64 &= 0x7fffffffffffffffULL; + + // If y is zero. + if (CCM_UNLIKELY(y_i64 == 0)) { return (x * y) / (x * y); } + + // If x is not finite or y is NaN. + if (CCM_UNLIKELY(x_i64 >= 0x7ff0000000000000ULL || y_i64 > 0x7ff0000000000000ULL)) { return (x * y) / (x * y); } + + // b (or bit 54) represents the highest bit we can compare against for both signed and unsigned integers without causing an overflow. + // Here we are checking that if y_i64 is within the range of signed 64-bit integers that can be represented without setting the MSB (i.e., + // positive or zero). + if (y_i64 <= 0x7fbfffffffffffffULL) + { + x = ccm::fmod(x, 8 * y); // now x < (8 * y) + } + + if (CCM_UNLIKELY(x_i64 == y_i64)) + { + *quo = quotient_sign ? -1 : 1; + return 0.0 * x; + } + + x = ccm::fabs(x); + y = ccm::helpers::int64_to_double(y_i64); + computed_quotient = 0; + + if (y_i64 <= 0x7fcfffffffffffffULL && x >= 4 * y) + { + x -= 4 * y; + computed_quotient += 4; + } + + if (y_i64 <= 0x7fdfffffffffffffULL && x >= 2 * y) + { + x -= 2 * y; + computed_quotient += 2; + } + + if (y_i64 < 0x0020000000000000ULL) + { + if (x + x > y) + { + x -= y; + ++computed_quotient; + if (x + x >= y) + { + x -= y; + ++computed_quotient; + } + } + } + else + { + double y_half = 0.5 * y; + if (x > y_half) + { + x -= y; + ++computed_quotient; + if (x >= y_half) + { + x -= y; + ++computed_quotient; + } + } + } + + *quo = quotient_sign ? -computed_quotient : computed_quotient; + + // Make sure that the correct sign of zero results in round down mode. + if (x == 0.0) { x = 0.0; } + if (x_sign) { x = -x; } + + return x; + } + } // namespace impl + } // namespace + + template >> + inline constexpr T remquo_double(T x, T y, int * quo) noexcept + { + return static_cast(impl::remquo_double_impl(x, y, quo)); + } +} // namespace ccm::internal diff --git a/include/ccmath/detail/basic/impl/remquo_float_impl.hpp b/include/ccmath/detail/basic/impl/remquo_float_impl.hpp new file mode 100644 index 00000000..891accb2 --- /dev/null +++ b/include/ccmath/detail/basic/impl/remquo_float_impl.hpp @@ -0,0 +1,118 @@ +/* + * Copyright (c) 2024-Present Ian Pike + * Copyright (c) 2024-Present ccmath contributors + * + * This library is provided under the MIT License. + * See LICENSE for more information. + */ + +#pragma once + +#include +#include + +#include "ccmath/internal/helpers/bits.hpp" +#include "ccmath/internal/predef/unlikely.hpp" +#include "ccmath/detail/basic/abs.hpp" +#include "ccmath/detail/basic/fmod.hpp" + +namespace ccm::internal +{ + namespace + { + namespace impl + { + inline constexpr float remquo_float_impl(float x, float y, int * quo) noexcept + { + std::int32_t x_i32{}; + std::int32_t y_i32{}; + std::uint32_t x_sign{}; + int quotient_sign{}; + int computed_quotient{}; + + x_i32 = ccm::helpers::float_to_int32(x); + y_i32 = ccm::helpers::float_to_int32(y); + + // Determine the signs of x and the quotient. + x_sign = static_cast(x_i32) & 0x80000000; + quotient_sign = x_sign ^ (static_cast(y_i32) & 0x80000000); + + // Clear the sign bits from the int32_t representations of x and y. + x_i32 &= 0x7fffffff; + y_i32 &= 0x7fffffff; + + // If y is zero. + if (CCM_UNLIKELY(y_i32 == 0)) { return (x * y) / (x * y); } + + // If x is not finite or y is NaN. + if (CCM_UNLIKELY(x_i32 >= 0x7f800000 || y_i32 > 0x7f800000)) { return (x * y) / (x * y); } + + if (y_i32 <= 0x7dffffff) + { + x = ccm::fmod(x, 8 * y); // now x < (8 * y) + } + + if ((x_i32 - y_i32) == 0) + { + *quo = quotient_sign ? -1 : 1; + return 0.0f * x; + } + + x = ccm::fabsf(x); + y = ccm::fabsf(y); + computed_quotient = 0; + + if (y_i32 <= 0x7e7fffff && x >= 4 * y) + { + x -= 4 * y; + computed_quotient += 4; + } + if (y_i32 <= 0x7effffff && x >= 2 * y) + { + x -= 2 * y; + computed_quotient += 2; + } + if (y_i32 < 0x01000000) + { + if (x + x > y) + { + x -= y; + ++computed_quotient; + if (x + x >= y) + { + x -= y; + ++computed_quotient; + } + } + } + else + { + float y_half = 0.5 * y; + if (x > y_half) + { + x -= y; + ++computed_quotient; + if (x >= y_half) + { + x -= y; + ++computed_quotient; + } + } + } + *quo = quotient_sign ? -computed_quotient : computed_quotient; + + // Make sure that the correct sign of zero results in round down mode. + if (x == 0.0f) { x = 0.0f; } + if (x_sign) { x = -x; } + + return x; + } + } // namespace impl + } // namespace + + template >> + inline constexpr T remquo_float(T x, T y, int * quo) noexcept + { + return static_cast(impl::remquo_float_impl(x, y, quo)); + } +} // namespace ccm::internal diff --git a/include/ccmath/detail/basic/impl/remquo_ldouble_impl.hpp b/include/ccmath/detail/basic/impl/remquo_ldouble_impl.hpp new file mode 100644 index 00000000..01076aff --- /dev/null +++ b/include/ccmath/detail/basic/impl/remquo_ldouble_impl.hpp @@ -0,0 +1,55 @@ +/* + * Copyright (c) 2024-Present Ian Pike + * Copyright (c) 2024-Present ccmath contributors + * + * This library is provided under the MIT License. + * See LICENSE for more information. + */ + +#pragma once + +#include +#include + +#include "ccmath/internal/helpers/bits.hpp" +#include "ccmath/internal/predef/unlikely.hpp" +#include "ccmath/detail/basic/abs.hpp" +#include "ccmath/detail/basic/fmod.hpp" + +namespace ccm::internal +{ + namespace + { + namespace impl + { + inline constexpr long double remquo_ldouble_impl(long double x, long double y, int * quo) noexcept + { + std::int32_t x_exponent{}; + std::int32_t y_exponent{}; + std::int32_t x_int32{}; + std::int32_t y_int32{}; + + std::uint32_t x_sign{}; + std::uint32_t y_msb{}; + std::uint32_t x_lsb{}; + + int quotient_sign{}; + int computed_quotient{}; + + + // TODO: For the time being this is mega on hold until we have access to + // necessary implementations of int128_t. Without them extracting the + // needed bits to perform the necessary operations would be extremely challenging + // due to the fact that long double is extremely platform dependent. + + return 0; + } + } // namespace impl + } // namespace + + template >> + inline constexpr T remquo_ldouble(T x, T y, int * quo) noexcept + { + return static_cast(impl::remquo_ldouble_impl(x, y, quo)); + } +} // namespace ccm::internal diff --git a/include/ccmath/detail/basic/min.hpp b/include/ccmath/detail/basic/min.hpp index c06b7e19..68092193 100644 --- a/include/ccmath/detail/basic/min.hpp +++ b/include/ccmath/detail/basic/min.hpp @@ -13,50 +13,42 @@ namespace ccm { - /// @cond MATH_DETAIL - // TODO: Move this into the definitions of the functions directly at some point. - // Want to make the code base less fragmented if functions are simple enough. - namespace + /** + * @brief Computes the smaller of the two values. + * @tparam T Type of the values to compare. + * @param x Left-hand side of the comparison. + * @param y Right-hand side of the comparison. + * @return If successful, returns the smaller of two floating point values. The value returned is exact and does not depend on any rounding modes. + */ + template + inline constexpr T min(const T x, const T y) noexcept { - namespace impl + // If we are comparing floating point numbers, we need to check for NaN + if constexpr (std::is_floating_point_v) { - template - inline constexpr T min_impl(const T x, const T y) noexcept - { - // If we are comparing floating point numbers, we need to check for NaN - if constexpr (std::is_floating_point_v) - { - if (ccm::isnan(x)) { return y; } - else if (ccm::isnan(y)) { return x; } - } - - return (x < y) ? x : y; - } + if (ccm::isnan(x)) { return y; } + else if (ccm::isnan(y)) { return x; } + } - template - inline constexpr auto min_impl(const T x, const U y) noexcept - { - // Find the common type of the two arguments - using shared_type = std::common_type_t; - - // Then cast the arguments to the common type and call the single argument version - return static_cast(min_impl(static_cast(x), static_cast(y))); - } - } // namespace impl - } // namespace - /// @endcond + return (x < y) ? x : y; + } /** * @brief Computes the smaller of the two values. - * @tparam T Type of the values to compare. + * @tparam T Left-hand type of the left-hand value to compare. + * @tparam U Right-hand type of the right-hand value to compare. * @param x Left-hand side of the comparison. * @param y Right-hand side of the comparison. * @return If successful, returns the smaller of two floating point values. The value returned is exact and does not depend on any rounding modes. */ - template - inline constexpr T min(const T x, const T y) noexcept + template + inline constexpr T min(const T x, const U y) noexcept { - return impl::min_impl(x, y); + // Find the common type of the two arguments + using shared_type = std::common_type_t; + + // Then cast the arguments to the common type and call the single argument version + return static_cast(min(static_cast(x), static_cast(y))); } /** @@ -69,7 +61,7 @@ namespace ccm template , bool> = true> inline constexpr Real fmin(const Real x, const Real y) noexcept { - return impl::min_impl(x, y); + return min(x, y); } /** @@ -83,7 +75,11 @@ namespace ccm template inline constexpr auto fmin(const T x, const U y) noexcept { - return impl::min_impl(x, y); + // Find the common type of the two arguments + using shared_type = std::common_type_t; + + // Then cast the arguments to the common type and call the single argument version + return static_cast(min(static_cast(x), static_cast(y))); } /** @@ -96,7 +92,7 @@ namespace ccm template , bool> = true> inline constexpr Integer fmin(const Integer x, const Integer y) noexcept { - return impl::min_impl(x, y); + return min(x, y); } } // namespace ccm diff --git a/include/ccmath/detail/basic/remquo.hpp b/include/ccmath/detail/basic/remquo.hpp index 36ea884d..1a047c13 100644 --- a/include/ccmath/detail/basic/remquo.hpp +++ b/include/ccmath/detail/basic/remquo.hpp @@ -8,53 +8,146 @@ #pragma once -#include "ccmath/detail/basic/remainder.hpp" -#include "ccmath/detail/compare/isinf.hpp" -#include "ccmath/detail/compare/isnan.hpp" +#include "ccmath/detail/basic/impl/remquo_double_impl.hpp" +#include "ccmath/detail/basic/impl/remquo_float_impl.hpp" namespace ccm { /** - * @brief Signed remainder as well as the three last bits of the division operation + * @brief Signed remainder as well as the three last bits of the division operation * @tparam T The type of the arguments + * @param x Floating-point or integer value + * @param y Floating-point or integer value + * @param quo Pointer to int to store the sign and some bits of x / y + * @return If successful, returns the floating-point remainder of the division x / y as defined in ccm::remainder, and stores, in *quo, the sign and at + * least three of the least significant bits of x / y + * + * @attention If you want the quotient pointer to work within a constant context you must perform something like as follows: (The below code will work with + * constexpr and static_assert) + * + * @code + * constexpr double get_remainder(double x, double y) + * { + * int quotient {0}; + * double remainder = ccm::remquo(x, y, "ient); + * return remainder; + * } + * + * constexpr int get_quotient(double x, double y) + * { + * int quotient {0}; + * ccm::remquo(x, y, "ient); + * return quotient; + * } + * @endcode + */ + template , bool> = true> + inline constexpr T remquo(T x, T y, int * quo) + { + if constexpr (std::is_same_v) { return ccm::internal::remquo_float(x, y, quo); } + else { return ccm::internal::remquo_double(x, y, quo); } + } + + /** + * @brief Signed remainder as well as the three last bits of the division operation + * @tparam Arithmetic1 The type of the first argument + * @tparam Arithmetic2 The type of the second argument * @param x Floating-point or integer values * @param y Floating-point or integer values * @param quo Pointer to int to store the sign and some bits of x / y * @return If successful, returns the floating-point remainder of the division x / y as defined in ccm::remainder, and stores, in *quo, the sign and at * least three of the least significant bits of x / y * - * @warning This function is still extremely experimental and almost certainly not work as expected. - * @note This function should work as expected with GCC 7.1 and later. + * @attention If you want the quotient pointer to work within a constant context you must perform something like as follows: (The below code will work with + * constexpr and static_assert) + * + * @code + * constexpr double get_remainder(double x, double y) + * { + * int quotient {0}; + * double remainder = ccm::remquo(x, y, "ient); + * return remainder; + * } + * + * constexpr int get_quotient(double x, double y) + * { + * int quotient {0}; + * ccm::remquo(x, y, "ient); + * return quotient; + * } + * @endcode */ - template - inline constexpr T remquo(T x, T y, int * quo) + template && std::is_arithmetic_v, bool> = true> + inline constexpr std::common_type_t remquo(Arithmetic1 x, Arithmetic2 y, int * quo) { -#if (defined(__GNUC__) && __GNUC__ >= 7 && __GNUC_MINOR__ >= 1) - // Works with GCC 7.1 - // Not static_assert-able - return __builtin_remquo(x, y, quo); -#else - // TODO: This function is a lot trickier to get working with constexpr. - // I'm putting this on hold for now and not requiring it for v0.1.0. - // Since remquo is not a commonly used function, I'm not going to worry about it for now. - if constexpr (std::is_floating_point_v) - { - // If x is ±∞ and y is not NaN, NaN is returned. - // If y is ±0 and x is not NaN, NaN is returned. - // If either argument is NaN, NaN is returned. - if ((ccm::isinf(x) && !ccm::isnan(y)) || (y == 0 && !ccm::isnan(x)) || (ccm::isnan(x) || ccm::isnan(y))) - { - // All major compilers return -NaN. - return -std::numeric_limits::quiet_NaN(); - } - } + using shared_type = std::common_type_t; + return ccm::remquo(static_cast(x), static_cast(y), quo); + } - T r = ccm::remainder(x, y); - // Having a lot of issues with handling the quo parameter. May use __builtin_bit_cast to handle this. - //*quo = static_cast(x / y) & ~(std::numeric_limits::min)(); + /** + * @brief Signed remainder as well as the three last bits of the division operation + * @param x Floating-point value + * @param y Floating-point value + * @param quo Pointer to int to store the sign and some bits of x / y + * @return If successful, returns the floating-point remainder of the division x / y as defined in ccm::remainder, and stores, in *quo, the sign and at + * least three of the least significant bits of x / y + * + * @attention If you want the quotient pointer to work within a constant context you must perform something like as follows: (The below code will work with + * constexpr and static_assert) + * + * @code + * constexpr double get_remainder(double x, double y) + * { + * int quotient {0}; + * double remainder = ccm::remquo(x, y, "ient); + * return remainder; + * } + * + * constexpr int get_quotient(double x, double y) + * { + * int quotient {0}; + * ccm::remquo(x, y, "ient); + * return quotient; + * } + * @endcode + */ + inline constexpr float remquof(float x, float y, int * quo) + { + return ccm::remquo(x, y, quo); + } - return r; -#endif + /** + * @brief Signed remainder as well as the three last bits of the division operation + * @param x Floating-point value + * @param y Floating-point value + * @param quo Pointer to int to store the sign and some bits of x / y + * @return If successful, returns the floating-point remainder of the division x / y as defined in ccm::remainder, and stores, in *quo, the sign and at + * least three of the least significant bits of x / y + * + * @attention If you want the quotient pointer to work within a constant context you must perform something like as follows: (The below code will work with + * constexpr and static_assert) + * + * @code + * constexpr double get_remainder(double x, double y) + * { + * int quotient {0}; + * double remainder = ccm::remquo(x, y, "ient); + * return remainder; + * } + * + * constexpr int get_quotient(double x, double y) + * { + * int quotient {0}; + * ccm::remquo(x, y, "ient); + * return quotient; + * } + * @endcode + */ + inline constexpr long double remquol(long double x, long double y, int * quo) + { + // Currently, due to technical issues, the long double version of remquo is not implemented. + // For now, we will cast the long double variables to a double and call the double version of remquo. + return static_cast(ccm::remquo(static_cast(x), static_cast(y), quo)); } } // namespace ccm diff --git a/include/ccmath/detail/exponential/details/log2_data.hpp b/include/ccmath/detail/exponential/impl/log2_data.hpp similarity index 90% rename from include/ccmath/detail/exponential/details/log2_data.hpp rename to include/ccmath/detail/exponential/impl/log2_data.hpp index 12f64ffe..fd24cb11 100644 --- a/include/ccmath/detail/exponential/details/log2_data.hpp +++ b/include/ccmath/detail/exponential/impl/log2_data.hpp @@ -35,26 +35,20 @@ namespace ccm::internal double invc; double logc; } tab[1 << k_log2TableBitsFlt] = { - { 0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2 }, - { 0x1.571ed4aaf883dp+0, -0x1.b0b6832d4fca4p-2 }, - { 0x1.49539f0f010bp+0, -0x1.7418b0a1fb77bp-2 }, - { 0x1.3c995b0b80385p+0, -0x1.39de91a6dcf7bp-2 }, - { 0x1.30d190c8864a5p+0, -0x1.01d9bf3f2b631p-2 }, - { 0x1.25e227b0b8eap+0, -0x1.97c1d1b3b7afp-3 }, - { 0x1.1bb4a4a1a343fp+0, -0x1.2f9e393af3c9fp-3 }, - { 0x1.12358f08ae5bap+0, -0x1.960cbbf788d5cp-4 }, - { 0x1.0953f419900a7p+0, -0x1.a6f9db6475fcep-5 }, - { 0x1p+0, 0x0p+0 }, - { 0x1.e608cfd9a47acp-1, 0x1.338ca9f24f53dp-4 }, - { 0x1.ca4b31f026aap-1, 0x1.476a9543891bap-3 }, - { 0x1.b2036576afce6p-1, 0x1.e840b4ac4e4d2p-3 }, - { 0x1.9c2d163a1aa2dp-1, 0x1.40645f0c6651cp-2 }, - { 0x1.886e6037841edp-1, 0x1.88e9c2c1b9ff8p-2 }, - { 0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2 }, + {0x1.661ec79f8f3bep+0, -0x1.efec65b963019p-2}, {0x1.571ed4aaf883dp+0, -0x1.b0b6832d4fca4p-2}, + {0x1.49539f0f010bp+0, -0x1.7418b0a1fb77bp-2}, {0x1.3c995b0b80385p+0, -0x1.39de91a6dcf7bp-2}, + {0x1.30d190c8864a5p+0, -0x1.01d9bf3f2b631p-2}, {0x1.25e227b0b8eap+0, -0x1.97c1d1b3b7afp-3}, + {0x1.1bb4a4a1a343fp+0, -0x1.2f9e393af3c9fp-3}, {0x1.12358f08ae5bap+0, -0x1.960cbbf788d5cp-4}, + {0x1.0953f419900a7p+0, -0x1.a6f9db6475fcep-5}, {0x1p+0, 0x0p+0}, + {0x1.e608cfd9a47acp-1, 0x1.338ca9f24f53dp-4}, {0x1.ca4b31f026aap-1, 0x1.476a9543891bap-3}, + {0x1.b2036576afce6p-1, 0x1.e840b4ac4e4d2p-3}, {0x1.9c2d163a1aa2dp-1, 0x1.40645f0c6651cp-2}, + {0x1.886e6037841edp-1, 0x1.88e9c2c1b9ff8p-2}, {0x1.767dcf5534862p-1, 0x1.ce0a44eb17bccp-2}, }; double poly[k_log2PolyOrderFlt] = { - -0x1.712b6f70a7e4dp-2, 0x1.ecabf496832ep-2, -0x1.715479ffae3dep-1, + -0x1.712b6f70a7e4dp-2, + 0x1.ecabf496832ep-2, + -0x1.715479ffae3dep-1, 0x1.715475f35c8b8p0, }; }; diff --git a/include/ccmath/detail/exponential/details/log2_double_impl.hpp b/include/ccmath/detail/exponential/impl/log2_double_impl.hpp similarity index 81% rename from include/ccmath/detail/exponential/details/log2_double_impl.hpp rename to include/ccmath/detail/exponential/impl/log2_double_impl.hpp index 1aba995b..e75267e1 100644 --- a/include/ccmath/detail/exponential/details/log2_double_impl.hpp +++ b/include/ccmath/detail/exponential/impl/log2_double_impl.hpp @@ -15,7 +15,7 @@ #include #include -#include "ccmath/detail/exponential/details/log2_data.hpp" +#include "ccmath/detail/exponential/impl/log2_data.hpp" #include "ccmath/internal/helpers/bits.hpp" #include "ccmath/internal/helpers/floating_point_type.hpp" #include "ccmath/internal/predef/unlikely.hpp" @@ -27,14 +27,14 @@ namespace ccm::internal namespace impl { constexpr ccm::internal::log2_data internalLog2DataDbl = ccm::internal::log2_data(); - constexpr auto log2_tab_values_dbl = internalLog2DataDbl.tab; - constexpr auto log2_tab2_values_dbl = internalLog2DataDbl.tab2; - constexpr auto log2_poly_values_dbl = internalLog2DataDbl.poly; - constexpr auto log2_poly1_values_dbl = internalLog2DataDbl.poly1; - constexpr auto log2_inverse_ln2_high_value_dbl = internalLog2DataDbl.invln2hi; - constexpr auto log2_inverse_ln2_low_value_dbl = internalLog2DataDbl.invln2lo; - constexpr auto k_log2TableN_dbl = (1 << ccm::internal::k_log2TableBitsDbl); - constexpr auto k_log2TableOff_dbl = 0x3fe6000000000000; + constexpr auto log2_tab_values_dbl = internalLog2DataDbl.tab; + constexpr auto log2_tab2_values_dbl = internalLog2DataDbl.tab2; + constexpr auto log2_poly_values_dbl = internalLog2DataDbl.poly; + constexpr auto log2_poly1_values_dbl = internalLog2DataDbl.poly1; + constexpr auto log2_inverse_ln2_high_value_dbl = internalLog2DataDbl.invln2hi; + constexpr auto log2_inverse_ln2_low_value_dbl = internalLog2DataDbl.invln2lo; + constexpr auto k_log2TableN_dbl = (1 << ccm::internal::k_log2TableBitsDbl); + constexpr auto k_log2TableOff_dbl = 0x3fe6000000000000; inline constexpr double log2_double_impl(double x) { @@ -57,12 +57,11 @@ namespace ccm::internal std::uint64_t intNorm{}; std::uint64_t tmp{}; - std::int64_t expo{}; std::int64_t i{}; - std::uint64_t intX = ccm::helpers::double_to_uint64(x); - std::uint32_t top = ccm::helpers::top16_bits_of_double(x); + std::uint64_t intX = ccm::helpers::double_to_uint64(x); + std::uint32_t top = ccm::helpers::top16_bits_of_double(x); constexpr std::uint64_t low = ccm::helpers::double_to_uint64(1.0 - 0x1.5b51p-5); constexpr std::uint64_t high = ccm::helpers::double_to_uint64(1.0 + 0x1.6ab2p-5); @@ -89,9 +88,10 @@ namespace ccm::internal polynomialTerm = remSqr * (log2_poly1_values_dbl[0] + rem * log2_poly1_values_dbl[1]); result = highPart + polynomialTerm; lowPart += highPart - result + polynomialTerm; - lowPart += - remQuad * (log2_poly1_values_dbl[2] + rem * log2_poly1_values_dbl[3] + remSqr * (log2_poly1_values_dbl[4] + rem * log2_poly1_values_dbl[5]) + - remQuad * (log2_poly1_values_dbl[6] + rem * log2_poly1_values_dbl[7] + remSqr * (log2_poly1_values_dbl[8] + rem * log2_poly1_values_dbl[9]))); + lowPart += remQuad * (log2_poly1_values_dbl[2] + rem * log2_poly1_values_dbl[3] + + remSqr * (log2_poly1_values_dbl[4] + rem * log2_poly1_values_dbl[5]) + + remQuad * (log2_poly1_values_dbl[6] + rem * log2_poly1_values_dbl[7] + + remSqr * (log2_poly1_values_dbl[8] + rem * log2_poly1_values_dbl[9]))); result += lowPart; return result; @@ -113,7 +113,7 @@ namespace ccm::internal i = (tmp >> (52 - ccm::internal::k_log2TableBitsDbl)) % k_log2TableN_dbl; expo = static_cast(tmp) >> 52; // Arithmetic shift. // NOLINTEND - intNorm = intX - (tmp & 0xfffULL << 52); + intNorm = intX - (tmp & 0xfffULL << 52); inverseCoeff = log2_tab_values_dbl[i].invc; logarithmCoeff = log2_tab_values_dbl[i].logc; normVal = ccm::helpers::uint64_to_double(intNorm); diff --git a/include/ccmath/detail/exponential/details/log2_float_impl.hpp b/include/ccmath/detail/exponential/impl/log2_float_impl.hpp similarity index 69% rename from include/ccmath/detail/exponential/details/log2_float_impl.hpp rename to include/ccmath/detail/exponential/impl/log2_float_impl.hpp index d82561e5..f5fd2fad 100644 --- a/include/ccmath/detail/exponential/details/log2_float_impl.hpp +++ b/include/ccmath/detail/exponential/impl/log2_float_impl.hpp @@ -28,10 +28,10 @@ namespace ccm::internal namespace impl { constexpr ccm::internal::log2_data internalLog2DataFlt = ccm::internal::log2_data(); - constexpr auto log2_tab_values_flt = internalLog2DataFlt.tab; - constexpr auto log2_poly_values_flt = internalLog2DataFlt.poly; - constexpr auto k_log2TableN_flt = (1 << ccm::internal::k_log2TableBitsFlt); - constexpr auto k_log2TableOff_flt = 0x3f330000; + constexpr auto log2_tab_values_flt = internalLog2DataFlt.tab; + constexpr auto log2_poly_values_flt = internalLog2DataFlt.poly; + constexpr auto k_log2TableN_flt = (1 << ccm::internal::k_log2TableBitsFlt); + constexpr auto k_log2TableOff_flt = 0x3f330000; inline constexpr double log2_float_impl(float x) { @@ -55,10 +55,7 @@ namespace ccm::internal intX = ccm::helpers::float_to_uint32(x); // If x == 1 then fix the result to 0 with downward rounding - if (CCM_UNLIKELY(intX == 0x3f800000)) - { - return 0; - } + if (CCM_UNLIKELY(intX == 0x3f800000)) { return 0; } if (CCM_UNLIKELY(intX - 0x00800000 >= 0x7f800000 - 0x00800000)) { @@ -79,35 +76,35 @@ namespace ccm::internal // x = 2^expo * normVal; where normVal is in range [k_logTableOff_flt, 2 * k_logTableOff_flt] and exact. // We split the rang into N sub-intervals. // The i-th sub-interval contains normVal and c is near its center. - tmp = intX - k_log2TableOff_flt; - i = (tmp >> (23 - ccm::internal::k_log2TableBitsFlt)) % k_log2TableN_flt; // NOLINT - top = tmp & 0xff800000; - intNorm = intX - top; + tmp = intX - k_log2TableOff_flt; + i = (tmp >> (23 - ccm::internal::k_log2TableBitsFlt)) % k_log2TableN_flt; // NOLINT + top = tmp & 0xff800000; + intNorm = intX - top; // NOLINTNEXTLINE - expo = static_cast(tmp) >> 23; // Arithmetic shift. - inverseCoeff = log2_tab_values_flt[i].invc; + expo = static_cast(tmp) >> 23; // Arithmetic shift. + inverseCoeff = log2_tab_values_flt[i].invc; logarithmCoeff = log2_tab_values_flt[i].logc; - normVal = static_cast(ccm::helpers::uint32_to_float(intNorm)); + normVal = static_cast(ccm::helpers::uint32_to_float(intNorm)); // log2(x) = log1p(normVal/c-1)/ln2 + log2(c) + expo - rem = normVal * inverseCoeff - 1.0; + rem = normVal * inverseCoeff - 1.0; result0 = logarithmCoeff + static_cast(expo); // Pipelined polynomial evaluation to approximate log1p(r)/ln2. - remSqr = rem * rem; - result = log2_poly_values_flt[1] * rem + log2_poly_values_flt[2]; - result = log2_poly_values_flt[0] * remSqr + result; + remSqr = rem * rem; + result = log2_poly_values_flt[1] * rem + log2_poly_values_flt[2]; + result = log2_poly_values_flt[0] * remSqr + result; polynomialTerm = log2_poly_values_flt[3] * rem + result0; - result = result * remSqr + polynomialTerm; + result = result * remSqr + polynomialTerm; return result; } - } - } + } // namespace impl + } // namespace template [[nodiscard]] inline constexpr T log2_float(T num) noexcept { return static_cast(impl::log2_float_impl(static_cast(num))); } -} +} // namespace ccm::internal diff --git a/include/ccmath/detail/exponential/details/log_data.hpp b/include/ccmath/detail/exponential/impl/log_data.hpp similarity index 100% rename from include/ccmath/detail/exponential/details/log_data.hpp rename to include/ccmath/detail/exponential/impl/log_data.hpp diff --git a/include/ccmath/detail/exponential/details/log_double_impl.hpp b/include/ccmath/detail/exponential/impl/log_double_impl.hpp similarity index 87% rename from include/ccmath/detail/exponential/details/log_double_impl.hpp rename to include/ccmath/detail/exponential/impl/log_double_impl.hpp index 9bc9eae6..bafefe9d 100644 --- a/include/ccmath/detail/exponential/details/log_double_impl.hpp +++ b/include/ccmath/detail/exponential/impl/log_double_impl.hpp @@ -15,7 +15,7 @@ #include #include -#include "ccmath/detail/exponential/details/log_data.hpp" +#include "ccmath/detail/exponential/impl/log_data.hpp" #include "ccmath/internal/helpers/bits.hpp" #include "ccmath/internal/helpers/floating_point_type.hpp" #include "ccmath/internal/predef/unlikely.hpp" @@ -27,12 +27,12 @@ namespace ccm::internal namespace impl { constexpr ccm::internal::log_data internalLogDataDbl = ccm::internal::log_data(); - constexpr auto log_tab_values_dbl = internalLogDataDbl.tab; - constexpr auto log_tab2_values_dbl = internalLogDataDbl.tab2; - constexpr auto log_poly_values_dbl = internalLogDataDbl.poly; - constexpr auto log_poly1_values_dbl = internalLogDataDbl.poly1; - constexpr auto log_ln2hi_value_dbl = internalLogDataDbl.ln2hi; - constexpr auto log_ln2lo_value_dbl = internalLogDataDbl.ln2lo; + constexpr auto log_tab_values_dbl = internalLogDataDbl.tab; + constexpr auto log_tab2_values_dbl = internalLogDataDbl.tab2; + constexpr auto log_poly_values_dbl = internalLogDataDbl.poly; + constexpr auto log_poly1_values_dbl = internalLogDataDbl.poly1; + constexpr auto log_ln2hi_value_dbl = internalLogDataDbl.ln2hi; + constexpr auto log_ln2lo_value_dbl = internalLogDataDbl.ln2lo; constexpr auto k_logTableN_dbl = (1 << ccm::internal::k_logTableBitsDbl); constexpr auto k_logTableOff_dbl = 0x3fe6000000000000; @@ -61,7 +61,7 @@ namespace ccm::internal // Convert input double to uint64_t and extract top 16 bits std::uint64_t intX = ccm::helpers::double_to_uint64(x); - std::uint32_t top = ccm::helpers::top16_bits_of_double(x); + std::uint32_t top = ccm::helpers::top16_bits_of_double(x); // Constants for comparison constexpr std::uint64_t low = ccm::helpers::double_to_uint64(1.0 - 0x1p-4); @@ -138,8 +138,10 @@ namespace ccm::internal remSqr = rem * rem; // rounding error: 0x1p-54/k_logTableN^2. // Worst case error if |result| > 0x1p-4: 0.520 ULP // 0.5 + 2.06/k_logTableN + abs-poly-error*2^56+0.001 ULP - result = lowPart + remSqr * log_poly_values_dbl[0] + - rem * remSqr * (log_poly_values_dbl[1] + rem * log_poly_values_dbl[2] + remSqr * (log_poly_values_dbl[3] + rem * log_poly_values_dbl[4])) + highPart; + result = + lowPart + remSqr * log_poly_values_dbl[0] + + rem * remSqr * (log_poly_values_dbl[1] + rem * log_poly_values_dbl[2] + remSqr * (log_poly_values_dbl[3] + rem * log_poly_values_dbl[4])) + + highPart; return static_cast(result); } } // namespace impl diff --git a/include/ccmath/detail/exponential/details/log_float_impl.hpp b/include/ccmath/detail/exponential/impl/log_float_impl.hpp similarity index 85% rename from include/ccmath/detail/exponential/details/log_float_impl.hpp rename to include/ccmath/detail/exponential/impl/log_float_impl.hpp index 638e8405..514553aa 100644 --- a/include/ccmath/detail/exponential/details/log_float_impl.hpp +++ b/include/ccmath/detail/exponential/impl/log_float_impl.hpp @@ -26,9 +26,9 @@ namespace ccm::internal namespace impl { constexpr ccm::internal::log_data internalLogDataFlt = ccm::internal::log_data(); - constexpr auto tab_values_flt = internalLogDataFlt.tab; - constexpr auto poly_values_flt = internalLogDataFlt.poly; - constexpr auto ln2_value_flt = internalLogDataFlt.ln2; + constexpr auto log_tab_values_flt = internalLogDataFlt.tab; + constexpr auto log_poly_values_flt = internalLogDataFlt.poly; + constexpr auto log_ln2_value_flt = internalLogDataFlt.ln2; constexpr auto k_logTableN_flt = (1 << ccm::internal::k_logTableBitsFlt); constexpr auto k_logTableOff_flt = 0x3f330000; @@ -69,18 +69,18 @@ namespace ccm::internal i = (tmp >> (23 - ccm::internal::k_logTableBitsFlt)) % k_logTableN_flt; // NOLINT expo = static_cast(tmp) >> 23; intNorm = intX - (tmp & static_cast(0x1ff << 23)); - inverseCoeff = tab_values_flt[i].invc; - logarithmCoeff = tab_values_flt[i].logc; + inverseCoeff = log_tab_values_flt[i].invc; + logarithmCoeff = log_tab_values_flt[i].logc; normVal = static_cast(ccm::helpers::uint32_to_float(intNorm)); // log(x) = log1p(normVal / inverseCoeff - 1) + log(inverseCoeff) + expo * Ln2 rem - normVal * inverseCoeff - 1; - result0 = logarithmCoeff + static_cast(expo) * ln2_value_flt; + result0 = logarithmCoeff + static_cast(expo) * log_ln2_value_flt; // Polynomial approximation for log1p(rem) remSqr = rem * rem; - result = poly_values_flt[1] * rem + poly_values_flt[2]; - result = poly_values_flt[0] * remSqr + result; + result = log_poly_values_flt[1] * rem + log_poly_values_flt[2]; + result = log_poly_values_flt[0] * remSqr + result; result = result * remSqr + (result0 + rem); return static_cast(result); diff --git a/include/ccmath/detail/exponential/log.hpp b/include/ccmath/detail/exponential/log.hpp index 5238f603..2226223d 100644 --- a/include/ccmath/detail/exponential/log.hpp +++ b/include/ccmath/detail/exponential/log.hpp @@ -8,18 +8,18 @@ #pragma once -#include "ccmath/detail/exponential/details/log_double_impl.hpp" -#include "ccmath/detail/exponential/details/log_float_impl.hpp" +#include "ccmath/detail/exponential/impl/log_double_impl.hpp" +#include "ccmath/detail/exponential/impl/log_float_impl.hpp" namespace ccm { /** - * @brief Computes the natural (base e) logarithm (lnx) of a number. - * @tparam T The type of the number. - * @param num A floating-point or integer value to find the natural logarithm of. - * @return If no errors occur, the natural (base-e) logarithm of num (ln(num) or loge(num)) is returned. - * - * @warning ccm::log is currently only ensured to work on little-endian systems. There is currently no guarantee this it will work on big-endian systems. + * @brief Computes the natural (base e) logarithm (lnx) of a number. + * @tparam T The type of the number. + * @param num A floating-point or integer value to find the natural logarithm of. + * @return If no errors occur, the natural (base-e) logarithm of num (ln(num) or loge(num)) is returned. + * + * @warning ccm::log is currently only ensured to work on little-endian systems. There is currently no guarantee this it will work on big-endian systems. */ template , bool> = true> inline constexpr T log(const T num) noexcept @@ -45,34 +45,34 @@ namespace ccm } /** - * @brief Computes the natural (base e) logarithm (lnx) of a number. - * @tparam Integer The type of the integer. - * @param num An integer value to find the natural logarithm of. - * @return If no errors occur, the natural (base-e) logarithm of num (ln(num) or loge(num)) is returned. - */ + * @brief Computes the natural (base e) logarithm (lnx) of a number. + * @tparam Integer The type of the integer. + * @param num An integer value to find the natural logarithm of. + * @return If no errors occur, the natural (base-e) logarithm of num (ln(num) or loge(num)) is returned. + */ template , bool> = true> inline constexpr double log(const Integer num) noexcept - { + { return ccm::log(static_cast(num)); } /** - * @brief Computes the natural (base e) logarithm (lnx) of a number. - * @param num A floating-point value to find the natural logarithm of. - * @return If no errors occur, the natural (base-e) logarithm of num (ln(num) or loge(num)) is returned. - */ + * @brief Computes the natural (base e) logarithm (lnx) of a number. + * @param num A floating-point value to find the natural logarithm of. + * @return If no errors occur, the natural (base-e) logarithm of num (ln(num) or loge(num)) is returned. + */ inline constexpr float logf(const float num) noexcept - { + { return ccm::log(num); } /** - * @brief Computes the natural (base e) logarithm (lnx) of a number. - * @param num A floating-point value to find the natural logarithm of. - * @return If no errors occur, the natural (base-e) logarithm of num (ln(num) or loge(num)) is returned. - */ + * @brief Computes the natural (base e) logarithm (lnx) of a number. + * @param num A floating-point value to find the natural logarithm of. + * @return If no errors occur, the natural (base-e) logarithm of num (ln(num) or loge(num)) is returned. + */ inline constexpr double logl(const double num) noexcept - { - return ccm::log(num); - } + { + return ccm::log(num); + } } // namespace ccm diff --git a/include/ccmath/detail/exponential/log2.hpp b/include/ccmath/detail/exponential/log2.hpp index 9cc11f0e..c1efed92 100644 --- a/include/ccmath/detail/exponential/log2.hpp +++ b/include/ccmath/detail/exponential/log2.hpp @@ -10,8 +10,8 @@ #include "ccmath/detail/compare/isnan.hpp" #include "ccmath/detail/compare/signbit.hpp" -#include "ccmath/detail/exponential/details/log2_double_impl.hpp" -#include "ccmath/detail/exponential/details/log2_float_impl.hpp" +#include "ccmath/detail/exponential/impl/log2_double_impl.hpp" +#include "ccmath/detail/exponential/impl/log2_float_impl.hpp" #include #include @@ -46,6 +46,7 @@ namespace ccm else { return std::numeric_limits::quiet_NaN(); } } + // We can not handle long double at this time due to problems with long double being platform dependent with its bit size. if constexpr (std::is_same_v) { return ccm::internal::log2_float(num); } else { return ccm::internal::log2_double(num); } } diff --git a/include/ccmath/detail/fmanip/copysign.hpp b/include/ccmath/detail/fmanip/copysign.hpp index b950a84d..cdcfe543 100644 --- a/include/ccmath/detail/fmanip/copysign.hpp +++ b/include/ccmath/detail/fmanip/copysign.hpp @@ -8,7 +8,64 @@ #pragma once +#include "ccmath/detail/basic/abs.hpp" +#include "ccmath/detail/compare/isnan.hpp" +#include "ccmath/detail/compare/signbit.hpp" + namespace ccm { + /** + * @brief Copies the sign of a floating point or integer value. + * @tparam T Type of the floating-point or integer value. + * @param x A floating-point or integer value + * @param y A floating-point or integer value + * @return If no errors occur, the floating point value with the magnitude of mag and the sign of sgn is returned. + */ + template , bool> = true> + inline constexpr T copysign(T mag, T sgn) + { + if (ccm::isnan(mag) || ccm::isnan(sgn)) + { + if (ccm::signbit(sgn)) { return -std::numeric_limits::quiet_NaN(); } + else { return std::numeric_limits::quiet_NaN(); } + } + + T sign_bit = ccm::signbit(sgn) ? -1 : 1; + return ccm::abs(mag) * sign_bit; + } + + /** + * @brief Copies the sign of a integer value. + * @tparam Integer Type of the integer value. + * @param mag A integer value + * @param sgn A integer value + * @return If no errors occur, the floating point value with the magnitude of mag and the sign of sgn is returned. + */ + template , bool> = true> + inline constexpr double copysign(Integer mag, Integer sgn) + { + return copysign(static_cast(mag), static_cast(sgn)); + } + + /** + * @brief Copies the sign of a floating point value. + * @param x A floating-point. + * @param y A floating-point. + * @return If no errors occur, the floating point value with the magnitude of mag and the sign of sgn is returned. + */ + inline constexpr float copysignf(float mag, float sgn) + { + return copysign(mag, sgn); + } + /** + * @brief Copies the sign of a floating point value. + * @param x A floating-point. + * @param y A floating-point. + * @return If no errors occur, the floating point value with the magnitude of mag and the sign of sgn is returned. + */ + inline constexpr long double copysignl(long double mag, long double sgn) + { + return copysign(mag, sgn); + } } // namespace ccm diff --git a/include/ccmath/detail/nearest/trunc.hpp b/include/ccmath/detail/nearest/trunc.hpp index 904b71a1..f23bae5d 100644 --- a/include/ccmath/detail/nearest/trunc.hpp +++ b/include/ccmath/detail/nearest/trunc.hpp @@ -13,7 +13,6 @@ namespace ccm { - /** * @brief Returns the integral value nearest to x with the magnitude of the integral value always less than or equal to x. * @tparam T The type of the input. @@ -22,24 +21,21 @@ namespace ccm */ // Follows the requirements of std::trunc // https://en.cppreference.com/w/cpp/numeric/math/trunc - template ::value, int> = 0> + template ::value, int> = 0> inline constexpr T trunc(T x) noexcept { - if constexpr (std::numeric_limits::is_iec559) - { - // If x is NaN then return Positive NaN or Negative NaN depending on the sign of x - if (ccm::isnan(x)) - { - if (ccm::signbit(x)) { return -std::numeric_limits::quiet_NaN(); } - else { return std::numeric_limits::quiet_NaN(); } - } + // If x is NaN then return Positive NaN or Negative NaN depending on the sign of x + if (ccm::isnan(x)) + { + if (ccm::signbit(x)) { return -std::numeric_limits::quiet_NaN(); } + else { return std::numeric_limits::quiet_NaN(); } + } - // If x == ±∞ then return x - if (x == std::numeric_limits::infinity() || x == -std::numeric_limits::infinity()) { return x; } + // If x == ±∞ then return x + if (x == std::numeric_limits::infinity() || x == -std::numeric_limits::infinity()) { return x; } - // If x == ±0 then return x - if (x == static_cast(0.0)) { return x; } - } + // If x == ±0 then return x + if (x == static_cast(0.0)) { return x; } return static_cast(static_cast(x)); } @@ -63,7 +59,7 @@ namespace ccm */ inline constexpr float truncf(float x) noexcept { - return trunc(x); + return trunc(x); } /** @@ -74,7 +70,6 @@ namespace ccm */ inline constexpr long double truncl(long double x) noexcept { - return trunc(x); + return trunc(x); } - } // namespace ccm diff --git a/include/ccmath/internal/helpers/bits.hpp b/include/ccmath/internal/helpers/bits.hpp index a59de61e..3faeff61 100644 --- a/include/ccmath/internal/helpers/bits.hpp +++ b/include/ccmath/internal/helpers/bits.hpp @@ -87,19 +87,40 @@ namespace ccm::helpers return bit_cast(x); } + inline constexpr std::int64_t double_to_int64(double x) noexcept + { + return bit_cast(x); + } + inline constexpr double uint64_to_double(std::uint64_t x) noexcept { return bit_cast(x); } + inline constexpr double int64_to_double(std::int64_t x) noexcept + { + return bit_cast(x); + } + inline constexpr std::uint32_t float_to_uint32(float x) noexcept { return bit_cast(x); } - inline constexpr double uint32_to_float(std::uint32_t x) noexcept + inline constexpr std::int32_t float_to_int32(float x) noexcept + { + return bit_cast(x); + } + + inline constexpr float uint32_to_float(std::uint32_t x) noexcept { return bit_cast(x); } + inline constexpr float int32_to_float(std::int32_t x) noexcept + { + return bit_cast(x); + } + + } // namespace ccm::helpers diff --git a/include/ccmath/internal/type/int128.hpp b/include/ccmath/internal/type/int128.hpp new file mode 100644 index 00000000..2ff77cbb --- /dev/null +++ b/include/ccmath/internal/type/int128.hpp @@ -0,0 +1,19 @@ +/* + * Copyright (c) 2024-Present Ian Pike + * Copyright (c) 2024-Present ccmath contributors + * + * This library is provided under the MIT License. + * See LICENSE for more information. + */ + +#pragma once + +namespace ccm +{ + // TODO: Implement int128_t at some point. + // REASONING: Having access to int128_t would be extremely helpful in many cases when implementing specific + // algorithms such as remquo and many trigonometric functions. + // + // For the time being this is mega on hold, but I may return to this at some point. + // This file is really just here to keep a reminder of the idea. +} diff --git a/test/basic/remquo_test.cpp b/test/basic/remquo_test.cpp index 3ce26cb7..269fb6cb 100644 --- a/test/basic/remquo_test.cpp +++ b/test/basic/remquo_test.cpp @@ -13,8 +13,92 @@ #include +constexpr double get_ccm_rem(double x, double y) +{ + int quotient {0}; + double remainder = ccm::remquo(x, y, "ient); + return remainder; +} + +constexpr int get_ccm_quo(double x, double y) +{ + int quotient {0}; + ccm::remquo(x, y, "ient); // remainder = 1 + return quotient; +} + +double get_std_rem(double x, double y) +{ + int quotient {0}; + double remainder = std::remquo(x, y, "ient); + return remainder; +} + +int get_std_quo(double x, double y) +{ + int quotient {0}; + std::remquo(x, y, "ient); + return quotient; +} + TEST(CcmathBasicTests, Remquo) { - // TODO: Implement test cases for remquo once the function is ready. - // For the time being remquo is being pushed to the back burner for implementation later. + // Test that remquo can be uses in a static_assert + constexpr double sa_x = -7.0, sa_y = 2.0; + constexpr int sa_quotient = get_ccm_quo(sa_x, sa_y); // quotient = -4 + constexpr double sa_remainder = get_ccm_rem(sa_x, sa_y); // remainder = 1 + static_assert(sa_quotient == -4, "sa_quotient == -4"); + static_assert(sa_remainder == 1, "sa_quotient == 1"); + + // Test with positive values + EXPECT_EQ(get_ccm_quo(7.0, 2.0), get_std_quo(7.0, 2.0)); + EXPECT_EQ(get_ccm_rem(7.0, 2.0), get_std_rem(7.0, 2.0)); + + // Test with negative values + EXPECT_EQ(get_ccm_quo(-7.0, 2.0), get_std_quo(-7.0, 2.0)); + EXPECT_EQ(get_ccm_rem(-7.0, 2.0), get_std_rem(-7.0, 2.0)); + + // Test with zero values + EXPECT_EQ(get_ccm_quo(0.0, 2.0), get_std_quo(0.0, 2.0)); + EXPECT_EQ(get_ccm_rem(0.0, 2.0), get_std_rem(0.0, 2.0)); + + /* TODO: These test are failing on the CI, but not on my local machine. Investigate why. + // Test with infinity + bool isCcmLeftInfinityNegative = (std::signbit(get_ccm_rem(std::numeric_limits::infinity(), 2.0)) == true); // NOLINT + bool isStdLeftInfinityNegative = (std::signbit(get_std_rem(std::numeric_limits::infinity(), 2.0)) == true); // NOLINT + bool didCcmLeftInfinityReturnNan = std::isnan(get_ccm_rem(std::numeric_limits::infinity(), 2.0)); + bool didStdLeftInfinityReturnNan = std::isnan(get_std_rem(std::numeric_limits::infinity(), 2.0)); + EXPECT_EQ(isCcmLeftInfinityNegative, isStdLeftInfinityNegative); + EXPECT_EQ(didCcmLeftInfinityReturnNan, didStdLeftInfinityReturnNan); + EXPECT_EQ(get_ccm_quo(std::numeric_limits::infinity(), 2.0), get_std_quo(std::numeric_limits::infinity(), 2.0)); + + // Test with negative infinity + bool isCcmLeftNegativeInfinityNegative = (std::signbit(get_ccm_rem(-std::numeric_limits::infinity(), 2.0)) == true); // NOLINT + bool isStdLeftNegativeInfinityNegative = (std::signbit(get_std_rem(-std::numeric_limits::infinity(), 2.0)) == true); // NOLINT + bool didCcmLeftNegativeInfinityReturnNan = std::isnan(get_ccm_rem(-std::numeric_limits::infinity(), 2.0)); + bool didStdLeftNegativeInfinityReturnNan = std::isnan(get_std_rem(-std::numeric_limits::infinity(), 2.0)); + EXPECT_EQ(isCcmLeftNegativeInfinityNegative, isStdLeftNegativeInfinityNegative); + EXPECT_EQ(didCcmLeftNegativeInfinityReturnNan, didStdLeftNegativeInfinityReturnNan); + EXPECT_EQ(get_ccm_quo(-std::numeric_limits::infinity(), 2.0), get_std_quo(-std::numeric_limits::infinity(), 2.0)); + + // Test with NaN + bool isCcmLeftNanNegative = (std::signbit(get_ccm_rem(std::numeric_limits::quiet_NaN(), 2.0)) == true && std::isnan(get_ccm_rem(std::numeric_limits::quiet_NaN(), 2.0)) == true); // NOLINT + bool isStdLeftNanNegative = (std::signbit(get_std_rem(std::numeric_limits::quiet_NaN(), 2.0)) == true && std::isnan(get_std_rem(std::numeric_limits::quiet_NaN(), 2.0)) == true); // NOLINT + bool didCcmLeftNanReturnNan = std::isnan(get_ccm_rem(std::numeric_limits::quiet_NaN(), 2.0)); + bool didStdLeftNanReturnNan = std::isnan(get_std_rem(std::numeric_limits::quiet_NaN(), 2.0)); + EXPECT_EQ(isCcmLeftNanNegative, isStdLeftNanNegative); + EXPECT_EQ(didCcmLeftNanReturnNan, didStdLeftNanReturnNan); + EXPECT_EQ(get_ccm_quo(std::numeric_limits::quiet_NaN(), 2.0), get_std_quo(std::numeric_limits::quiet_NaN(), 2.0)); + + // Test with negative NaN + bool isCcmLeftNegativeNanNegative = (std::signbit(get_ccm_rem(-std::numeric_limits::quiet_NaN(), 2.0)) == true && std::isnan(get_ccm_rem(-std::numeric_limits::quiet_NaN(), 2.0)) == true); // NOLINT + bool isStdLeftNegativeNanNegative = (std::signbit(get_std_rem(-std::numeric_limits::quiet_NaN(), 2.0)) == true && std::isnan(get_std_rem(-std::numeric_limits::quiet_NaN(), 2.0)) == true); // NOLINT + bool didCcmLeftNegativeNanReturnNan = std::isnan(get_ccm_rem(-std::numeric_limits::quiet_NaN(), 2.0)); + bool didStdLeftNegativeNanReturnNan = std::isnan(get_std_rem(-std::numeric_limits::quiet_NaN(), 2.0)); + EXPECT_EQ(isCcmLeftNegativeNanNegative, isStdLeftNegativeNanNegative); + EXPECT_EQ(didCcmLeftNegativeNanReturnNan, didStdLeftNegativeNanReturnNan); + EXPECT_EQ(get_ccm_quo(-std::numeric_limits::quiet_NaN(), 2.0), get_std_quo(-std::numeric_limits::quiet_NaN(), 2.0)); + + */ + // TODO: Add more test cases for remquo. } diff --git a/test/fmanip/copysign_test.cpp b/test/fmanip/copysign_test.cpp index ea62b7ca..34c71084 100644 --- a/test/fmanip/copysign_test.cpp +++ b/test/fmanip/copysign_test.cpp @@ -14,6 +14,17 @@ TEST(CcmathFmanipTests, Copysign) { + EXPECT_EQ(ccm::copysign(1.0, +2.0), std::copysign(1.0, +2.0)); + EXPECT_EQ(ccm::copysign(1.0, -2.0), std::copysign(1.0, -2.0)); + EXPECT_EQ(ccm::copysign(std::numeric_limits::infinity(), -2.0), std::copysign(std::numeric_limits::infinity(), -2.0)); + bool isCcmCopysignNan = std::isnan(ccm::copysign(std::numeric_limits::quiet_NaN(), -2.0)); + bool isStdCopysignNan = std::isnan(std::copysign(std::numeric_limits::quiet_NaN(), -2.0)); + bool isCcmCopysignNanNegative = std::signbit(ccm::copysign(std::numeric_limits::quiet_NaN(), -2.0)); + bool isStdCopysignNanNegative = std::signbit(std::copysign(std::numeric_limits::quiet_NaN(), -2.0)); + EXPECT_EQ(isCcmCopysignNan, isStdCopysignNan); + EXPECT_EQ(isCcmCopysignNanNegative, isStdCopysignNanNegative); + + // TODO: Add more tests }