ポラードのロー法（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個