結果

問題 No.36 素数が嫌い！
コンテスト
ユーザー GOTKAKO
提出日時 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
権限があれば一括ダウンロードができます

ソースコード

diff # 

#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";
}
0