#include <cstdint>
#include <iostream>
std::uint64_t mul_mod(std::uint64_t a, std::uint64_t b,
std::uint64_t mod)
{
std::uint64_t res = 0;
while
(b > 0) {
if
(b & 1) {
res = (res + a) % mod;
}
a = (a * 2) % mod;
b >>= 1;
}
return
res;
}
std::uint64_t pow_mod(std::uint64_t a, std::uint64_t b,
std::uint64_t mod)
{
std::uint64_t res = 1;
a %= mod;
while
(b > 0) {
if
(b & 1) {
res = mul_mod(res, a, mod);
}
a = mul_mod(a, a, mod);
b >>= 1;
}
return
res;
}
std::uint64_t barrett_reduce(std::uint64_t x,
std::uint64_t mod,
std::uint64_t mu)
{
std::uint64_t q1 = x >> 32;
std::uint64_t q2 = ((mu * q1) >> 32) * mod;
std::uint64_t r1 = x & 0xffffffff;
std::uint64_t r2 = mu * r1;
std::uint64_t q3 = r2 >> 32;
std::uint64_t r3 = r1 - q3 * mod;
std::uint64_t res = r3 + ((r3 >> 63) & mod);
if
(res >= mod) {
res -= mod;
}
return
res;
}
int
main()
{
std::uint64_t x = 123456789;
std::uint64_t mod = 1000000007;
std::uint64_t mu = ((1ULL << 64) + mod - 1) / mod;
std::uint64_t x_squared = mul_mod(x, x, mod);
std::uint64_t x_squared_barrett
= barrett_reduce(x_squared, mod, mu);
std::cout <<
"x^2 % mod = "
<< x_squared_barrett
<< std::endl;
std::uint64_t x_pow_1234 = pow_mod(x, 1234, mod);
std::uint64_t x_pow_1234_barrett
= barrett_reduce(x_pow_1234, mod, mu);
std::cout <<
"x^1234 % mod = "
<< x_pow_1234_barrett
<< std::endl;
return
0;
}