Skip to content

Commit eb9c197

Browse files
RinziiRinzii
committed
Finish implementation of remquo
1 parent 0c8b23e commit eb9c197

File tree

6 files changed

+418
-39
lines changed

6 files changed

+418
-39
lines changed

‎ccmath_headers.cmake‎

Lines changed: 28 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,7 @@
1+
##########################################
2+
# Internal headers
3+
##########################################
4+
15
set(ccmath_internal_helpers_headers
26
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/internal/helpers/bits.hpp
37
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/internal/helpers/endian.hpp
@@ -39,7 +43,18 @@ set(ccmath_internal_headers
3943
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/internal/version.hpp
4044
)
4145

46+
47+
##########################################
48+
# Detail headers
49+
##########################################
50+
51+
set(ccmath_detail_basic_impl_headers
52+
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/impl/remquo_float_impl.hpp
53+
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/impl/remquo_double_impl.hpp
54+
)
55+
4256
set(ccmath_detail_basic_headers
57+
${ccmath_detail_basic_impl_headers}
4358
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/abs.hpp
4459
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/fdim.hpp
4560
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/basic/fma.hpp
@@ -65,17 +80,17 @@ set(ccmath_detail_compare_headers
6580
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/compare/signbit.hpp
6681
)
6782

68-
set(ccmath_detail_exponential_details_headers
69-
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log_float_impl.hpp
70-
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log_double_impl.hpp
71-
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log_data.hpp
72-
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log2_float_impl.hpp
73-
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log2_double_impl.hpp
74-
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/details/log2_data.hpp
83+
set(ccmath_detail_exponential_impl_headers
84+
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log_float_impl.hpp
85+
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log_double_impl.hpp
86+
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log_data.hpp
87+
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log2_float_impl.hpp
88+
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log2_double_impl.hpp
89+
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/impl/log2_data.hpp
7590
)
7691

7792
set(ccmath_detail_exponential_headers
78-
${ccmath_detail_exponential_details_headers}
93+
${ccmath_detail_exponential_impl_headers}
7994
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/exp.hpp
8095
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/exp2.hpp
8196
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/detail/exponential/expm1.hpp
@@ -175,6 +190,11 @@ set(ccmath_detail_headers
175190
${ccmath_detail_root_headers}
176191
)
177192

193+
194+
##########################################
195+
# Root headers
196+
##########################################
197+
178198
set(ccmath_root_headers
179199
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/basic.hpp
180200
${CMAKE_CURRENT_SOURCE_DIR}/include/ccmath/compare.hpp

‎include/ccmath/detail/basic/impl/remquo_double_impl.hpp‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
/*
2+
* Copyright (c) 2024-Present Ian Pike
3+
* Copyright (c) 2024-Present ccmath contributors
4+
*
5+
* This library is provided under the MIT License.
6+
* See LICENSE for more information.
7+
*/
8+
9+
#pragma once
10+
11+
#include <cstdint>
12+
#include <limits>
13+
14+
#include "ccmath/internal/helpers/bits.hpp"
15+
#include "ccmath/internal/predef/unlikely.hpp"
16+
17+
#include "ccmath/detail/compare/isfinite.hpp"
18+
#include "ccmath/detail/compare/isnan.hpp"
19+
20+
#include "ccmath/detail/basic/abs.hpp"
21+
#include "ccmath/detail/basic/fmod.hpp"
22+
23+
namespace ccm::internal
24+
{
25+
namespace
26+
{
27+
namespace impl
28+
{
29+
inline constexpr double remquo_double_impl(double x, double y, int * quo) noexcept
30+
{
31+
// Assume all edge case checking is done before calling this function.
32+
std::int64_t x_i64{};
33+
std::int64_t y_i64{};
34+
std::uint64_t x_sign{};
35+
std::uint64_t quotient_sign{};
36+
int computed_quotient{};
37+
38+
x_i64 = ccm::helpers::double_to_int64(x);
39+
y_i64 = ccm::helpers::double_to_int64(y);
40+
41+
// Determine the signs of x and the quotient.
42+
x_sign = static_cast<std::uint64_t>(x_i64) & 0x8000000000000000ULL;
43+
quotient_sign = x_sign ^ (static_cast<std::uint64_t>(y_i64) & 0x8000000000000000ULL);
44+
45+
// Clear the sign bits from the int64_t representations of x and y.
46+
y_i64 &= 0x7fffffffffffffffULL;
47+
x_i64 &= 0x7fffffffffffffffULL;
48+
49+
// b (or bit 54) represents the highest bit we can compare against for both signed and unsigned integers without causing an overflow.
50+
// 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).
51+
if (y_i64 <= 0x7fbfffffffffffffULL)
52+
{
53+
x = ccm::fmod(x, 8 * y); // now x < (8 * y)
54+
}
55+
56+
if (CCM_UNLIKELY(x_i64 == y_i64))
57+
{
58+
*quo = quotient_sign ? -1 : 1;
59+
return 0.0 * x;
60+
}
61+
62+
x = ccm::fabs(x);
63+
y = ccm::helpers::int64_to_double(y_i64);
64+
computed_quotient = 0;
65+
66+
67+
if (y_i64 <= 0x7fcfffffffffffffULL && x >= 4 * y)
68+
{
69+
x -= 4 * y;
70+
computed_quotient += 4;
71+
}
72+
73+
if (y_i64 <= 0x7fdfffffffffffffULL && x >= 2 * y)
74+
{
75+
x -= 2 * y;
76+
computed_quotient += 2;
77+
}
78+
79+
if (y_i64 < 0x0020000000000000ULL)
80+
{
81+
if (x + x > y)
82+
{
83+
x -= y;
84+
++computed_quotient;
85+
if (x + x >= y)
86+
{
87+
x -= y;
88+
++computed_quotient;
89+
}
90+
}
91+
}
92+
else
93+
{
94+
double y_half = 0.5 * y;
95+
if (x > y_half)
96+
{
97+
x -= y;
98+
++computed_quotient;
99+
if (x >= y_half)
100+
{
101+
x -= y;
102+
++computed_quotient;
103+
}
104+
}
105+
}
106+
107+
*quo = quotient_sign ? -computed_quotient : computed_quotient;
108+
109+
// Make sure that the correct sign of zero results in round down mode.
110+
if (x == 0.0) { x = 0.0; }
111+
if (x_sign) { x = -x; }
112+
113+
return x;
114+
}
115+
} // namespace impl
116+
} // namespace
117+
118+
template <typename T, typename = std::enable_if_t<std::is_floating_point_v<T>>>
119+
inline constexpr T remquo_double(T x, T y, int* quo) noexcept
120+
{
121+
return static_cast<T>(impl::remquo_double_impl(x, y, quo));
122+
}
123+
} // namespace ccm::internal

‎include/ccmath/detail/basic/impl/remquo_float_impl.hpp‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,117 @@
1+
/*
2+
* Copyright (c) 2024-Present Ian Pike
3+
* Copyright (c) 2024-Present ccmath contributors
4+
*
5+
* This library is provided under the MIT License.
6+
* See LICENSE for more information.
7+
*/
8+
9+
#pragma once
10+
11+
#include <cstdint>
12+
#include <limits>
13+
14+
#include "ccmath/internal/helpers/bits.hpp"
15+
#include "ccmath/internal/predef/unlikely.hpp"
16+
17+
#include "ccmath/detail/compare/isfinite.hpp"
18+
#include "ccmath/detail/compare/isnan.hpp"
19+
20+
#include "ccmath/detail/basic/abs.hpp"
21+
#include "ccmath/detail/basic/fmod.hpp"
22+
23+
namespace ccm::internal
24+
{
25+
namespace
26+
{
27+
namespace impl
28+
{
29+
inline constexpr float remquo_float_impl(float x, float y, int * quo) noexcept
30+
{
31+
// Assume all edge case checking is done before calling this function.
32+
std::int32_t x_i32{};
33+
std::int32_t y_i32{};
34+
std::uint32_t x_sign{};
35+
int quotient_sign{};
36+
int computed_quotient{};
37+
38+
x_i32 = ccm::helpers::float_to_int32(x);
39+
y_i32 = ccm::helpers::float_to_int32(y);
40+
41+
// Determine the signs of x and the quotient.
42+
x_sign = static_cast<std::uint32_t>(x_i32) & 0x80000000;
43+
quotient_sign = x_sign ^ (static_cast<std::uint32_t>(y_i32) & 0x80000000);
44+
45+
// Clear the sign bits from the int32_t representations of x and y.
46+
y_i32 &= 0x7fffffff;
47+
x_i32 &= 0x7fffffff;
48+
49+
if (y_i32 <= 0x7dffffff)
50+
{
51+
x = ccm::fmod(x, 8 * y); // now x < (8 * y)
52+
}
53+
54+
if ((x_i32 - y_i32) == 0)
55+
{
56+
*quo = quotient_sign ? -1 : 1;
57+
return 0.0f * x;
58+
}
59+
60+
x = ccm::fabsf(x);
61+
y = ccm::fabsf(y);
62+
computed_quotient = 0;
63+
64+
if (y_i32 <= 0x7e7fffff && x >= 4 * y)
65+
{
66+
x -= 4 * y;
67+
computed_quotient += 4;
68+
}
69+
if (y_i32 <= 0x7effffff && x >= 2 * y)
70+
{
71+
x -= 2 * y;
72+
computed_quotient += 2;
73+
}
74+
if (y_i32 < 0x01000000)
75+
{
76+
if (x + x > y)
77+
{
78+
x -= y;
79+
++computed_quotient;
80+
if (x + x >= y)
81+
{
82+
x -= y;
83+
++computed_quotient;
84+
}
85+
}
86+
}
87+
else
88+
{
89+
float y_half = 0.5 * y;
90+
if (x > y_half)
91+
{
92+
x -= y;
93+
++computed_quotient;
94+
if (x >= y_half)
95+
{
96+
x -= y;
97+
++computed_quotient;
98+
}
99+
}
100+
}
101+
*quo = quotient_sign ? -computed_quotient : computed_quotient;
102+
103+
// Make sure that the correct sign of zero results in round down mode.
104+
if (x == 0.0f) { x = 0.0f; }
105+
if (x_sign) { x = -x; }
106+
107+
return x;
108+
}
109+
} // namespace impl
110+
} // namespace
111+
112+
template <typename T, typename = std::enable_if_t<std::is_floating_point_v<T>>>
113+
inline constexpr T remquo_float(T x, T y, int * quo) noexcept
114+
{
115+
return static_cast<T>(impl::remquo_float_impl(x, y, quo));
116+
}
117+
} // namespace ccm::internal

‎include/ccmath/detail/basic/remquo.hpp‎

Lines changed: 41 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -8,9 +8,8 @@
88

99
#pragma once
1010

11-
#include "ccmath/detail/basic/remainder.hpp"
12-
#include "ccmath/detail/compare/isinf.hpp"
13-
#include "ccmath/detail/compare/isnan.hpp"
11+
#include "ccmath/detail/basic/impl/remquo_double_impl.hpp"
12+
#include "ccmath/detail/basic/impl/remquo_float_impl.hpp"
1413

1514
namespace ccm
1615
{
@@ -23,38 +22,52 @@ namespace ccm
2322
* @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
2423
* least three of the least significant bits of x / y
2524
*
26-
* @warning This function is still extremely experimental and almost certainly not work as expected.
27-
* @note This function should work as expected with GCC 7.1 and later.
25+
* @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
26+
* constexpr and static_assert)
27+
*
28+
* @code
29+
* constexpr double get_remainder(double x, double y)
30+
* {
31+
* int quotient {0};
32+
* double remainder = ccm::remquo(x, y, &quotient);
33+
* return remainder;
34+
* }
35+
*
36+
* constexpr int get_quotient(double x, double y)
37+
* {
38+
* int quotient {0};
39+
* ccm::remquo(x, y, &quotient);
40+
* return quotient;
41+
* }
42+
* @endcode
2843
*/
29-
template <typename T>
44+
template <typename T, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
3045
inline constexpr T remquo(T x, T y, int * quo)
3146
{
32-
#if (defined(__GNUC__) && __GNUC__ >= 7 && __GNUC_MINOR__ >= 1)
33-
// Works with GCC 7.1
34-
// Not static_assert-able
35-
return __builtin_remquo(x, y, quo);
36-
#else
37-
// TODO: This function is a lot trickier to get working with constexpr.
38-
// I'm putting this on hold for now and not requiring it for v0.1.0.
39-
// Since remquo is not a commonly used function, I'm not going to worry about it for now.
40-
if constexpr (std::is_floating_point_v<T>)
47+
// We are not using __builtin_remquo with GCC due to a failure for it to work with static_assert
48+
// Our implementation does not have this issue.
49+
50+
// If x is ±∞ and y is not NaN, NaN is returned.
51+
if (CCM_UNLIKELY(ccm::isinf(x) && !ccm::isnan(y))) { return (x * y) / (x * y); }
52+
53+
// If y is ±0 and x is not NaN, NaN is returned.
54+
if (CCM_UNLIKELY(y == 0.0 && !ccm::isnan(x))) { return (x * y) / (x * y); }
55+
56+
// If either x or y is NaN, NaN is returned.
57+
if (CCM_UNLIKELY(ccm::isnan(x)))
4158
{
42-
// If x is ±∞ and y is not NaN, NaN is returned.
43-
// If y is ±0 and x is not NaN, NaN is returned.
44-
// If either argument is NaN, NaN is returned.
45-
if ((ccm::isinf(x) && !ccm::isnan(y)) || (y == 0 && !ccm::isnan(x)) || (ccm::isnan(x) || ccm::isnan(y)))
46-
{
47-
// All major compilers return -NaN.
48-
return -std::numeric_limits<T>::quiet_NaN();
49-
}
59+
if (ccm::signbit(x)) { return -std::numeric_limits<T>::quiet_NaN(); }
60+
else { return std::numeric_limits<T>::quiet_NaN(); }
5061
}
5162

52-
T r = ccm::remainder(x, y);
53-
// Having a lot of issues with handling the quo parameter. May use __builtin_bit_cast to handle this.
54-
//*quo = static_cast<int>(x / y) & ~(std::numeric_limits<int>::min)();
63+
if (CCM_UNLIKELY(ccm::isnan(y)))
64+
{
65+
if (ccm::signbit(y)) { return -std::numeric_limits<T>::quiet_NaN(); }
66+
else { return std::numeric_limits<T>::quiet_NaN(); }
67+
}
5568

56-
return r;
57-
#endif
69+
if constexpr (std::is_same_v<T, float>) { return ccm::internal::remquo_float(x, y, quo); }
70+
else { return ccm::internal::remquo_double(x, y, quo); }
5871
}
5972
} // namespace ccm
6073

0 commit comments

Comments
 (0)