結果
問題 | No.577 Prime Powerful Numbers |
ユーザー |
![]() |
提出日時 | 2017-10-13 22:25:59 |
言語 | C++14 (gcc 13.3.0 + boost 1.87.0) |
結果 |
AC
|
実行時間 | 187 ms / 2,000 ms |
コード長 | 7,800 bytes |
コンパイル時間 | 1,277 ms |
コンパイル使用メモリ | 109,960 KB |
実行使用メモリ | 6,820 KB |
最終ジャッジ日時 | 2024-11-17 17:55:37 |
合計ジャッジ時間 | 2,035 ms |
ジャッジサーバーID (参考情報) |
judge1 / judge5 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
sample | AC * 1 |
other | AC * 10 |
ソースコード
#include <cstdio>#include <cassert>#include <cmath>#include <cstring>#include <iostream>#include <algorithm>#include <vector>#include <map>#include <set>#include <functional>#include <stack>#include <queue>#include <tuple>#define getchar getchar_unlocked#define putchar putchar_unlocked#define _rep(_1, _2, _3, _4, name, ...) name#define rep2(i, n) rep3(i, 0, n)#define rep3(i, a, b) rep4(i, a, b, 1)#define rep4(i, a, b, c) for (int i = int(a); i < int(b); i += int(c))#define rep(...) _rep(__VA_ARGS__, rep4, rep3, rep2, _)(__VA_ARGS__)using namespace std;using i64 = long long;using u8 = unsigned char;using u32 = unsigned;using u64 = unsigned long long;using f80 = long double;int get_int() {int c, n;while ((c = getchar()) < '0');n = c - '0';while ((c = getchar()) >= '0') n = n * 10 + (c - '0');return n;}namespace factor {using u128 = __uint128_t;template <typename word, typename dword, typename sword>struct Mod {Mod() : n_(0) {}Mod(word n) : n_(init(n)) {}static constexpr int word_bits = sizeof(word) * 8;static void set_mod(word m) { mod = m, inv = mul_inv(m), r2 = -dword(m) % m; }static word mul_inv(word n) {word x = n;rep(_, 5) x *= 2 - n * x;return x;}static word reduce(dword w) {word y = (w >> word_bits) - ((dword(word(w) * inv) * mod) >> word_bits);return sword(y) < 0 ? y + mod : y;}static word init(word n) { return reduce(u128(n) * r2); }static word ilog2(word n) { return (n == 0) ? 0 : __lg(n); }Mod& operator += (Mod rhs) { if ((n_ += rhs.n_) >= mod) n_ -= mod; return *this; }Mod& operator -= (Mod rhs) { if (sword(n_ -= rhs.n_) < 0) n_ += mod; return *this; }Mod& operator *= (Mod rhs) { n_ = reduce(u128(n_) * rhs.n_); return *this; }bool operator == (Mod rhs) { return n_ == rhs.n_; }bool operator != (Mod rhs) { return !(*this == rhs); }Mod operator + (Mod rhs) { return Mod(*this) += rhs; }Mod operator - (Mod rhs) { return Mod(*this) -= rhs; }Mod operator * (Mod rhs) { return Mod(*this) *= rhs; }Mod operator - () { return Mod() - *this; };Mod pow(word e) {Mod ret = Mod(1), base = *this;for (; e; e >>= 1, base *= base) {if (e & 1) ret *= base;}return ret;}word val() { return reduce(n_); }friend ostream& operator << (ostream& os, Mod& m) { return os << m.val(); }word n_;static word mod, inv, r2;};using m64 = Mod<u64, u128, i64>;using m32 = Mod<u32, u64, int>;template <> u32 m32::mod = 0;template <> u32 m32::inv = 0;template <> u32 m32::r2 = 0;template <> u64 m64::mod = 0;template <> u64 m64::inv = 0;template <> u64 m64::r2 = 0;template <typename word>word gcd(word a, word b) {while (b) { word t = a % b; a = b; b = t; }return a;}template <typename word, typename mod>word brent(word n, word c) {// n must be composite and odd.const u64 s = 256;mod::set_mod(n);const mod one = mod(1), mc = mod(c);auto f = [&] (mod x) { return x * x + mc; };mod y = one;for (u64 l = 1; ; l <<= 1) {auto x = y;rep(_, l) y = f(y);mod p = one;rep(k, 0, l, s) {auto sy = y;rep(_, min(s, l - k)) y = f(y), p *= y - x;word g = gcd(n, p.n_);if (g == 1) continue;if (g == n) for (g = 1, y = sy; g == 1; ) y = f(y), g = gcd(n, (y - x).n_);return g;}}}u64 brent(u64 n, u64 c) {if (n < (u32(1) << 31)) {return brent<u32, m32>(n, c);} else if (n < (u64(1) << 63)) {return brent<u64, m64>(n, c);} else {assert(0);}}template <typename word, typename mod>bool composite(word n, word d, int s, const u32* bases, int base_size) {mod::set_mod(n);mod one = mod(1), minus_one = -one;rep(i, base_size) {mod b = mod(bases[i]).pow(d);if (b == one || b == minus_one) continue;int t = s - 1;for (; t > 0; --t) if ((b *= b) == minus_one) break;if (t == 0) return true;}return false;}bool miller_rabin(u64 n) {static const u32 bases[][7] {{2, 3},{2, 299417},{2, 7, 61},{15, 176006322, u32(4221622697)},{2, 2570940, 211991001, u32(3749873356)},{2, 2570940, 880937, 610386380, u32(4130785767)},{2, 325, 9375, 28178, 450775, 9780504, 1795265022},};if (n <= 1) return false;if (!(n & 1)) return n == 2;if (n <= 8) return true;u64 d = n - 1;u64 s = __builtin_ctzll(d);d >>= s;u32 base_id = 6, base_size = 7;if (n < 1373653) {base_id = 0, base_size = 2;} else if (n < 19471033) {base_id = 1, base_size = 2;} else if (n < u64(4759123141)) {base_id = 2, base_size = 3;} else if (n < u64(154639673381)) {base_id = 3, base_size = 3;} else if (n < u64(47636622961201)) {base_id = 4, base_size = 4;} else if (n < u64(3770579582154547)) {base_id = 5, base_size = 5;}if (n < (u32(1) << 31)) {return !composite<u32, m32>(n, d, s, bases[base_id], base_size);} else if (n < (u64(1) << 63)) {return !composite<u64, m64>(n, d, s, bases[base_id], base_size);} else assert(0);return true;}struct ExactDiv {ExactDiv() {}ExactDiv(u64 n) : n(n), i(m64::mul_inv(n)), t(u64(-1) / n) {}friend u64 operator / (u64 n, ExactDiv d) { return n * d.i; };bool divide(u64 n) { return n / *this <= t; }u64 n, i, t;};vector<ExactDiv> primes;void init(u32 n) {u32 sqrt_n = sqrt(n);vector<u8> isprime(n + 1, 1);rep(i, 2, sqrt_n + 1) if (isprime[i]) rep(j, i * i, n + 1, i) isprime[j] = 0;primes.clear();rep(i, 2, n + 1) if (isprime[i]) primes.push_back(ExactDiv(i));}u32 isqrt(u64 n) {return sqrtl(n);}u64 square(u64 n) {return n * n;}u64 ctz(u64 n) {return __builtin_ctzll(n);}bool is_prime(u64 n) {if (n <= 1) return false;if (!(n & 1)) return n == 2;rep(i, 1, primes.size()) {auto p = primes[i];if (square(p.n) > n) return true;if (p.divide(n)) return n == p.n;}return miller_rabin(n);}u64 next_prime(u64 n) {if (n <= 1) return 2;if (n <= 2) return 3;for (n = (n + 1) / 2 * 2 + 1; !is_prime(n); n += 2);return n;}using P = pair<u64, u32>;vector<P> factors(u64 n) {assert(n < (u64(1) << 63));auto ret = vector<P>();if (n <= 1) {return ret;}u32 v = isqrt(n);if (u64(v) * v == n) {auto res = factors(v);for (auto& pe : res) pe.second <<= 1;return res;}if (!(n & 1)) {u32 e = ctz(n);ret.emplace_back(2, e);n >>= e;}u64 lim = square(primes[primes.size() - 1].n);rep(pi, 1, primes.size()) {auto p = primes[pi];if (square(p.n) > n) break;if (p.divide(n)) {u32 e = 1; n = n / p;while (p.divide(n)) n = n / p, e++;ret.emplace_back(p.n, e);}}u32 s = ret.size();while (n > lim && !miller_rabin(n)) {for (u64 c = 1; ; ++c) {u64 p = brent(n, c);if (!miller_rabin(p)) continue;u32 e = 1; n /= p;while (n % p == 0) n /= p, e += 1;ret.emplace_back(p, e);break;}}if (n > 1) ret.emplace_back(n, 1);if (ret.size() - s >= 2) sort(ret.begin() + s, ret.end());return ret;}} // namespace factorvoid solve() {factor::init(512);int T = get_int();rep(_, T) {i64 n; scanf("%lld", &n);if (n == 2) {puts("No");} else if (n % 2 == 0) {puts("Yes");} else {i64 two = 2;bool found = false;while (two <= n) {i64 m = n - two;auto f = factor::factors(m);if (f.size() == 1) {found = true;break;}two <<= 1;}puts(found ? "Yes" : "No");}}}int main() {clock_t beg = clock();solve();clock_t end = clock();fprintf(stderr, "%.3f sec\n", double(end - beg) / CLOCKS_PER_SEC);return 0;}