Skip to content

Commit 5ac823a

Browse files
ckormanyosckormanyos
committed
Merge branch 'develop' into restore_promote
2 parents 30f16d3 + 3085424 commit 5ac823a

File tree

4 files changed

+170
-16
lines changed

4 files changed

+170
-16
lines changed

‎include/boost/decimal/detail/remove_trailing_zeros.hpp‎

Lines changed: 43 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -40,8 +40,13 @@ constexpr auto remove_trailing_zeros(std::uint32_t n) noexcept -> remove_trailin
4040
{
4141
std::size_t s {};
4242

43-
auto r = rotr<32>(n * UINT32_C(184254097), 4);
44-
auto b = r < UINT32_C(429497);
43+
auto r = rotr<32>(n * UINT32_C(15273505), 8);
44+
auto b = r < UINT32_C(43);
45+
s = s * 2U + static_cast<std::size_t>(b);
46+
n = b ? r : n;
47+
48+
r = rotr<32>(n * UINT32_C(184254097), 4);
49+
b = r < UINT32_C(429497);
4550
s = s * 2U + static_cast<std::size_t>(b);
4651
n = b ? r : n;
4752

@@ -62,8 +67,13 @@ constexpr auto remove_trailing_zeros(std::uint64_t n) noexcept -> remove_trailin
6267
{
6368
std::size_t s {};
6469

65-
auto r = rotr<64>(n * UINT64_C(28999941890838049), 8);
66-
auto b = r < UINT64_C(184467440738);
70+
auto r = rotr<64>(n * UINT64_C(230079197716545), 16);
71+
auto b = r < UINT64_C(1845);
72+
s = s * 2U + static_cast<std::size_t>(b);
73+
n = b ? r : n;
74+
75+
r = rotr<64>(n * UINT64_C(28999941890838049), 8);
76+
b = r < UINT64_C(184467440738);
6777
s = s * 2U + static_cast<std::size_t>(b);
6878
n = b ? r : n;
6979

@@ -85,22 +95,39 @@ constexpr auto remove_trailing_zeros(std::uint64_t n) noexcept -> remove_trailin
8595
return {n, s};
8696
}
8797

88-
// TODO(mborland): Make this better for the 2-word case
8998
constexpr auto remove_trailing_zeros(uint128 n) noexcept -> remove_trailing_zeros_return<uint128>
9099
{
91-
if (n.high == UINT64_C(0))
92-
{
93-
const auto temp {remove_trailing_zeros(n.low)};
94-
return {static_cast<uint128>(temp.trimmed_number), temp.number_of_removed_zeros};
95-
}
96-
97100
std::size_t s {};
98101

99-
while (n % 10 == 0)
100-
{
101-
n /= 10;
102-
++s;
103-
}
102+
auto r = rotr<128>(n * uint128(UINT64_C(0x62B42691AD836EB1), UINT64_C(0x16590F420A835081)), 32);
103+
auto b = r < uint128 {UINT64_C(0x0), UINT64_C(0x33EC48)};
104+
s = s * 2U + static_cast<std::size_t>(b);
105+
n = b ? r : n;
106+
107+
r = rotr<128>(n * uint128 {UINT64_C(0x3275305C1066), UINT64_C(0xE4A4D1417CD9A041)}, 16);
108+
b = r < uint128 {UINT64_C(0x734), UINT64_C(0xACA5F6226F0ADA62)};
109+
s = s * 2U + static_cast<std::size_t>(b);
110+
n = b ? r : n;
111+
112+
r = rotr<128>(n * uint128 {UINT64_C(0x6B7213EE9F5A78), UINT64_C(0xC767074B22E90E21)}, 8);
113+
b = r < uint128 {UINT64_C(0x2AF31DC461), UINT64_C(0x1873BF3F70834ACE)};
114+
s = s * 2U + static_cast<std::size_t>(b);
115+
n = b ? r : n;
116+
117+
r = rotr<128>(n * uint128 {UINT64_C(0x95182A9930BE0DE), UINT64_C(0xD288CE703AFB7E91)}, 4);
118+
b = r < uint128 {UINT64_C(0x68DB8BAC710CB), UINT64_C(0x295E9E1B089A0276)};
119+
s = s * 2U + static_cast<std::size_t>(b);
120+
n = b ? r : n;
121+
122+
r = rotr<128>(n * uint128 {UINT64_C(0x28F5C28F5C28F5C2), UINT64_C(0x8F5C28F5C28F5C29)}, 2);
123+
b = r < uint128 {UINT64_C(0x28F5C28F5C28F5C), UINT64_C(0x28F5C28F5C28F5C3)};
124+
s = s * 2U + static_cast<std::size_t>(b);
125+
n = b ? r : n;
126+
127+
r = rotr<128>(n * uint128 {UINT64_C(0xCCCCCCCCCCCCCCCC), UINT64_C(0xCCCCCCCCCCCCCCCD)}, 1);
128+
b = r < uint128 {UINT64_C(0x1999999999999999), UINT64_C(0x999999999999999A)};
129+
s = s * 2U + static_cast<std::size_t>(b);
130+
n = b ? r : n;
104131

105132
return {n, s};
106133
}

‎test/Jamfile‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -110,6 +110,7 @@ run test_log1p.cpp ;
110110
run test_pow.cpp ;
111111
run test_promotion.cpp ;
112112
run test_remainder_remquo.cpp ;
113+
run test_remove_trailing_zeros.cpp ;
113114
run test_sin_cos.cpp ;
114115
run test_sinh.cpp ;
115116
run test_snprintf.cpp ;

‎test/test_remove_trailing_zeros.cpp‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
// Copyright 2024 Matt Borland
2+
// Distributed under the Boost Software License, Version 1.0.
3+
// https://www.boost.org/LICENSE_1_0.txt
4+
5+
#include <boost/decimal.hpp>
6+
#include <boost/core/lightweight_test.hpp>
7+
#include <array>
8+
#include <limits>
9+
#include <cstdint>
10+
11+
template <typename T>
12+
void test()
13+
{
14+
constexpr std::array<std::uint64_t, 20> powers_of_10 =
15+
{{
16+
UINT64_C(1), UINT64_C(10), UINT64_C(100), UINT64_C(1000), UINT64_C(10000), UINT64_C(100000), UINT64_C(1000000),
17+
UINT64_C(10000000), UINT64_C(100000000), UINT64_C(1000000000), UINT64_C(10000000000), UINT64_C(100000000000),
18+
UINT64_C(1000000000000), UINT64_C(10000000000000), UINT64_C(100000000000000), UINT64_C(1000000000000000),
19+
UINT64_C(10000000000000000), UINT64_C(100000000000000000), UINT64_C(1000000000000000000), UINT64_C(10000000000000000000)
20+
}};
21+
22+
for (const auto& val : powers_of_10)
23+
{
24+
if (val < std::numeric_limits<T>::max())
25+
{
26+
const auto temp {boost::decimal::detail::remove_trailing_zeros(static_cast<T>(val))};
27+
if (!BOOST_TEST_EQ(temp.trimmed_number, T(1)))
28+
{
29+
// LCOV_EXCL_START
30+
std::cerr << "Input Number: " << val
31+
<< "\nOutput Number: " << temp.trimmed_number
32+
<< "\nZeros removed: " << temp.number_of_removed_zeros << std::endl;
33+
// LCOV_EXCL_STOP
34+
}
35+
}
36+
}
37+
}
38+
39+
void test_extended()
40+
{
41+
using namespace boost::decimal;
42+
constexpr std::array<detail::uint128, 18> powers_of_10 =
43+
{{
44+
detail::uint128 {UINT64_C(0x5), UINT64_C(0x6BC75E2D63100000)},
45+
detail::uint128 {UINT64_C(0x36), UINT64_C(0x35C9ADC5DEA00000)},
46+
detail::uint128 {UINT64_C(0x21E), UINT64_C(0x19E0C9BAB2400000)},
47+
detail::uint128 {UINT64_C(0x152D), UINT64_C(0x02C7E14AF6800000)},
48+
detail::uint128 {UINT64_C(0x84595), UINT64_C(0x161401484A000000)},
49+
detail::uint128 {UINT64_C(0x52B7D2), UINT64_C(0xDCC80CD2E4000000)},
50+
detail::uint128 {UINT64_C(0x33B2E3C), UINT64_C(0x9FD0803CE8000000)},
51+
detail::uint128 {UINT64_C(0x204FCE5E), UINT64_C(0x3E25026110000000)},
52+
detail::uint128 {UINT64_C(0x1431E0FAE), UINT64_C(0x6D7217CAA0000000)},
53+
detail::uint128 {UINT64_C(0xC9F2C9CD0), UINT64_C(0x4674EDEA40000000)},
54+
detail::uint128 {UINT64_C(0x7E37BE2022), UINT64_C(0xC0914B2680000000)},
55+
detail::uint128 {UINT64_C(0x4EE2D6D415B), UINT64_C(0x85ACEF8100000000)},
56+
detail::uint128 {UINT64_C(0x314DC6448D93), UINT64_C(0x38C15B0A00000000)},
57+
detail::uint128 {UINT64_C(0x1ED09BEAD87C0), UINT64_C(0x378D8E6400000000)},
58+
detail::uint128 {UINT64_C(0x13426172C74D82), UINT64_C(0x2B878FE800000000)},
59+
detail::uint128 {UINT64_C(0xC097CE7BC90715), UINT64_C(0xB34B9F1000000000)},
60+
detail::uint128 {UINT64_C(0x785EE10D5DA46D9), UINT64_C(0x00F436A000000000)},
61+
detail::uint128 {UINT64_C(0x4B3B4CA85A86C47A), UINT64_C(0x098A224000000000)}
62+
}};
63+
64+
for (const auto& val : powers_of_10)
65+
{
66+
const auto temp {boost::decimal::detail::remove_trailing_zeros(val)};
67+
if (!BOOST_TEST_EQ(temp.trimmed_number, detail::uint128(1)))
68+
{
69+
// LCOV_EXCL_START
70+
std::cerr << "Input Number: " << val
71+
<< "\nOutput Number: " << temp.trimmed_number
72+
<< "\nZeros removed: " << temp.number_of_removed_zeros << std::endl;
73+
// LCOV_EXCL_STOP
74+
}
75+
}
76+
}
77+
78+
int main()
79+
{
80+
test<std::uint32_t>();
81+
test<std::uint64_t>();
82+
test<boost::decimal::detail::uint128>();
83+
84+
test_extended();
85+
86+
return boost::report_errors();
87+
}

‎tools/granlund-montgomery.py‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
def extended_euclidean(a, b):
2+
if a == 0:
3+
return b, 0, 1
4+
gcd, x1, y1 = extended_euclidean(b % a, a)
5+
x = y1 - (b // a) * x1
6+
y = x1
7+
return gcd, x, y
8+
9+
def mod_inverse(a, m):
10+
gcd, x, _ = extended_euclidean(a, m)
11+
if gcd != 1:
12+
return 0 # Modular inverse doesn't exist
13+
else:
14+
return x % m
15+
16+
# Constants
17+
bits = int(128)
18+
t = int(32)
19+
20+
q = int(10)**t
21+
q0 = int(q / int(2)**t)
22+
print("Q0: ", q0)
23+
twobt_min_t = int(2**(bits - t))
24+
25+
# Calculate the modular inverse
26+
m0 = int(mod_inverse(q0, twobt_min_t))
27+
print("M0: ", m0)
28+
29+
p0 = int((q0 * m0 - 1) / twobt_min_t)
30+
print("P0: ", p0)
31+
32+
p = int(q0 + p0)
33+
print("P: ", p)
34+
35+
m = int((twobt_min_t * p + 1) / q0)
36+
print("M: ", m)
37+
38+
threshold_value = int(2**bits / q + 1)
39+
print("Threshold Value: ", threshold_value)

0 commit comments

Comments
 (0)