ポラードのロー法で64bit符号なし整数を素因数分解する

ポラードのロー法(Pollard's rho algorithm)は、与えられた数nの非自明な因数を一つ見つける乱択アルゴリズムです。この手法が答えを返した場合、解は常に正しいことが保証されます。ただし、このアルゴリズムは「失敗」を返して終了することもあります*1

ポラードのロー法の概要

ポラードのロー法では、与えられた数nの任意の因数pについて、乱数列のうちmod pで見た時に同じ値になる部分を、フロイドの循環検出アルゴリズムを用いて発見します。 発見された二つの乱数値をxyとすると、gcd(|x-y|, n)pの倍数になるため、これを非自明な因数として返却します。 完全にランダムな乱数列を使った場合、mod pで見た時に同じ値のペアを一つ見つけるのに必要な乱数生成回数は、 \Theta(\sqrt 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個

*1:乱数生成器のパラメータ変更のみによって常に答えが得られるのかは知りません。

*2:線形合同法の乱数の質が条件を満たしているのかは怪しいところですが、実験的にはうまくいっているようです。

*3:Wikipediaに載っていた例で、1.1×263程度の大きさ