結果
問題 |
No.36 素数が嫌い!
|
ユーザー |
|
提出日時 | 2025-05-22 18:59:39 |
言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
結果 |
AC
|
実行時間 | 2 ms / 5,000 ms |
コード長 | 5,403 bytes |
コンパイル時間 | 2,202 ms |
コンパイル使用メモリ | 207,988 KB |
実行使用メモリ | 7,844 KB |
最終ジャッジ日時 | 2025-05-22 18:59:43 |
合計ジャッジ時間 | 3,455 ms |
ジャッジサーバーID (参考情報) |
judge3 / judge1 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
sample | AC * 4 |
other | AC * 26 |
ソースコード
#include <bits/stdc++.h> using namespace std; struct Montgomery{ //2^62未満&奇数modのみ. //初めにsetmodする. using u64 = uint64_t; using u128 = __uint128_t; private: static u64 mod,N2,Rsq; //N*N2≡1(mod N); //Rsq = R^2modN; R=2^64. u64 v = 0; public: long long val(){return reduce(v);} u64 getmod(){return mod;} static void setmod(u64 m){ assert(m<(1LL<<62)&&(m&1)); mod = m; N2 = mod; for(int i=0; i<5; i++) N2 *= 2-N2*mod; Rsq = (-u128(mod))%mod; } //reduce = T*R^-1modNを求める. u64 reduce(const u128 &T){ //T*R^-1≡(T+(T*(-N2))modR*N)/R 2N未満なので-N必要かだけで良い. u64 ret = (T+u128(((u64)T)*(-N2))*mod)>>64; if(ret >= mod) ret -= mod; return ret; } //初期値<mod. 初めにw*R modN...->reduce(R^2)でok. Montgomery(){v = 0;} Montgomery(long long w):v(reduce(u128(w)*Rsq)){} Montgomery& operator=(const Montgomery &b) = default; Montgomery operator-()const{return Montgomery()-Montgomery(*this);} Montgomery operator+(const Montgomery &b)const{return Montgomery(*this)+=b;} Montgomery operator-(const Montgomery &b)const{return Montgomery(*this)-=b;} Montgomery operator*(const Montgomery &b)const{return Montgomery(*this)*=b;} Montgomery operator/(const Montgomery &b)const{return Montgomery(*this)/=b;} Montgomery& operator+=(const Montgomery &b){ v += b.v; if(v >= mod) v -= mod; return (*this); } Montgomery& operator-=(const Montgomery &b){ v += mod-b.v; if(v >= mod) v -= mod; return (*this); } Montgomery& operator*=(const Montgomery &b){ v = reduce(u128(v)*b.v); return (*this); } Montgomery& operator/=(const Montgomery &b){ (*this) *= b.inv(); return (*this); } Montgomery pow(u64 b)const{ Montgomery ret = 1,p = (*this); while(b){ if(b&1) ret *= p; p *= p; b >>= 1; } return ret; } Montgomery inv()const{return pow(mod-2);} bool operator!=(const Montgomery &b)const{return v!=b.v;} bool operator==(const Montgomery &b)const{return v==b.v;} }; typename Montgomery::u64 Montgomery::mod,Montgomery::N2,Montgomery::Rsq; using mont = Montgomery; bool MillerRabin(long long N,const vector<long long> &A){ mont::setmod(N); long long s = __builtin_ctzll(N-1),d = N-1; d >>= s; for(auto &a : A){ if(N <= a) break; mont x = mont(a).pow(d); if(x != 1){ long long t; for(t=0; t<s; t++){ if(x == N-1) break; x *= x; } if(t == s) return false; } } return true; } bool isprime(const long long N){ if(N <= 1) return false; else if(N == 2) return true; else if(N%2 == 0) return false; else if(N < 4759123141LL) return MillerRabin(N,{2,7,61}); else return MillerRabin(N, {2,325,9375,28178,450775,9780504,1795265022}); } long long stein_gcd(long long a,long long b){ if((!a)||(!b)) return a+b; int n = __builtin_ctzll(a); int m = __builtin_ctzll(b); auto f = [](auto f,long long a,long long b) -> long long { if(a == b) return a; long long s = a>b?a-b:b-a; int n = __builtin_ctzll(s); return f(f,s>>n,a>b?b:a); }; return f(f,a>>n,b>>m)<<(n>m?m:n); } template<typename T> vector<T> PollardsRho(T N,bool first = true){ if(N <= 1) return {}; vector<T> ret; while(N%2 == 0) N >>= 1,ret.push_back(2); if(N == 1) return ret; if(isprime(N)){ ret.push_back(N); return ret; } if(N <= 1024){ for(int i=3; i*i<=N; i++) while(N%i == 0) N /= i,ret.push_back(i); if(N != 1) ret.push_back(N); return ret; } mont::setmod(N); mont one = 1; for(int i=1; i<N; i++){ mont x1 = 0,y1 = 0,z1 = one; mont x2 = 0,y2 = 0,z2 = one; mont add1 = i*2-1,add2 = i*2; T g = 1; for(int r=512; ; r<<=1){ mont Y1 = y1,Y2 = y2; for(int t=0; t<r; t++){ y1 *= y1; y1 += add1; y2 *= y2; y2 += add2; z1 *= (x1-y1); z2 *= (x2-y2); } g = stein_gcd((z1*z2).val(),N); if(g == 1){x1 = y1; x2 = y2; continue;} if(g != N) break; T g1 = stein_gcd(z1.val(),N); if(g1 != 1 && g1 != N){g = g1; break;} T g2 = stein_gcd(z2.val(),N); if(g2 != 1 && g2 != N){g = g2; break;} g = 1; mont X = (g1==1)?x2:x1; mont Y = (g1==1)?y2:y1; mont a = (g1==1)?add2:add1; while(g == 1){ Y *= Y; Y += a; g = stein_gcd((X-Y).val(),N); } break; } if(g == N) continue; vector<T> P = PollardsRho(g,false),Q = PollardsRho(N/g,false); for(auto &p : P) ret.push_back(p); for(auto &q : Q) ret.push_back(q); break; } if(first) sort(ret.begin(),ret.end()); return ret; } int main(){ ios_base::sync_with_stdio(false); cin.tie(nullptr); long long N; cin >> N; auto P = PollardsRho(N); if(P.size() >= 3) cout << "YES\n"; else cout << "NO\n"; }