discrete_log.hpp¶
用 exBSGS 求最小非负 \(x\) 使 \(a^x\equiv b\pmod m\)，允许底数与模数不互质，并在无解时返回空值。
\[
\displaystyle a^x\equiv b\pmod m,\qquad x=\min\{t\ge 0:a^t\equiv b\pmod m\}
\]
Complexity: Time: \(O(\sqrt{m})\) expected time. Space: \(O(\sqrt{m})\).
AC 记录：discrete_logarithm_mod。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(sqrt(m)) expected time.
/// Space: O(sqrt(m)).
#include <cassert>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <optional>
#include <unordered_map>
namespace noya {
namespace discrete_log_internal {
using u64 = std::uint64_t;
using i128 = __int128;
inline u64 multiply_mod(u64 a, u64 b, u64 mod) {
return u64(static_cast<unsigned __int128>(a) * b % mod);
}
inline u64 inverse_mod(u64 val, u64 mod) {
assert(mod >= 2 && std::gcd(val, mod) == 1);
i128 c0 = 1;
i128 cf = 0;
u64 r0 = val;
u64 rem = mod;
while (rem != 0) {
u64 quo = r0 / rem;
u64 r1 = r0 - quo * rem;
r0 = rem;
rem = r1;
i128 c1 = c0 - i128(quo) * cf;
c0 = cf;
cf = c1;
}
c0 %= i128(mod);
if (c0 < 0) {
c0 += mod;
}
return u64(c0);
}
inline u64 ceil_sqrt(u64 val) {
u64 res = u64(std::sqrt(static_cast<long double>(val)));
while (res != 0 && res > val / res) {
res--;
}
while (res == 0 || res <= val / res) {
if (res * res == val) {
return res;
}
res++;
}
return res;
}
inline std::optional<u64> coprime_bsgs(u64 bas, u64 tar, u64 mod) {
assert(mod >= 2 && std::gcd(bas, mod) == 1);
// Powers of a unit cannot reach a nonunit.
if (std::gcd(tar, mod) != 1) {
return std::nullopt;
}
// Both operands here are reduced; through 2^32 their product fits u64.
auto mul = [mod](u64 a, u64 b) {
return mod <= (u64(1) << 32) ? a * b % mod : multiply_mod(a, b, mod);
};
u64 blk = ceil_sqrt(mod);
std::unordered_map<u64, u64> bs;
u64 val = 1;
for (u64 exp = 0; exp < blk; exp++) {
if (val == tar) {
return exp;
}
if (exp != 0 && val == 1) {
return std::nullopt;
}
// Keep short answers/cycles small; reserve the full table only on demand.
if (exp == 64) {
bs.reserve(std::size_t(blk * 2 + 1));
}
bs.emplace(val, exp);
val = mul(val, bas);
}
if (val == 1) {
return std::nullopt;
}
// val is already bas^blk; invert it instead of recomputing this power.
u64 ib = inverse_mod(val, mod);
val = tar;
for (u64 gia = 0; gia <= blk; gia++) {
auto it = bs.find(val);
if (it != bs.end()) {
return gia * blk + it->second;
}
val = mul(val, ib);
}
return std::nullopt;
}
} // namespace discrete_log_internal
/// @brief Find the smallest x >= 0 satisfying bas^x = tar (mod mod),
/// or nullopt if no solution exists. Requires mod >= 1; bas and tar are
/// reduced internally. Full BSGS can require O(sqrt(mod)) hash entries.
inline std::optional<std::uint64_t>
discrete_log(std::uint64_t bas, std::uint64_t tar, std::uint64_t mod) {
using discrete_log_internal::coprime_bsgs;
using discrete_log_internal::multiply_mod;
assert(mod >= 1);
if (mod == 1) {
return 0;
}
bas %= mod;
tar %= mod;
if (tar == 1) {
return 0;
}
std::uint64_t rmv = 0;
std::uint64_t acc = 1;
while (true) {
std::uint64_t div = std::gcd(bas, mod);
if (div == 1) {
break;
}
if (tar == acc) {
return rmv;
}
if (tar % div != 0) {
return std::nullopt;
}
tar /= div;
mod /= div;
acc = multiply_mod(acc, bas / div, mod);
rmv++;
}
if (tar == acc) {
return rmv;
}
std::uint64_t nrm =
multiply_mod(tar, discrete_log_internal::inverse_mod(acc, mod), mod);
auto res = coprime_bsgs(bas % mod, nrm, mod);
if (!res) {
return std::nullopt;
}
return *res + rmv;
}
/// @brief Extended BSGS: find the smallest x >= 0 satisfying
/// bas^x = tar (mod mod), without requiring gcd(bas, mod) = 1.
inline std::optional<std::uint64_t> exbsgs(std::uint64_t bas, std::uint64_t tar,
std::uint64_t mod) {
return discrete_log(bas, tar, mod);
}
} // namespace noya
// NOYA_DISCRETE_LOG_HPP
#ifndef NOYA_DISCRETE_LOG_HPP
#define NOYA_DISCRETE_LOG_HPP 1
/// @complexity Time: O(sqrt(m)) expected time.
/// Space: O(sqrt(m)).
#include <cassert>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <optional>
#include <unordered_map>
namespace noya {
namespace discrete_log_internal {
using u64 = std::uint64_t;
using i128 = __int128;
inline u64 multiply_mod(u64 a, u64 b, u64 mod) {
return u64(static_cast<unsigned __int128>(a) * b % mod);
}
inline u64 inverse_mod(u64 val, u64 mod) {
assert(mod >= 2 && std::gcd(val, mod) == 1);
i128 c0 = 1;
i128 cf = 0;
u64 r0 = val;
u64 rem = mod;
while (rem != 0) {
u64 quo = r0 / rem;
u64 r1 = r0 - quo * rem;
r0 = rem;
rem = r1;
i128 c1 = c0 - i128(quo) * cf;
c0 = cf;
cf = c1;
}
c0 %= i128(mod);
if (c0 < 0) {
c0 += mod;
}
return u64(c0);
}
inline u64 ceil_sqrt(u64 val) {
u64 res = u64(std::sqrt(static_cast<long double>(val)));
while (res != 0 && res > val / res) {
res--;
}
while (res == 0 || res <= val / res) {
if (res * res == val) {
return res;
}
res++;
}
return res;
}
inline std::optional<u64> coprime_bsgs(u64 bas, u64 tar, u64 mod) {
assert(mod >= 2 && std::gcd(bas, mod) == 1);
// Powers of a unit cannot reach a nonunit.
if (std::gcd(tar, mod) != 1) {
return std::nullopt;
}
// Both operands here are reduced; through 2^32 their product fits u64.
auto mul = [mod](u64 a, u64 b) {
return mod <= (u64(1) << 32) ? a * b % mod : multiply_mod(a, b, mod);
};
u64 blk = ceil_sqrt(mod);
std::unordered_map<u64, u64> bs;
u64 val = 1;
for (u64 exp = 0; exp < blk; exp++) {
if (val == tar) {
return exp;
}
if (exp != 0 && val == 1) {
return std::nullopt;
}
// Keep short answers/cycles small; reserve the full table only on demand.
if (exp == 64) {
bs.reserve(std::size_t(blk * 2 + 1));
}
bs.emplace(val, exp);
val = mul(val, bas);
}
if (val == 1) {
return std::nullopt;
}
// val is already bas^blk; invert it instead of recomputing this power.
u64 ib = inverse_mod(val, mod);
val = tar;
for (u64 gia = 0; gia <= blk; gia++) {
auto it = bs.find(val);
if (it != bs.end()) {
return gia * blk + it->second;
}
val = mul(val, ib);
}
return std::nullopt;
}
} // namespace discrete_log_internal
/// @brief Find the smallest x >= 0 satisfying bas^x = tar (mod mod),
/// or nullopt if no solution exists. Requires mod >= 1; bas and tar are
/// reduced internally. Full BSGS can require O(sqrt(mod)) hash entries.
inline std::optional<std::uint64_t>
discrete_log(std::uint64_t bas, std::uint64_t tar, std::uint64_t mod) {
using discrete_log_internal::coprime_bsgs;
using discrete_log_internal::multiply_mod;
assert(mod >= 1);
if (mod == 1) {
return 0;
}
bas %= mod;
tar %= mod;
if (tar == 1) {
return 0;
}
std::uint64_t rmv = 0;
std::uint64_t acc = 1;
while (true) {
std::uint64_t div = std::gcd(bas, mod);
if (div == 1) {
break;
}
if (tar == acc) {
return rmv;
}
if (tar % div != 0) {
return std::nullopt;
}
tar /= div;
mod /= div;
acc = multiply_mod(acc, bas / div, mod);
rmv++;
}
if (tar == acc) {
return rmv;
}
std::uint64_t nrm =
multiply_mod(tar, discrete_log_internal::inverse_mod(acc, mod), mod);
auto res = coprime_bsgs(bas % mod, nrm, mod);
if (!res) {
return std::nullopt;
}
return *res + rmv;
}
/// @brief Extended BSGS: find the smallest x >= 0 satisfying
/// bas^x = tar (mod mod), without requiring gcd(bas, mod) = 1.
inline std::optional<std::uint64_t> exbsgs(std::uint64_t bas, std::uint64_t tar,
std::uint64_t mod) {
return discrete_log(bas, tar, mod);
}
} // namespace noya
#endif // NOYA_DISCRETE_LOG_HPP
/// @complexity Time: O(sqrt(m)) expected time.
/// Space: O(sqrt(m)).
#include <cassert>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <optional>
#include <unordered_map>
namespace noya {
namespace discrete_log_internal {
using u64 = std::uint64_t;
using i128 = __int128;
inline u64 multiply_mod(u64 a, u64 b, u64 mod) {
return u64(static_cast<unsigned __int128>(a) * b % mod);
}
inline u64 inverse_mod(u64 val, u64 mod) {
assert(mod >= 2 && std::gcd(val, mod) == 1);
i128 c0 = 1;
i128 cf = 0;
u64 r0 = val;
u64 rem = mod;
while (rem != 0) {
u64 quo = r0 / rem;
u64 r1 = r0 - quo * rem;
r0 = rem;
rem = r1;
i128 c1 = c0 - i128(quo) * cf;
c0 = cf;
cf = c1;
}
c0 %= i128(mod);
if (c0 < 0) {
c0 += mod;
}
return u64(c0);
}
inline u64 ceil_sqrt(u64 val) {
u64 res = u64(std::sqrt(static_cast<long double>(val)));
while (res != 0 && res > val / res) {
res--;
}
while (res == 0 || res <= val / res) {
if (res * res == val) {
return res;
}
res++;
}
return res;
}
inline std::optional<u64> coprime_bsgs(u64 bas, u64 tar, u64 mod) {
assert(mod >= 2 && std::gcd(bas, mod) == 1);
// Powers of a unit cannot reach a nonunit.
if (std::gcd(tar, mod) != 1) {
return std::nullopt;
}
// Both operands here are reduced; through 2^32 their product fits u64.
auto mul = [mod](u64 a, u64 b) {
return mod <= (u64(1) << 32) ? a * b % mod : multiply_mod(a, b, mod);
};
u64 blk = ceil_sqrt(mod);
std::unordered_map<u64, u64> bs;
u64 val = 1;
for (u64 exp = 0; exp < blk; exp++) {
if (val == tar) {
return exp;
}
if (exp != 0 && val == 1) {
return std::nullopt;
}
// Keep short answers/cycles small; reserve the full table only on demand.
if (exp == 64) {
bs.reserve(std::size_t(blk * 2 + 1));
}
bs.emplace(val, exp);
val = mul(val, bas);
}
if (val == 1) {
return std::nullopt;
}
// val is already bas^blk; invert it instead of recomputing this power.
u64 ib = inverse_mod(val, mod);
val = tar;
for (u64 gia = 0; gia <= blk; gia++) {
auto it = bs.find(val);
if (it != bs.end()) {
return gia * blk + it->second;
}
val = mul(val, ib);
}
return std::nullopt;
}
} // namespace discrete_log_internal
/// @brief Find the smallest x >= 0 satisfying bas^x = tar (mod mod),
/// or nullopt if no solution exists. Requires mod >= 1; bas and tar are
/// reduced internally. Full BSGS can require O(sqrt(mod)) hash entries.
inline std::optional<std::uint64_t>
discrete_log(std::uint64_t bas, std::uint64_t tar, std::uint64_t mod) {
using discrete_log_internal::coprime_bsgs;
using discrete_log_internal::multiply_mod;
assert(mod >= 1);
if (mod == 1) {
return 0;
}
bas %= mod;
tar %= mod;
if (tar == 1) {
return 0;
}
std::uint64_t rmv = 0;
std::uint64_t acc = 1;
while (true) {
std::uint64_t div = std::gcd(bas, mod);
if (div == 1) {
break;
}
if (tar == acc) {
return rmv;
}
if (tar % div != 0) {
return std::nullopt;
}
tar /= div;
mod /= div;
acc = multiply_mod(acc, bas / div, mod);
rmv++;
}
if (tar == acc) {
return rmv;
}
std::uint64_t nrm =
multiply_mod(tar, discrete_log_internal::inverse_mod(acc, mod), mod);
auto res = coprime_bsgs(bas % mod, nrm, mod);
if (!res) {
return std::nullopt;
}
return *res + rmv;
}
/// @brief Extended BSGS: find the smallest x >= 0 satisfying
/// bas^x = tar (mod mod), without requiring gcd(bas, mod) = 1.
inline std::optional<std::uint64_t> exbsgs(std::uint64_t bas, std::uint64_t tar,
std::uint64_t mod) {
return discrete_log(bas, tar, mod);
}
} // namespace noya
// NOYA_DISCRETE_LOG_HPP