離散対数

math/discrete_log.hpp · GitHub で見る

共通ページ · 軽量実装なし

#include "lib/math/discrete_log.hpp"
#include <iostream>

auto x = discrete_log(2, 8, 13);
if (x) std::cout << *x << '\n'; // 3
else std::cout << -1 << '\n';

関数	内容・返り値	計算量
discrete_log(a,b,mod)	\(a^x\equiv b\pmod{\mathrm{mod}}\) を満たす最小の非負整数を optional<uint32_t> で返す。解がなければ nullopt	期待 \(\displaystyle O(\sqrt m+\log^2 m)\)

a と b は uint64_t、mod は uint32_t です。\(1\leq\mathrm{mod}<2^{32}\) に対応し、合成数や偶数の法、底と法が互いに素でない場合も扱えます。a と b は法以上でも構いません。負数は呼び出し前に非負の剰余へ正規化してください。

\(0^0=1\) とします。法が 1 の場合は常に 0 を返します。

例では \(2^3\equiv8\pmod{13}\) なので 3 を返します。戻り値は std::optional です。if (x) で解の有無を確認し、解があれば *x を使います。答えが 0 の場合も if (x) は真になります。

最大公約数で非周期部分を取り除き、底の位数を求めて Pohlig–Hellman 法を使います。素数位数の部分問題には baby-step giant-step 法を使い、中国剰余定理で指数を復元します。短い周期や小さい指数は先に確認します。固定した素数 mod と原始根に対する大量のクエリには、前計算型の ModPow も使えます。

法を \(m\) として、期待 \(\displaystyle O(\sqrt m+\log^2 m)\) 時間、追加空間 \(\displaystyle O(\sqrt m)\) です。位数が小さい素因数に分かれる場合は、法全体を探索する BSGS より探索量を減らせます。ハッシュ表には線形探索を使います。

依存ライブラリ

• 32 bit Barrett 演算
• 素数判定の基数表
• 素数判定の共通演算
• Montgomery 64
• Montgomery 32
• SplitMix64
• 素因数分解・素数判定・約数列挙
• 素数判定（高速版）

検証コード

• test/discrete_log.test.cpp

コード

#ifndef DISCRETE_LOG_HPP
#define DISCRETE_LOG_HPP

#include "../components/barrett32.hpp"
#include "factorize.hpp"
#include <bit>
#include <cmath>
#include <numeric>
#include <optional>
#include <vector>

namespace discrete_log_detail {
inline std::uint32_t inverse(std::uint32_t x, std::uint32_t mod) {
	std::int64_t u = 0, v = 1, r = mod, s = x;
	while (s) {
		std::int64_t q = r / s, t = r - q * s;
		r = s, s = t, t = u - q * v, u = v, v = t;
	}
	return std::uint32_t(u < 0 ? u + mod : u);
}

inline std::optional<std::uint32_t> search(std::uint32_t a, std::uint32_t b, const Barrett32 &m, std::uint32_t order) {
	std::uint32_t value = 1;
	for (std::uint32_t j = 0; j < 32; ++j, value = m.mul(value, a)) {
		if (value == b) return j;
		if (j && value == 1) return std::nullopt;
	}
	std::uint32_t step = std::sqrt(double(order));
	step += std::uint64_t(step) * step < order;
	std::uint32_t capacity = std::bit_ceil(2 * step), mask = capacity - 1;
	int shift = 64 - std::countr_zero(capacity);
	std::vector<std::uint64_t> table(capacity);
	auto hash = [&](std::uint32_t x) { return std::uint32_t((x * 11995408973635179863ULL) >> shift); };
	value = 1;
	for (std::uint32_t j = 0; j < step; ++j, value = m.mul(value, a)) {
		if (value == b) return j;
		if (j && value == 1) return std::nullopt;
		auto slot = hash(value);
		while (table[slot]) slot = (slot + 1) & mask;
		table[slot] = (std::uint64_t(value) << 32) | j;
	}
	std::uint32_t stride = inverse(value, m.mod), target = b;
	for (std::uint64_t i = step; i < order; i += step) {
		target = m.mul(target, stride);
		auto slot = hash(target);
		while (table[slot] && std::uint32_t(table[slot] >> 32) != target) slot = (slot + 1) & mask;
		if (table[slot]) {
			std::uint64_t answer = i + std::uint32_t(table[slot]);
			if (answer < order) return std::uint32_t(answer);
			break;
		}
	}
	return std::nullopt;
}
}

inline std::optional<std::uint32_t> discrete_log(std::uint64_t a, std::uint64_t b, std::uint32_t mod) {
	assert(mod > 0);
	a %= mod, b %= mod;
	std::uint32_t offset = 0, start = 1 % mod;
	for (;;) {
		if (b == start) return offset;
		std::uint32_t g = std::gcd(std::uint32_t(a), mod);
		if (g == 1) break;
		if (b % g) return std::nullopt;
		b /= g, mod /= g, ++offset;
		start = std::uint64_t(start) * (a / g) % mod;
	}
	if (std::gcd(std::uint32_t(b), mod) != 1) return std::nullopt;
	Barrett32 m(mod);
	a %= mod;
	b = m.mul(b, discrete_log_detail::inverse(start, mod));
	std::uint32_t value = 1;
	for (std::uint32_t j = 0; j < 32; ++j, value = m.mul(value, a)) {
		if (value == b) return offset + j;
		if (j && value == 1) return std::nullopt;
	}
	std::uint32_t order = mod;
	for (auto [p, e] : factorize_count(mod)) order -= order / p;
	auto factors = factorize_count(order);
	for (auto &[p, e] : factors) while (e && m.pow(a, order / p) == 1) order /= p, --e;
	if (m.pow(b, order) != 1) return std::nullopt;
	std::uint64_t answer = 0, period = 1;
	std::uint32_t inverse_a = m.pow(a, order - 1);
	for (auto [p, e] : factors) if (e) {
		std::uint32_t base = m.pow(a, order / p), target = b, residue = 0, power = 1, exponent = order;
		for (int i = 0; i < e; ++i) {
			exponent /= p;
			auto digit = discrete_log_detail::search(base, m.pow(target, exponent), m, p);
			if (!digit) return std::nullopt;
			std::uint32_t term = *digit * power;
			residue += term, power *= p;
			target = m.mul(target, m.pow(inverse_a, term));
		}
		std::uint32_t delta = (std::uint64_t(residue) + power - answer % power) % power;
		delta = std::uint64_t(delta) * discrete_log_detail::inverse(period % power, power) % power;
		answer += period * delta, period *= power;
	}
	return std::uint32_t(answer + offset);
}

#endif

Back to top page