Library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub ret2home/Library

:x: Fast Prime Factorization
(math/factor.cpp)

概要

ミラーラビン素数判定法と、ロー法による高速な素因数分解

計算量

詳しくは知らないけど、$10^{18}$ の素因数分解 $1000$ 個くらいは余裕ってレベル

Depends on

Verified with

Code

#pragma once
#include "../template/template.cpp"

ll pow128(ll x, ll y, ll m) {
    ll res = 1;
    while (y > 0) {
        if (y & 1) res = (__int128_t(res) * x) % m;
        y >>= 1;
        x = (__int128_t(x) * x) % m;
    }
    return res;
}

bool prime(ll n) {
    if (n < 2 || n % 2 == 0) return n == 2;
    ll d = n - 1;
    while (d % 2 == 0) d >>= 1;
    for (ll x : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
        if (n <= x) break;
        ll t = d, y = pow128(x, t, n);
        while (t != n - 1 && y != 1 && y != n - 1) {
            y = (__int128_t(y) * y) % n;
            t <<= 1;
        }
        if (y != n - 1 && t % 2 == 0) return false;
    }
    return true;
}
ll rho(ll n) {
    if (n % 2 == 0) return 2;
    ll MOD = n;
    auto f = [&](ll x) -> ll { return ((__int128_t)x * x + 1) % MOD; };
    auto g = [](ll x, ll y) -> ll {
        if (x < y) return y - x;
        return x - y;
    };
    for (ll x1 = 0;; x1++) {
        for (ll x = x1, y = f(x);; x = f(x), y = f(f(y))) {
            ll d = gcd(g(y, x), n);
            if (1 < d && d < n) return d;
            if (d == n) break;
        }
    }
}
vector<ll> factor(ll n) {
    if (n == 1) return {};
    if (prime(n)) return {n};
    ll d = rho(n);
    vector<ll> res = factor(d);
    for (ll i : factor(n / d)) res.push_back(i);
    return res;
}
/*
@brief Fast Prime Factorization
@docs docs/factor.md
*/
#line 2 "template/template.cpp"
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define rep(i, n) for (int i = 0; i < n; i++)
#define REP(i, n) for (int i = 1; i < n; i++)
#define rev(i, n) for (int i = n - 1; i >= 0; i--)
#define REV(i, n) for (int i = n - 1; i > 0; i--)
#define all(v) v.begin(), v.end()
#define PL pair<ll, ll>
#define PI pair<int, int>
#define pi acos(-1)
#define len(s) (int)s.size()
#define compress(v) \
    sort(all(v));   \
    v.erase(unique(all(v)), v.end());
#define comid(v, x) lower_bound(all(v), x) - v.begin()

template<class T>
using prique=priority_queue<T,vector<T>,greater<>>;

template <class T, class U>
inline bool chmin(T &a, U b) {
    if (a > b) {
        a = b;
        return true;
    }
    return false;
}
template <class T, class U>
inline bool chmax(T &a, U b) {
    if (a < b) {
        a = b;
        return true;
    }
    return false;
}
constexpr ll inf = 3e18;
#line 3 "math/factor.cpp"

ll pow128(ll x, ll y, ll m) {
    ll res = 1;
    while (y > 0) {
        if (y & 1) res = (__int128_t(res) * x) % m;
        y >>= 1;
        x = (__int128_t(x) * x) % m;
    }
    return res;
}

bool prime(ll n) {
    if (n < 2 || n % 2 == 0) return n == 2;
    ll d = n - 1;
    while (d % 2 == 0) d >>= 1;
    for (ll x : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
        if (n <= x) break;
        ll t = d, y = pow128(x, t, n);
        while (t != n - 1 && y != 1 && y != n - 1) {
            y = (__int128_t(y) * y) % n;
            t <<= 1;
        }
        if (y != n - 1 && t % 2 == 0) return false;
    }
    return true;
}
ll rho(ll n) {
    if (n % 2 == 0) return 2;
    ll MOD = n;
    auto f = [&](ll x) -> ll { return ((__int128_t)x * x + 1) % MOD; };
    auto g = [](ll x, ll y) -> ll {
        if (x < y) return y - x;
        return x - y;
    };
    for (ll x1 = 0;; x1++) {
        for (ll x = x1, y = f(x);; x = f(x), y = f(f(y))) {
            ll d = gcd(g(y, x), n);
            if (1 < d && d < n) return d;
            if (d == n) break;
        }
    }
}
vector<ll> factor(ll n) {
    if (n == 1) return {};
    if (prime(n)) return {n};
    ll d = rho(n);
    vector<ll> res = factor(d);
    for (ll i : factor(n / d)) res.push_back(i);
    return res;
}
/*
@brief Fast Prime Factorization
@docs docs/factor.md
*/
Back to top page