fork(1) download
copy 
  1. #pragma GCC optimize("O3,unroll-loops")
  2. #include <iostream>
  3. #include <vector>
  4. #include <string>
  5. #include <algorithm>
  6.  
  7. using namespace std;
  8.  
  9. // ================= 大數類別 (Base 2^32) =================
  10. struct BigInt {
  11.     vector<uint32_t> d;
  12.  
  13.     BigInt() {}
  14.     BigInt(unsigned long long v) {
  15.         while (v > 0) {
  16.             d.push_back(v & 0xFFFFFFFF);
  17.             v >>= 32;
  18.         }
  19.     }
  20.  
  21.     bool is_zero() const { return d.empty(); }
  22.  
  23.     void trim() {
  24.         while (!d.empty() && d.back() == 0) d.pop_back();
  25.     }
  26.  
  27.     bool is_even() const { return d.empty() || (d[0] & 1) == 0; }
  28.     bool is_odd() const { return !d.empty() && (d[0] & 1) != 0; }
  29.  
  30.     int bits() const {
  31.         if (d.empty()) return 0;
  32.         int res = (d.size() - 1) * 32;
  33.         uint32_t top = d.back();
  34.         while (top) { res++; top >>= 1; }
  35.         return res;
  36.     }
  37.  
  38.     bool get_bit(int i) const {
  39.         int word = i / 32;
  40.         int bit = i % 32;
  41.         if (word >= d.size()) return false;
  42.         return (d[word] >> bit) & 1;
  43.     }
  44.  
  45.     void set_bit(int i, bool val) {
  46.         int word = i / 32;
  47.         int bit = i % 32;
  48.         while (d.size() <= word) d.push_back(0);
  49.         if (val) d[word] |= (1u << bit);
  50.         else d[word] &= ~(1u << bit);
  51.         trim();
  52.     }
  53.  
  54.     bool operator==(const BigInt& o) const { return d == o.d; }
  55.     bool operator!=(const BigInt& o) const { return d != o.d; }
  56.     bool operator<(const BigInt& o) const {
  57.         if (d.size() != o.d.size()) return d.size() < o.d.size();
  58.         for (int i = (int)d.size() - 1; i >= 0; --i) {
  59.             if (d[i] != o.d[i]) return d[i] < o.d[i];
  60.         }
  61.         return false;
  62.     }
  63.     bool operator<=(const BigInt& o) const { return !(o < *this); }
  64.     bool operator>(const BigInt& o) const { return o < *this; }
  65.     bool operator>=(const BigInt& o) const { return !(*this < o); }
  66.  
  67.     BigInt operator+(const BigInt& o) const {
  68.         BigInt res;
  69.         uint64_t carry = 0;
  70.         int n = max(d.size(), o.d.size());
  71.         for (int i = 0; i < n || carry; ++i) {
  72.             uint64_t sum = carry;
  73.             if (i < d.size()) sum += d[i];
  74.             if (i < o.d.size()) sum += o.d[i];
  75.             res.d.push_back(sum & 0xFFFFFFFF);
  76.             carry = sum >> 32;
  77.         }
  78.         res.trim();
  79.         return res;
  80.     }
  81.  
  82.     BigInt operator-(const BigInt& o) const {
  83.         BigInt res;
  84.         int64_t borrow = 0;
  85.         for (int i = 0; i < d.size(); ++i) {
  86.             int64_t diff = d[i] - borrow;
  87.             if (i < o.d.size()) diff -= o.d[i];
  88.             if (diff < 0) {
  89.                 diff += 0x100000000LL;
  90.                 borrow = 1;
  91.             } else {
  92.                 borrow = 0;
  93.             }
  94.             res.d.push_back(diff);
  95.         }
  96.         res.trim();
  97.         return res;
  98.     }
  99.  
  100.     void shift_left_1() {
  101.         if (d.empty()) return;
  102.         uint32_t carry = 0;
  103.         for (int i = 0; i < d.size(); ++i) {
  104.             uint32_t next_carry = d[i] >> 31;
  105.             d[i] = (d[i] << 1) | carry;
  106.             carry = next_carry;
  107.         }
  108.         if (carry) d.push_back(carry);
  109.     }
  110.  
  111.     void shift_right_1() {
  112.         if (d.empty()) return;
  113.         uint32_t carry = 0;
  114.         for (int i = (int)d.size() - 1; i >= 0; --i) {
  115.             uint32_t next_carry = (d[i] & 1) << 31;
  116.             d[i] = (d[i] >> 1) | carry;
  117.             carry = next_carry;
  118.         }
  119.         trim();
  120.     }
  121.  
  122.     BigInt operator*(const BigInt& o) const {
  123.         if (is_zero() || o.is_zero()) return BigInt();
  124.         BigInt res;
  125.         res.d.assign(d.size() + o.d.size(), 0);
  126.         for (size_t i = 0; i < d.size(); ++i) {
  127.             uint64_t carry = 0;
  128.             for (size_t j = 0; j < o.d.size() || carry; ++j) {
  129.                 uint64_t cur = res.d[i+j] + d[i] * 1ULL * (j < o.d.size() ? o.d[j] : 0) + carry;
  130.                 res.d[i+j] = cur & 0xFFFFFFFF;
  131.                 carry = cur >> 32;
  132.             }
  133.         }
  134.         res.trim();
  135.         return res;
  136.     }
  137.  
  138.     // 二進位長除法 (僅用於 Barrett 預處理)
  139.     friend void divmod(const BigInt& a, const BigInt& b, BigInt& q, BigInt& r) {
  140.         q.d.clear(); r.d.clear();
  141.         if (a < b) { r = a; return; }
  142.         for (int i = a.bits() - 1; i >= 0; --i) {
  143.             r.shift_left_1();
  144.             if (a.get_bit(i)) {
  145.                 if (r.d.empty()) r.d.push_back(1);
  146.                 else r.d[0] |= 1;
  147.             }
  148.             if (r >= b) {
  149.                 r = r - b;
  150.                 q.set_bit(i, true);
  151.             }
  152.         }
  153.     }
  154.  
  155.     long long mod_small(long long b) const {
  156.         long long rem = 0;
  157.         long long base = (1ULL << 32) % b;
  158.         for (int i = (int)d.size() - 1; i >= 0; --i) {
  159.             rem = (rem * base + d[i]) % b;
  160.         }
  161.         return rem;
  162.     }
  163. };
  164.  
  165. // ================= Barrett Reduction =================
  166. struct Barrett {
  167.     BigInt N, mu;
  168.     int k;
  169.  
  170.     Barrett(BigInt n) : N(n) {
  171.         k = N.d.size();
  172.         BigInt b2k;
  173.         b2k.set_bit(2 * k * 32, true);
  174.         BigInt rem;
  175.         divmod(b2k, N, mu, rem);
  176.     }
  177.  
  178.     BigInt reduce(BigInt A) const {
  179.         BigInt q1;
  180.         if (A.d.size() > k - 1) q1.d.assign(A.d.begin() + k - 1, A.d.end());
  181.         BigInt q2 = q1 * mu;
  182.         BigInt q3;
  183.         if (q2.d.size() > k + 1) q3.d.assign(q2.d.begin() + k + 1, q2.d.end());
  184.  
  185.         BigInt r1;
  186.         int r1_size = min((int)A.d.size(), k + 1);
  187.         r1.d.assign(A.d.begin(), A.d.begin() + r1_size);
  188.         r1.trim();
  189.  
  190.         BigInt r2 = q3 * N;
  191.         int r2_size = min((int)r2.d.size(), k + 1);
  192.         r2.d.resize(r2_size);
  193.         r2.trim();
  194.  
  195.         BigInt r;
  196.         if (r1 >= r2) r = r1 - r2;
  197.         else {
  198.             BigInt bk1;
  199.             bk1.set_bit((k + 1) * 32, true);
  200.             r = (r1 + bk1) - r2;
  201.         }
  202.         while (r >= N) r = r - N;
  203.         return r;
  204.     }
  205.  
  206.     BigInt mul(const BigInt& a, const BigInt& b) const {
  207.         return reduce(a * b);
  208.     }
  209.  
  210.     BigInt power(BigInt base, BigInt exp) const {
  211.         BigInt res = 1;
  212.         base = reduce(base);
  213.         for (int i = 0; i < exp.bits(); ++i) {
  214.             if (exp.get_bit(i)) res = mul(res, base);
  215.             base = mul(base, base);
  216.         }
  217.         return res;
  218.     }
  219. };
  220.  
  221. // ================= BPSW 演算法 =================
  222. BigInt parse(const string& s) {
  223.     BigInt res;
  224.     int i = 0;
  225.     while (i < s.length()) {
  226.         int len = min((int)s.length() - i, 9);
  227.         long long val = stoll(s.substr(i, len));
  228.         long long multiplier = 1;
  229.         for(int j = 0; j < len; ++j) multiplier *= 10;
  230.         res = res * BigInt(multiplier) + BigInt(val);
  231.         i += len;
  232.     }
  233.     return res;
  234. }
  235.  
  236. int jacobi(long long a, const BigInt& n) {
  237.     int t = 1;
  238.     long long n8 = n.mod_small(8);
  239.     if (a < 0) {
  240.         a = -a;
  241.         if (n8 % 4 == 3) t = -t;
  242.     }
  243.     long long n_val = n.mod_small(a);
  244.     while (a != 0) {
  245.         while (a % 2 == 0) {
  246.             a /= 2;
  247.             if (n8 == 3 || n8 == 5) t = -t;
  248.         }
  249.         swap(a, n_val);
  250.         if (a % 4 == 3 && n_val % 4 == 3) t = -t;
  251.         a %= n_val;
  252.     }
  253.     if (n_val == 1) return t;
  254.     return 0;
  255. }
  256.  
  257. bool is_prime(BigInt n) {
  258.     if (n == BigInt(2)) return true;
  259.     if (n < BigInt(2) || n.is_even()) return false;
  260.  
  261.     // 1. 小質數試除 (加速過濾)
  262.     static const int primes[] = {3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97};
  263.     for (int p : primes) {
  264.         if (n == BigInt(p)) return true;
  265.         if (n.mod_small(p) == 0) return false;
  266.     }
  267.  
  268.     Barrett barrett(n);
  269.  
  270.     // 2. Miller-Rabin (Base 2)
  271.     BigInt d = n - BigInt(1);
  272.     int s = 0;
  273.     while (d.is_even()) {
  274.         d.shift_right_1();
  275.         s++;
  276.     }
  277.  
  278.     BigInt x = barrett.power(BigInt(2), d);
  279.     bool mr_pass = false;
  280.     if (x == BigInt(1) || x == n - BigInt(1)) {
  281.         mr_pass = true;
  282.     } else {
  283.         for (int i = 1; i < s; i++) {
  284.             x = barrett.mul(x, x);
  285.             if (x == n - BigInt(1)) {
  286.                 mr_pass = true;
  287.                 break;
  288.             }
  289.         }
  290.     }
  291.     if (!mr_pass) return false;
  292.  
  293.     // 3. 尋找 Lucas 參數 D (Selfridge 方法)
  294.     long long D = 5;
  295.     int sign = 1;
  296.     while (true) {
  297.         long long curD = D * sign;
  298.         int j = jacobi(curD, n);
  299.         if (j == -1) {
  300.             D = curD;
  301.             break;
  302.         }
  303.         if (j == 0) return false;
  304.         D += 2;
  305.         sign = -sign;
  306.     }
  307.  
  308.     // 4. Strong Lucas Test
  309.     long long Q = (1 - D) / 4;
  310.     BigInt Q_mod = (Q < 0) ? n - BigInt(-Q) : BigInt(Q);
  311.     BigInt D_mod = (D < 0) ? n - BigInt(-D) : BigInt(D);
  312.  
  313.     BigInt U = BigInt(1), V = BigInt(1);
  314.     BigInt Qk = Q_mod;
  315.  
  316.     auto div2_mod = [&](BigInt A) {
  317.         if (A.is_odd()) A = A + n;
  318.         A.shift_right_1();
  319.         if (A >= n) A = A - n;
  320.         return A;
  321.     };
  322.  
  323.     BigInt K = n + BigInt(1);
  324.     BigInt d_lucas = K;
  325.     int s_lucas = 0;
  326.     while (d_lucas.is_even()) {
  327.         d_lucas.shift_right_1();
  328.         s_lucas++;
  329.     }
  330.  
  331.     for (int i = d_lucas.bits() - 2; i >= 0; --i) {
  332.         BigInt U_new = barrett.mul(U, V);
  333.         BigInt V_new = barrett.mul(V, V) - barrett.mul(BigInt(2), Qk);
  334.         if (V_new >= n) V_new = V_new + n; // 處理負數
  335.         Qk = barrett.mul(Qk, Qk);
  336.  
  337.         U = U_new;
  338.         V = V_new;
  339.  
  340.         if (d_lucas.get_bit(i)) {
  341.             BigInt DU = barrett.mul(D_mod, U);
  342.             U = div2_mod(U + V);
  343.             V = div2_mod(DU + V);
  344.             Qk = barrett.mul(Qk, Q_mod);
  345.         }
  346.     }
  347.  
  348.     if (U.is_zero() || V.is_zero()) return true;
  349.  
  350.     for (int r = 1; r < s_lucas; ++r) {
  351.         V = barrett.mul(V, V) - barrett.mul(BigInt(2), Qk);
  352.         if (V >= n) V = V + n;
  353.         Qk = barrett.mul(Qk, Qk);
  354.         if (V.is_zero()) return true;
  355.     }
  356.  
  357.     return false;
  358. }
  359.  
  360. int main() {
  361.     ios_base::sync_with_stdio(false);
  362.     cin.tie(NULL);
  363.     string s;
  364.     while (cin >> s) {
  365.         if (is_prime(parse(s))) {
  366.             cout << "YES\n";
  367.         } else {
  368.             cout << "NO\n";
  369.         }
  370.     }
  371.     return 0;
  372. }
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555151555
stdout
copy 
NO
https://ideone.com/mDN2QU
language:
C++14 (gcc 8.3)
created:
10 hours ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!