int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 21, 23, 29, 31, 37}; IL boolMiller_Rabin(ll n){ if (n <= 1) returnfalse; for (int i = 0; i < 9; ++i) if (n % pri[i] == 0) return n == pri[i]; ll r, x, y; int t; for (r = n - 1, t = 0; ~r & 1; r >>= 1, t++); for (int i = 0; i < 9; ++i) { x = fpw(pri[i], r, n); for (int j = 0; j < t && x > 1; ++j) { y = mul(x, x, n); if (y == 1 && x != n - 1) returnfalse; x = y; } if (x != 1) returnfalse; } returntrue; }
IL ll func(ll x, ll mod, ll a){ return (mul(x, x, mod) + a) % mod; } IL ll Find(ll n){ staticconstint Test_Limit = 150000; ll a = rand(), x, y; int tim = 0; for (int i = 0; i < 12; ++i) if (n % pri[i] == 0) return pri[i]; x = func(rand(), n, a), y = x; do { ull g = gcd((x - y + n) % n, n); if (g != 1 && g != n) return g; x = func(x, n, a), y = func(func(y, n, a), n, a); ++tim; } while (y != x && tim <= Test_Limit); return-1; } voidPollard_s_Rho(ll n, ll &a){ ll d; if (n <= a) return; while (!(n & 1)) n >>= 1, a = 2; if (n == 1 || Miller_Rabin(n)) { a = max(a, n); return; } for (d = Find(n); d == -1; d = Find(n)); if (d > n / d) d = n / d; Pollard_s_Rho(d, a); Pollard_s_Rho(n / d, a); }
Pollard’s Rho Update
但这并不是能达到的最优复杂度,可以优化为。
首先的写法可以加上一点常数优化。
1 2 3 4 5 6 7 8 9 10
IL ll gcd(ll x, ll y){ int Shift = ctzll(x | y); y >>= ctzll(y); while (x) { x >>= ctzll(x); if (x < y) swap(x, y); x -= y; } return y << Shift; }
#define ctzll __builtin_ctzll namespace Rho { IL ll gcd(ll x, ll y){ int Shift = ctzll(x | y); y >>= ctzll(y); while (x) { x >>= ctzll(x); if (x < y) swap(x, y); x -= y; } return y << Shift; } IL ll func(ll x, ll Mod, ll c){ return (mul(x, x, Mod) + c) % Mod; } IL ll Find(ll n){ staticint Step = 1 << 9; int c = rand(); ll x, y, temp; x = y = rand() % n; for (int l = 1; ; l <<= 1) { x = y; for (int i = 0; i < l; ++i) y = func(y, n, c); for (int k = 0; k < l; k += Step) { ll Prod_g = 1; temp = y; int rem = min(l - k, Step); for (int i = 0; i < rem; ++i) y = func(y, n, c), Prod_g = mul(Prod_g, (y + n - x) % n, n); Prod_g = gcd(Prod_g, n); if (Prod_g == 1) continue; if (Prod_g == n) for (y = temp, Prod_g = 1; Prod_g == 1; ) y = func(y, n, c), Prod_g = gcd(y + n - x, n); return Prod_g; } } } IL voidPollard_s_Rho(ll n){} // 同上 }