ポラードのロー法(Pollard's rho algorithm)は、与えられた数n
の非自明な因数を一つ見つける乱択アルゴリズムです。この手法が答えを返した場合、解は常に正しいことが保証されます。ただし、このアルゴリズムは「失敗」を返して終了することもあります*1。
ポラードのロー法の概要
ポラードのロー法では、与えられた数n
の任意の因数p
について、乱数列のうちmod p
で見た時に同じ値になる部分を、フロイドの循環検出アルゴリズムを用いて発見します。
発見された二つの乱数値をx
、y
とすると、gcd(|x-y|, n)
はp
の倍数になるため、これを非自明な因数として返却します。
完全にランダムな乱数列を使った場合、mod p
で見た時に同じ値のペアを一つ見つけるのに必要な乱数生成回数は、回です。
この手法の肝は、n
の因数p
について、疑似乱数列の中にあるmod p
で見た時に同じ値になるペアを、p
を知らないままに検出することです。
そのためには、mod p
で循環する疑似乱数列を使用しないといけません。ここで、p
はまだわからない点が難しいですが、幸いmod n
での線形合同法が条件を満たします*2。
素因数分解プログラム
先週紹介した、高速gcd
プログラムを使用しています。
clang++
でコンパイルする場合に高速化するようなコードになっています。
__builtin_ctzll
や__uint128_t
など非標準の機能を使っていますが、勘弁してください。
#include <cstdint> #include <iostream> std::uint64_t gcd_impl( std::uint64_t x, std::uint64_t y ) noexcept { if( x == y ) { return x; } const std::uint64_t a = y - x; const std::uint64_t b = x - y; const int n = __builtin_ctzll( b ); const std::uint64_t s = x < y ? a : b; const std::uint64_t t = x < y ? x : y; return ::gcd_impl( s >> n, t ); } std::uint64_t gcd( std::uint64_t x, std::uint64_t y ) noexcept { const int n = __builtin_ctzll( x ); const int m = __builtin_ctzll( y ); return ::gcd_impl( x >> n, y >> m ) << ( n < m ? n : m ); } int main() { std::uint64_t n; std::cin >> n; for( std::uint64_t x = 2, y = 2; ; ) { x = ((__uint128_t)x * x + 1) % n; y = ((__uint128_t)y * y + 1) % n; y = ((__uint128_t)y * y + 1) % n; if( std::uint64_t z = ::gcd( x < y ? y - x : x - y, n ); z != 1 ) { std::cout << n << " = " << z << " * " << n / z << std::endl; return 0; } } }
このソースコードをclang++-10 -std=c++17 -O3
でコンパイルした場合、10023859281455311421*3の素因数分解を0.01秒で完了することができました。
ちなみに、21,25,95,125, ... などは、合成数にもかかわらず、このプログラムで素因数分解できません。 そういった数の存在数を以下の表にまとめてみました。
範囲 | 個数 |
---|---|
102以下 | 3個 |
103以下 | 14個 |
104以下 | 62個 |
105以下 | 213個 |
106以下 | 919個 |
107以下 | 3660個 |