verify/LC-PrimalityTest.test.cpp

• View this file on GitHub
• Last update: 2026-06-01 13:41:50+09:00
• Problem: https://judge.yosupo.jp/problem/primality_test

Depends on

• Library/Common.hpp
• Library/Math/MillerRabin.hpp
• Library/Math/Montgomery.hpp
• Library/Template.hpp

Code

#define PROBLEM "https://judge.yosupo.jp/problem/primality_test"

#include "../Library/Template.hpp"
#include "../Library/Math/MillerRabin.hpp"

int main(){
    cin.tie(0)->sync_with_stdio(false);
    int Q; cin >> Q;
    while(Q--){
        long long N; cin >> N;
        YnPrint(MillerRabin(N));
    }
}

#line 1 "verify/LC-PrimalityTest.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/primality_test"

#line 2 "Library/Template.hpp"

#line 2 "Library/Common.hpp"

/**
 * @file Common.hpp
 */

#include <algorithm>
#include <array>
#include <bit>
#include <bitset>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <deque>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <stack>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
using namespace std;

using ll = int64_t;
using ull = uint64_t;

constexpr const ll INF = (1LL << 62) - (3LL << 30) - 1;
#line 4 "Library/Template.hpp"

inline bool YnPrint(bool flag){cout << (flag ? "Yes" : "No") << '\n'; return flag;}
inline bool YNPrint(bool flag){cout << (flag ? "YES" : "NO") << '\n'; return flag;}
template<typename Container>
inline void Sort(Container &container){sort(container.begin(), container.end());}
template<typename Container>
inline void ReverseSort(Container &container){sort(container.rbegin(), container.rend());}
template<typename Container>
inline void Reverse(Container &container){reverse(container.begin(), container.end());}
template<typename Value>
inline int PopCount(const Value &value){return __builtin_popcount(value);}
template<typename Value>
inline Value Floor(Value numerator, Value denominator){if(denominator < 0) numerator *= -1, denominator *= -1; return numerator < 0 ? (numerator + 1) / denominator - 1 : numerator / denominator;}
template<typename Value>
inline Value Ceil(Value numerator, Value denominator){if(denominator < 0) numerator *= -1, denominator *= -1; return numerator > 0 ? (numerator - 1) / denominator + 1 : numerator / denominator;}
template<typename Value>
inline int LowerBoundIndex(const vector<Value> &container, const Value &value){return distance(container.begin(), lower_bound(container.begin(), container.end(), value));}
template<typename Value>
inline int UpperBoundIndex(const vector<Value> &container, const Value &value){return distance(container.begin(), upper_bound(container.begin(), container.end(), value));}
template<typename Value>
inline bool Between(const Value &lower, const Value &x, const Value &higher){return lower <= x && x <= higher;}
template<typename Value>
inline bool InGrid(const Value &y, const Value &x, const Value &ymax, const Value &xmax){return Between(0, y, ymax - 1) && Between(0, x, xmax - 1);}
template<typename Value>
inline Value Median(const Value &a, const Value &b, const Value &c){return Between(b, a, c) || Between(c, a, b) ? a : (Between(a, b, c) || Between(c, b, a) ? b : c);}
template<typename Value>
inline Value Except(Value &src, Value &cond, Value &excp){return (src == cond ? excp : src);}

template<class Value>
bool chmin(Value &src, const Value &cmp){if(src > cmp){src = cmp; return true;} return false;}
template<class Value>
bool chmax(Value &src, const Value &cmp){if(src < cmp){src = cmp; return true;} return false;}
template<typename Value>
inline Value min(vector<Value> &v){return *min_element((v).begin(), (v).end());}
template<typename Value>
inline Value max(vector<Value> &v){return *max_element((v).begin(), (v).end());}

const int dx4[4] = {1, 0, -1, 0};
const int dy4[4] = {0, -1, 0, 1};
const int dx8[8] = {1, 1, 0, -1, -1, -1, 0, 1};
const int dy8[8] = {0, -1, -1, -1, 0, 1, 1, 1};

vector<pair<int, int>> adjacent(int current_y, int current_x, int max_y, int max_x, bool dir_8 = false){
    vector<pair<int, int>> ret;
    for(int d = 0; d < 4 * (1 + dir_8); ++d){
        int next_y = current_y + (dir_8 ? dy8[d] : dy4[d]);
        int next_x = current_x + (dir_8 ? dx8[d] : dx4[d]);
        if(InGrid(next_y, next_x, max_y, max_x)){
            ret.emplace_back(next_y, next_x);
        }
    }
    return ret;
}

template <typename T1, typename T2>
ostream &operator<<(ostream &os, const pair<T1, T2> &p){
    os << p.first << " " << p.second;
    return os;
}

template <typename T1, typename T2>
istream &operator>>(istream &is, pair<T1, T2> &p){
    is >> p.first >> p.second;
    return is;
}

template <typename T>
ostream &operator<<(ostream &os, vector<T> &v){
    for (int i = 0; i < v.size(); ++i){
        os << v[i] << (i + 1 != v.size() ? " " : "");
    }
    return os;
}

template <typename T>
ostream &operator<<(ostream &os, vector<vector<T>> &v){
    for (int i = 0; i < v.size(); ++i){
        os << v[i] << (i + 1 != v.size() ? "\n" : "");
    }
    return os;
}

template <typename T>
istream &operator>>(istream &is, vector<T> &v){
    for (int i = 0; i < v.size(); ++i) is >> v[i];
    return is;
}

template <typename T>
ostream &operator<<(ostream &os, set<T> &v){
    for (auto &u : v){
        os << u << " ";
    }
    return os;
}

template<typename T1, typename T2>
vector<pair<T1, T2>> AssembleVectorPair(vector<T1> &v1, vector<T2> &v2){
    assert(v1.size() == v2.size());
    vector<pair<T1, T2>> v;
    for(int i = 0; i < v1.size(); ++i) v.push_back({v1[i], v2[i]});
    return v;
}

template<typename T1, typename T2>
pair<vector<T1>, vector<T2>> DisassembleVectorPair(vector<pair<T1, T2>> &v){
    vector<T1> v1;
    vector<T2> v2;
    transform(v.begin(), v.end(), back_inserter(v1), [](auto p){return p.first;});
    transform(v.begin(), v.end(), back_inserter(v2), [](auto p){return p.second;});
    return {v1, v2};
}

template<typename T1, typename T2, typename T3>
tuple<vector<T1>, vector<T2>, vector<T3>> DisassembleVectorTuple(vector<tuple<T1, T2, T3>> &v){
    vector<T1> v1;
    vector<T2> v2;
    vector<T3> v3;
    transform(v.begin(), v.end(), back_inserter(v1), [](auto p){return get<0>(p);});
    transform(v.begin(), v.end(), back_inserter(v2), [](auto p){return get<1>(p);});
    transform(v.begin(), v.end(), back_inserter(v3), [](auto p){return get<2>(p);});
    return {v1, v2, v3};
}

template<typename T1 = int, typename T2 = T1>
pair<vector<T1>, vector<T2>> InputVectorPair(int size){
    vector<pair<T1, T2>> v(size);
    for(auto &[p, q] : v) cin >> p >> q;
    return DisassembleVectorPair(v);
}

template<typename T1 = int, typename T2 = T1, typename T3 = T1>
tuple<vector<T1>, vector<T2>, vector<T3>> InputVectorTuple(int size){
    vector<tuple<T1, T2, T3>> v(size);
    for(auto &[p, q, r] : v) cin >> p >> q >> r;
    return DisassembleVectorTuple(v);
}
#line 2 "Library/Math/Montgomery.hpp"

struct Montgomery{
    using mont = Montgomery;
    using u64 = uint64_t;
    using u128 = __uint128_t;

    static void SetMod(const u64 m){
        assert(m < (1LL << 62));
        assert(m & 1);
        mod = m;
        inv_mod = m;
        for(int i = 0; i < 5; ++i){
            inv_mod *= 2 - m * inv_mod;
        }
        r2 = -u128(m) % m;
    }

    u64 val;

    Montgomery() : val(0){}

    Montgomery(long long x) : val(Reduct((u128(x) + mod) * r2)){}

    inline u64 Get() const {
        u64 res = Reduct(val);
        return res >= mod ? res - mod : res;
    }

    mont &operator+=(const mont &r){
        if((val += r.val) >= 2 * mod) val -= 2 * mod;
        return (*this);
    }

    mont &operator-=(const mont &r){
        if((val += 2 * mod - r.val) >= 2 * mod) val -= 2 * mod;
        return (*this);
    }

    mont &operator*=(const mont &r){
        val = Reduct(u128(val) * r.val);
        return (*this);
    }

    mont &operator/=(const mont &r){
        (*this) *= r.Inverse();
        return (*this);
    }

    mont Inverse() const {
        return Power(mod - 2);
    }

    mont Power(long long k) const {
        mont ret(1), mul(*this);
        while(k){
            if(k & 1) ret *= mul;
            mul *= mul;
            k >>= 1;
        }
        return ret;
    }

    mont operator-() const {
        return mont() - mont(*this);
    }

    mont operator+(const mont &r) const {
        return mont(*this) += r;
    }

    mont operator-(const mont &r) const {
        return mont(*this) -= r;
    }

    mont operator*(const mont &r) const {
        return mont(*this) *= r;
    }

    bool operator==(const mont &r) const {
        return (val >= mod ? val - mod : val) == (r.val >= mod ? r.val - mod : r.val);
    }

    bool operator!=(const mont &r) const {
        return (val >= mod ? val - mod : val) != (r.val >= mod ? r.val - mod : r.val);
    }

    private:
    static u64 mod;
    static u64 inv_mod;
    static u64 r2;

    static u64 Reduct(const u128 x){
        return (x + u128(u64(x) * u64(-inv_mod)) * mod) >> 64;
    }
};

typename Montgomery::u64 Montgomery::mod, Montgomery::inv_mod, Montgomery::r2;
#line 2 "Library/Math/MillerRabin.hpp"

bool MillerRabin(const long long N){
    if(N <= 1) return false;
    if(N == 2) return true;
    if(!(N & 1)) return false;
    vector<long long> as;
    if(N < 4759123141LL) as = {2, 7, 61};
    else as = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    Montgomery::SetMod(N);
    long long s = 0, d = N - 1;
    while(!(d & 1)) ++s, d >>= 1;
    for(const auto &a : as){
        Montgomery x = Montgomery(a).Power(d);
        if(x != 1 && x != 0){
            long long t;
            for(t = 0; t < s; ++t){
                if(x == N - 1) break;
                x *= x;
            }
            if(t == s) return false;
        }
    }
    return true;
}
#line 5 "verify/LC-PrimalityTest.test.cpp"

int main(){
    cin.tie(0)->sync_with_stdio(false);
    int Q; cin >> Q;
    while(Q--){
        long long N; cin >> N;
        YnPrint(MillerRabin(N));
    }
}

Back to top page