素因数分解・素数判定・約数列挙

math/factorize.hpp · GitHub で見る

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

#include "lib/math/factorize.hpp"

auto factors = factorize(360); // {2,2,2,3,3,5}
auto counts = factorize_count(360); // { {2,3},{3,2},{5,1} }
bool prime = is_prime(998244353); // true
auto ds = divisors(12); // {1,2,3,4,6,12}
auto unsorted = divisors(12, false); // 順序を指定せず列挙
auto reused = divisors(counts); // 360 の約数。再度の素因数分解は不要

操作	内容・返り値	計算量
is_prime(n)	素数かを判定	\(O(\log n)\)
factorize(n)	素因数を重複込み・昇順で返す。n >= 1	期待 \(O(n^{1/4}\log n)\)
factorize_count(n)	{素因数, 指数} の昇順列を返す。n >= 1	期待 \(O(n^{1/4}\log n)\)
divisors(n, sorted = true)	正の約数を vector<uint64_t> で返す。n >= 1	素因数分解に加え、昇順は \(O(d\log d)\)、false は \(O(d)\)
divisors(factors, sorted = true)	factorize_count 形式の分解結果から正の約数を返す	昇順は \(O(d\log d)\)、false は \(O(d)\)

引数は uint64_t です。factorize(n) は素因数の昇順列、factorize_count(n) は {素因数, 指数} の昇順列を返します。素因数分解は \(n \ge 1\) とします。

例えば counts の {2,3} は素因数 2 の指数が 3、すなわち \(2^3\) を表します。1 の分解結果はどちらも空配列です。0 と負数の素因数分解は扱いません。前処理を呼ぶ必要はなく、そのまま各関数を使えます。

divisors(n) は 1 と n を含む全ての正の約数を、重複なしの昇順で返します。divisors(1) は [1] です。第 2 引数に false を渡すとソートを省略します。この場合も全ての約数を 1 回ずつ返しますが、順序は保証しません。0 と負数の約数列挙は扱いません。

divisors(factors) の factors は vector<pair<uint64_t,int>> です。異なる素数と正の指数を渡し、その積は uint64_t に収まる必要があります。素数の並び順は問いません。空配列は 1 の分解結果として扱い、[1] を返します。渡した分解結果の素数判定や正当性の検査は行いません。

約数と同時にメビウス関数・オイラー関数などの値も列挙する場合は、約数ごとの乗法的関数の値 を利用できます。

64 bit 整数演算を \(O(1)\) として、is_prime(n) は \(O(\log n)\) 時間、\(O(1)\) 空間です。factorize(n)・factorize_count(n) は Pollard Rho の期待計算量で \(O(n^{1/4} \log n)\) 時間、\(O(\log n)\) 空間です。

約数の個数を \(d\) とすると、divisors(n) は期待 \(O(n^{1/4}\log n+d\log d)\) 時間、divisors(n,false) は期待 \(O(n^{1/4}\log n+d)\) 時間、いずれも返り値を含めて \(O(d+\log n)\) 空間です。分解結果を渡す版は素因数分解を省き、\(O(d)\) 空間です。

依存ライブラリ

• 素数判定の基数表
• 素数判定の共通演算
• Montgomery 64
• Montgomery 32
• SplitMix64
• 素数判定（高速版）

検証コード

• test/discrete_log.test.cpp
• test/factorize.test.cpp
• test/multiplicative_convolution.test.cpp

コード

#ifndef FACTORIZE_HPP
#define FACTORIZE_HPP

#include "is_prime.hpp"
#include "../components/splitmix64.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <limits>
#include <numeric>
#include <utility>
#include <vector>

namespace factorize_detail {
using U = std::uint64_t;
using W = __uint128_t;

struct Factorizer {
	SplitMix64 rng;
	std::vector<U> result;
	explicit Factorizer(U n) : rng(n ^ 0x9e3779b97f4a7c15ULL) { result.reserve(64); }
	U random() { return rng(); }

	struct RhoArithmetic {
		U mod, inverse, one;
		explicit RhoArithmetic(const Montgomery64 &m) : mod(m.mod), inverse(m.inverse), one(m.one) {}
		U mul(U a, U b) const {
			W x = W(a) * b;
			return U(x >> 64) + mod - U((W(U(x) * inverse) * mod) >> 64);
		}
		U step(U x, U c) const {
			return mul(x, x) + c;
		}
	};

	U rho_parallel(const Montgomery64 &original) {
		RhoArithmetic m(original);
		const U n = m.mod;
		bool first = true;
		for (;;) {
			U y1 = first ? 1 : random() % n, y2 = first ? 2 : random() % n;
			const U c1 = first ? 1 : (random() & 31) + 1, c2 = first ? 2 : (random() & 31) + 1;
			first = false;
			bool restart = false;
			for (int length = 64; length <= (1 << 18) && !restart; length <<= 1) {
				const U x1 = y1 + 2 * n, x2 = y2 + 2 * n;
				for (int k = 0; k < length; k += 1024) {
					const U saved1 = y1, saved2 = y2;
					U q1 = m.one, q2 = m.one;
					const int count = std::min(1024, length - k);
					y1 = m.step(y1, c1), y2 = m.step(y2, c2);
					for (int i = 0; i < count; ++i) {
						const U d1 = x1 - y1, d2 = x2 - y2;
						y1 = m.step(y1, c1), y2 = m.step(y2, c2);
						q1 = m.mul(q1, d1), q2 = m.mul(q2, d2);
					}
					q1 = m.mul(q1, x1 - y1), q2 = m.mul(q2, x2 - y2);
					const U g = std::gcd(m.mul(q1, q2), n);
					if (g == 1) continue;
					if (g != n) return g;
					const U g1 = std::gcd(q1, n), g2 = std::gcd(q2, n);
					if (g1 > 1 && g1 < n) return g1;
					if (g2 > 1 && g2 < n) return g2;
					U y = g1 == n ? saved1 : saved2;
					const U x = g1 == n ? x1 : x2, c = g1 == n ? c1 : c2;
					for (int i = 0; i <= count; ++i) {
						y = m.step(y, c);
						const U d = std::gcd(x - y, n);
						if (d > 1 && d < n) return d;
						if (d == n) break;
					}
					restart = true;
					break;
				}
			}
		}
	}
	U rho(const Montgomery64 &m) {
		if ((m.mod >> 32) && m.mod < (U(1) << 61)) return rho_parallel(m);
		U n = m.mod;
		for (;;) {
			U y = m.in(random() % (n - 1) + 1), c = m.in(random() % (n - 1) + 1);
			U x = 0, saved = 0, g = 1;
			auto step = [&](U v) { return m.add(m.mul(v, v), c); };
			auto diff = [](U a, U b) { return a < b ? b - a : a - b; };
			for (int length = 1; g == 1 && length <= (1 << 18); length <<= 1) {
				x = y;
				for (int repeat_index = 0; repeat_index < length; ++repeat_index) y = step(y);
				for (int k = 0; k < length && g == 1; k += 128) {
					saved = y;
					U q = m.one;
					for (int repeat_index = 0; repeat_index < std::min(128, length - k); ++repeat_index) {
						y = step(y);
						q = m.mul(q, diff(x, y));
					}
					g = std::gcd(q, n);
				}
			}
			if (g == n)
				for (int repeat_index = 0; repeat_index < 128; ++repeat_index) {
					saved = step(saved);
					g = std::gcd(diff(x, saved), n);
					if (g != 1) break;
				}
			if (g > 1 && g < n) return g;
		}
	}
	void split(U n) {
		if (n < 1681) {
			result.push_back(n);
			return;
		}
		if ((n & 7) == 1) {
			U r = std::sqrt(static_cast<long double>(n));
			r -= W(r) * r > n;
			if (W(r) * r == n) {
				int begin = static_cast<int>(result.size());
				split(r);
				int end = static_cast<int>(result.size());
				for (int i = begin; i < end; ++i) result.push_back(result[i]);
				return;
			}
		}
		Montgomery64 m(n);
		if (is_prime_detail::fast_prime(m)) result.push_back(n);
		else {
			U d = rho(m);
			split(d);
			split(n / d);
		}
	}
};
}

inline std::vector<std::uint64_t> factorize(std::uint64_t n) {
	assert(n > 0);
	if (n < 2) return {};
	factorize_detail::Factorizer f(n);
	for (auto p : is_prime_detail::small_primes)
		while (n % p == 0) f.result.push_back(p), n /= p;
	if (n > 1) f.split(n);
	std::sort(f.result.begin(), f.result.end());
	return std::move(f.result);
}

inline std::vector<std::pair<std::uint64_t, int>> factorize_count(std::uint64_t n) {
	std::vector<std::pair<std::uint64_t, int>> result;
	for (auto p : factorize(n)) {
		if (result.empty() || result.back().first != p) result.emplace_back(p, 1);
		else ++result.back().second;
	}
	return result;
}

inline std::vector<std::uint64_t> divisors(const std::vector<std::pair<std::uint64_t, int>> &factors, bool sorted = true) {
	std::size_t count = 1;
	for (auto [p, e] : factors) {
		assert(p >= 2 && e > 0);
		count *= std::size_t(e) + 1;
	}
	std::vector<std::uint64_t> result{1};
	result.reserve(count);
	for (auto [p, e] : factors) {
		std::size_t size = result.size();
		std::uint64_t power = 1;
		for (int i = 0; i < e; ++i) {
			assert(power <= std::numeric_limits<std::uint64_t>::max() / p);
			power *= p;
			for (std::size_t j = 0; j < size; ++j) {
				assert(result[j] <= std::numeric_limits<std::uint64_t>::max() / power);
				result.push_back(result[j] * power);
			}
		}
	}
	if (sorted) std::sort(result.begin(), result.end());
	return result;
}

inline std::vector<std::uint64_t> divisors(std::uint64_t n, bool sorted = true) {
	return divisors(factorize_count(n), sorted);
}

#endif

Back to top page