:heavy_check_mark: 高速素因数分解 (Math/fast_factorizer.hpp)

FastFactorizer

Pollard’s rho アルゴリズムを用いた素因数分解を行います。 素数判定にはMiller-Rabin素数判定法を使用しています。

FastFactorizer::factorize

static vector<long long> FastFactorizer::factorize(long long N)

整数 N を素因数分解し、因数を昇順に並べたリストを返します。

制約

  • $0 \lt N$

計算量

  • $O(\sqrt[4]{N})$

Depends on

Verified with

Code

#include <bits/stdc++.h>
using namespace std;

#include "Math/primality_test.hpp"

struct FastFactorizer
{
private:
    static unsigned long long abs_diff(unsigned long long a, unsigned long long b)
    {
        return (a > b) ? (a - b) : (b - a);
    }

    // pollards_rho法によりNの約数を返す
    static unsigned long long pollards_rho(unsigned long long N)
    {
        if (PrimalityTest::is_prime(N))
        {
            return N;
        }

        unsigned long long m = 1ll << (bit_width(N) / 8);

        // 100回試して見つからなかったら失敗
        for (unsigned long long c = 1; c < 100; c++)
        {
            auto f = [&](unsigned long long x) -> unsigned long long
            { return ((__uint128_t(x) * x) + c) % N; };

            unsigned long long x0 = c % (N + 1);
            unsigned long long y = x0, q = 1, g = 1, ys = 0, x = 0;

            for (int r = 1; g <= 1; r <<= 1)
            {
                x = y;
                for (int i = 0; i < r; i++)
                {
                    y = f(y);
                }

                unsigned long long k = 0;
                while (k < r && g <= 1)
                {
                    ys = y;
                    for (int i = 0; i < min(m, r - k); i++)
                    {
                        y = f(y);
                        q = __uint128_t(q) * abs_diff(x, y) % N;
                    }
                    g = gcd(q, N);
                    k += m;
                }
            }

            if (g == N)
            {
                while (g <= 1)
                {
                    ys = f(ys);
                    g = gcd(abs_diff(x, ys), N);
                }
            }

            if (g != N)
            {
                return g;
            }
        }

        assert(false && "fail to factorize by pollards rho algorithm");
    }

public:
    static vector<long long> factorize(long long N)
    {
        assert(0 < N);

        if (N == 1)
        {
            return {};
        }

        long long p = pollards_rho(N);
        if (p == N)
        {
            return {p};
        }

        vector<long long> left = factorize(p);
        vector<long long> right = factorize(N / p);
        vector<long long> result;
        merge(left.begin(), left.end(), right.begin(), right.end(), back_inserter(result));
        return result;
    }
};
#line 1 "Math/fast_factorizer.hpp"
#include <bits/stdc++.h>
using namespace std;

#line 2 "Math/primality_test.hpp"
using namespace std;

#line 2 "Math/mod_power.hpp"
using namespace std;

#line 2 "Other/fast_power.hpp"
using namespace std;

template <class S>
S fast_pow(S x, long long n, function<S(S, S)> mul, function<S()> e)
{
    assert(0 <= n);

    S ans = e();

    while (n)
    {
        if (n & 1)
        {
            ans = mul(ans, x);
        }
        x = mul(x, x);
        n >>= 1;
    }

    return ans;
}
#line 5 "Math/mod_power.hpp"

template <typename T>
enable_if_t<is_integral_v<T> || is_same_v<T, __int128_t>, T>
mod_pow(T x, T n, T mod)
{
    assert(0 <= n);
    assert(0 < mod);
    assert(mod <= numeric_limits<T>::max() / mod);

    x %= mod;
    if (x < 0)
    {
        x += mod;
    }

    auto mul = [&](T a, T b) -> T
    {
        return (a * b) % mod;
    };

    auto e = [&]() -> T
    {
        return 1;
    };

    return fast_pow<T>(x, n, mul, e);
}
#line 5 "Math/primality_test.hpp"

struct PrimalityTest
{
private:
    static bool miller_rabin(long long N, const vector<long long> &A)
    {
        const int s = countr_zero((unsigned long long)N - 1);
        const long long d = (N - 1) >> s;

        auto is_composite = [&](long long a) -> bool
        {
            // Fermat の小定理が成り立たない
            if (a % N == 0)
            {
                return false;
            }

            long long x = mod_pow<__int128_t>(a, d, N);
            if (x == 1 || x == N - 1)
            {
                return false;
            }

            for (int t = 0; t < s - 1; ++t)
            {
                x = __int128_t(x) * x % N;
                if (x == N - 1)
                {
                    return false;
                }
            }

            return true;
        };

        return none_of(A.begin(), A.end(), is_composite);
    }

public:
    static bool is_prime(long long N)
    {
        assert(0 <= N);
        if (N <= 2)
        {
            return N == 2;
        }
        if (N % 2 == 0)
        {
            return false;
        }

        return miller_rabin(N, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
    }
};
#line 5 "Math/fast_factorizer.hpp"

struct FastFactorizer
{
private:
    static unsigned long long abs_diff(unsigned long long a, unsigned long long b)
    {
        return (a > b) ? (a - b) : (b - a);
    }

    // pollards_rho法によりNの約数を返す
    static unsigned long long pollards_rho(unsigned long long N)
    {
        if (PrimalityTest::is_prime(N))
        {
            return N;
        }

        unsigned long long m = 1ll << (bit_width(N) / 8);

        // 100回試して見つからなかったら失敗
        for (unsigned long long c = 1; c < 100; c++)
        {
            auto f = [&](unsigned long long x) -> unsigned long long
            { return ((__uint128_t(x) * x) + c) % N; };

            unsigned long long x0 = c % (N + 1);
            unsigned long long y = x0, q = 1, g = 1, ys = 0, x = 0;

            for (int r = 1; g <= 1; r <<= 1)
            {
                x = y;
                for (int i = 0; i < r; i++)
                {
                    y = f(y);
                }

                unsigned long long k = 0;
                while (k < r && g <= 1)
                {
                    ys = y;
                    for (int i = 0; i < min(m, r - k); i++)
                    {
                        y = f(y);
                        q = __uint128_t(q) * abs_diff(x, y) % N;
                    }
                    g = gcd(q, N);
                    k += m;
                }
            }

            if (g == N)
            {
                while (g <= 1)
                {
                    ys = f(ys);
                    g = gcd(abs_diff(x, ys), N);
                }
            }

            if (g != N)
            {
                return g;
            }
        }

        assert(false && "fail to factorize by pollards rho algorithm");
    }

public:
    static vector<long long> factorize(long long N)
    {
        assert(0 < N);

        if (N == 1)
        {
            return {};
        }

        long long p = pollards_rho(N);
        if (p == N)
        {
            return {p};
        }

        vector<long long> left = factorize(p);
        vector<long long> right = factorize(N / p);
        vector<long long> result;
        merge(left.begin(), left.end(), right.begin(), right.end(), back_inserter(result));
        return result;
    }
};
Back to top page