結果
| 問題 |
No.8030 ミラー・ラビン素数判定法のテスト
|
| ユーザー |
 Rac
|
| 提出日時 | 2024-12-30 21:08:46 |
| 言語 | C++23 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
RE
|
| 実行時間 | - |
| コード長 | 34,507 bytes |
| コンパイル時間 | 1,477 ms |
| コンパイル使用メモリ | 70,660 KB |
| 実行使用メモリ | 5,248 KB |
| 最終ジャッジ日時 | 2024-12-30 21:08:49 |
| 合計ジャッジ時間 | 2,764 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge5 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 7 WA * 1 RE * 2 |
ソースコード
/**
* libcpprime https://github.com/sortA0329/libcpprime
*
* MIT License
*
* Copyright (c) 2024 Rac
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
**/
#ifndef LIBCPPRIME_INCLUDED_IS_PRIME_NO_TABLE
#define LIBCPPRIME_INCLUDED_IS_PRIME_NO_TABLE
/**
* libcpprime https://github.com/sortA0329/libcpprime
*
* MIT License
*
* Copyright (c) 2024 Rac
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
**/
/**
*
* The algorithm in this library is based on Bradley Berg's method.
* See this page for more information: https://www.techneon.com/download/is.prime.32.base.data
*
* Copyright 2018 Bradley Berg < (My last name) @ t e c h n e o n . c o m >
*
* Permission to use, copy, modify, and distribute this software for any
* purpose with or without fee is hereby granted, provided that the above
* copyright notice and this permission notice appear in all copies.
*
* This algorithm is deliberately unpatented. The license above also
* lets you even freely use it in commercial code.
*
* Primality testing using a hash table of bases originated with Steven Worley.
*
**/
#ifndef LIBCPPRIME_INCLUDED_IS_PRIME_COMMON
#define LIBCPPRIME_INCLUDED_IS_PRIME_COMMON
#include <cstdint>
#ifdef __has_include
#if __has_include(<type_traits>)
#include <type_traits>
#endif
#if __has_include(<bit>)
#include <bit>
#endif
#endif
#ifdef _MSC_VER
#include <intrin.h>
#include <immintrin.h>
#pragma intrinsic(_tzcnt_u64, _lzcnt_u64, _umul128, __umulh)
#if _MSC_VER >= 1900
#pragma intrinsic(_udiv128)
#endif
#endif
#if defined(__cpp_lib_is_constant_evaluated) && defined(__cpp_lib_bitops)
#define LIBCPPRIME_CONSTEXPR constexpr
#else
#define LIBCPPRIME_CONSTEXPR
#endif
#ifdef __cpp_if_constexpr
#define LIBCPPRIME_IF_CONSTEXPR constexpr
#else
#define LIBCPPRIME_IF_CONSTEXPR
#endif
namespace cppr {
namespace internal {
struct Int64Pair {
std::uint64_t high, low;
};
LIBCPPRIME_CONSTEXPR void Assume(const bool f) noexcept {
#ifdef __cpp_lib_is_constant_evaluated
if (std::is_constant_evaluated()) return;
#endif
#if defined(__clang__)
#if __has_builtin(__builtin_assume)
__builtin_assume(f);
#endif
#elif defined(__GNUC__)
if (!f) __builtin_unreachable();
#elif defined(_MSC_VER)
__assume(f);
#endif
}
LIBCPPRIME_CONSTEXPR std::int32_t CountrZero(std::uint64_t n) noexcept {
#ifdef __cpp_lib_bitops
Assume(n != 0);
return std::countr_zero(n);
#elif defined(__GNUC__)
return __builtin_ctzll(n);
#elif defined(_MSC_VER)
return _tzcnt_u64(n);
#else
std::uint64_t x = (n ^ (n - 1)) >> 1;
x = x - ((x >> 1) & 0x5555555555555555);
x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333);
x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f;
x = x + (x >> 8);
x = x + (x >> 16);
x = x + (x >> 32);
return x & 0x0000007f;
#endif
}
#if !defined(__cpp_lib_bitops) && !defined(__GNUC__) && !defined(_MSC_VER)
LIBCPPRIME_CONSTEXPR const char CountlZeroTable[32] = { 0, 31, 9, 30, 3, 8, 13, 29, 2, 5, 7, 21, 12, 24, 28, 19, 1, 10, 4, 14, 6, 22, 25, 20, 11, 15, 23, 26, 16, 27, 17, 18 };
#endif
LIBCPPRIME_CONSTEXPR std::int32_t CountlZero(std::uint64_t n) noexcept {
#ifdef __cpp_lib_bitops
Assume(n != 0);
return std::countl_zero(n);
#elif defined(__GNUC__)
return __builtin_clzll(n);
#elif defined(_MSC_VER)
return _lzcnt_u64(n);
#else
std::int32_t offset = 0;
std::uint32_t x = n >> 32;
if (n <= 0xffffffff) {
offset = 32;
x = n;
}
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
x++;
return offset + CountlZeroTable[(x * 0x076be629u) >> 27];
#endif
}
LIBCPPRIME_CONSTEXPR Int64Pair Mulu128(std::uint64_t muler, std::uint64_t mulnd) noexcept {
#if defined(__SIZEOF_INT128__)
__uint128_t tmp = static_cast<__uint128_t>(muler) * mulnd;
return { static_cast<std::uint64_t>(tmp >> 64), static_cast<std::uint64_t>(tmp) };
#else
#if defined(_MSC_VER)
#ifdef __cpp_lib_is_constant_evaluated
if (!std::is_constant_evaluated())
#endif
{
std::uint64_t high;
std::uint64_t low = _umul128(muler, mulnd, &high);
return { high, low };
}
#endif
std::uint64_t u1 = (muler & 0xffffffff);
std::uint64_t v1 = (mulnd & 0xffffffff);
std::uint64_t t = (u1 * v1);
std::uint64_t w3 = (t & 0xffffffff);
std::uint64_t k = (t >> 32);
muler >>= 32;
t = (muler * v1) + k;
k = (t & 0xffffffff);
std::uint64_t w1 = (t >> 32);
mulnd >>= 32;
t = (u1 * mulnd) + k;
k = (t >> 32);
return { (muler * mulnd) + w1 + k, (t << 32) + w3 };
#endif
}
LIBCPPRIME_CONSTEXPR std::uint64_t Mulu128High(std::uint64_t muler, std::uint64_t mulnd) noexcept {
#if defined(__SIZEOF_INT128__)
return static_cast<std::uint64_t>((static_cast<__uint128_t>(muler) * mulnd) >> 64);
#else
#if defined(_MSC_VER)
#ifdef __cpp_lib_is_constant_evaluated
if (!std::is_constant_evaluated())
#endif
return __umulh(muler, mulnd);
#endif
return Mulu128(muler, mulnd).high;
#endif
}
LIBCPPRIME_CONSTEXPR Int64Pair Divu128(std::uint64_t high, std::uint64_t low, std::uint64_t div) noexcept {
#if (defined(__GNUC__) || defined(__ICC)) && defined(__x86_64__)
if LIBCPPRIME_IF_CONSTEXPR (sizeof(void*) == 8) {
#ifdef __cpp_lib_is_constant_evaluated
if (!std::is_constant_evaluated())
#endif
{
std::uint64_t res = 0, rem = 0;
__asm__("divq %[v]" : "=a"(res), "=d"(rem) : [v] "r"(div), "a"(low), "d"(high));
return { res, rem };
}
}
#elif defined(_MSC_VER) && _MSC_VER >= 1900
#ifdef __cpp_lib_is_constant_evaluated
if (!std::is_constant_evaluated())
#endif
{
std::uint64_t rem = 0;
std::uint64_t res = _udiv128(high, low, div, &rem);
return { res, rem };
}
#endif
#if defined(__SIZEOF_INT128__)
__uint128_t n = (static_cast<__uint128_t>(high) << 64 | low);
__uint128_t res = n / div;
return { static_cast<std::uint64_t>(res), static_cast<std::uint64_t>(n - res * div) };
#else
std::uint64_t res = 0;
std::uint64_t cur = high;
for (std::uint64_t i = 0; i != 64; ++i) {
std::uint64_t large = cur >> 63;
cur = cur << 1 | (low >> 63);
low <<= 1;
large |= (cur >= div);
res = res << 1 | large;
cur -= div & (0 - large);
}
return { res, cur };
#endif
}
LIBCPPRIME_CONSTEXPR std::uint32_t GCD32(std::uint32_t x, std::uint32_t y) noexcept {
if (x == 0) return 0;
const std::int32_t n = CountrZero(x);
const std::int32_t m = CountrZero(y);
const std::int32_t l = n < m ? n : m;
x >>= n;
y >>= m;
while (x != y) {
const std::uint32_t a = y - x, b = x - y;
const std::int32_t m = CountrZero(a), n = CountrZero(b);
Assume(m == n);
const std::uint32_t s = y < x ? b : a;
const std::uint32_t t = x < y ? x : y;
x = s >> m;
y = t;
}
return x << l;
}
template<std::int32_t Strict = 0> class MontgomeryModint64Impl {
std::uint64_t mod_ = 0, rs = 0, nr = 0, np = 0;
LIBCPPRIME_CONSTEXPR std::uint64_t reduce(const std::uint64_t n) const noexcept {
std::uint64_t q = n * nr;
if LIBCPPRIME_IF_CONSTEXPR (Strict == 1) {
std::uint64_t m = Mulu128High(q, mod_);
return m <= mod_ ? mod_ - m : 0 - m;
} else if LIBCPPRIME_IF_CONSTEXPR (Strict == 2) {
auto m = Mulu128(q, mod_);
std::uint64_t t = m.high + std::uint64_t(m.low + n < m.low);
return t - mod_ * (t >= mod_);
} else {
std::uint64_t m = Mulu128High(q, mod_);
return mod_ - m;
}
}
LIBCPPRIME_CONSTEXPR std::uint64_t reduce(const std::uint64_t a, const std::uint64_t b) const noexcept {
auto tmp = Mulu128(a, b);
std::uint64_t d = tmp.high, c = tmp.low;
std::uint64_t q = c * nr;
if LIBCPPRIME_IF_CONSTEXPR (Strict == 1) {
std::uint64_t m = Mulu128High(q, mod_);
std::uint64_t t = d - m;
return t >= mod_ ? t + mod_ : t;
} else if LIBCPPRIME_IF_CONSTEXPR (Strict == 2) {
auto m = Mulu128(q, mod_);
std::uint64_t t = m.high + std::uint64_t(m.low + c < m.low);
if (t >= mod_ - d) return t + d - mod_;
else return t + d;
} else {
std::uint64_t m = Mulu128High(q, mod_);
return d + mod_ - m;
}
}
public:
LIBCPPRIME_CONSTEXPR MontgomeryModint64Impl() noexcept {}
LIBCPPRIME_CONSTEXPR void set(std::uint64_t n) noexcept {
Assume(n > 2 && n % 2 != 0);
mod_ = n;
rs = Divu128(0xffffffffffffffff % n, 0 - n, n).low;
nr = n;
for (std::uint32_t i = 0; i != 5; ++i) nr *= 2 - n * nr;
if LIBCPPRIME_IF_CONSTEXPR (Strict == 2) nr = 0 - nr;
np = reduce(rs);
}
LIBCPPRIME_CONSTEXPR std::uint64_t build(std::uint32_t x) const noexcept { return reduce(x % mod_, rs); }
LIBCPPRIME_CONSTEXPR std::uint64_t build(std::uint64_t x) const noexcept { return reduce(x % mod_, rs); }
LIBCPPRIME_CONSTEXPR std::uint64_t raw(std::uint64_t x) const noexcept {
Assume(x < mod_);
return reduce(x, rs);
}
LIBCPPRIME_CONSTEXPR std::uint64_t one() const noexcept {
if LIBCPPRIME_IF_CONSTEXPR (Strict) {
Assume(np < mod_);
return np;
} else {
Assume(np < 2 * mod_);
return np;
}
}
LIBCPPRIME_CONSTEXPR std::uint64_t neg(std::uint64_t x) const noexcept {
if LIBCPPRIME_IF_CONSTEXPR (Strict) {
Assume(x < mod_);
return (mod_ - x) * (x != 0);
} else {
Assume(x < 2 * mod_);
return (2 * mod_ - x) * (x != 0);
}
}
LIBCPPRIME_CONSTEXPR std::uint64_t mul(std::uint64_t x, std::uint64_t y) const noexcept {
if LIBCPPRIME_IF_CONSTEXPR (Strict) {
Assume(x < mod_ && y < mod_);
return reduce(x, y);
} else {
Assume(x < 2 * mod_ && y < 2 * mod_);
return reduce(x, y);
}
}
LIBCPPRIME_CONSTEXPR bool same(std::uint64_t x, std::uint64_t y) const noexcept {
if LIBCPPRIME_IF_CONSTEXPR (Strict) {
Assume(x < mod_ && y < mod_);
return x == y;
} else {
Assume(x < 2 * mod_ && y < 2 * mod_);
std::uint64_t tmp = x - y;
return (tmp == 0) || (tmp == mod_) || (tmp == 0 - mod_);
}
}
LIBCPPRIME_CONSTEXPR bool is_zero(std::uint64_t x) const noexcept {
if LIBCPPRIME_IF_CONSTEXPR (Strict) {
Assume(x < mod_);
return x == 0;
} else {
Assume(x < 2 * mod_);
return x == 0 || x == mod_;
}
}
LIBCPPRIME_CONSTEXPR std::uint64_t add(std::uint64_t x, std::uint64_t y) const noexcept {
if LIBCPPRIME_IF_CONSTEXPR (Strict) {
Assume(x < mod_ && y < mod_);
return x + y - (x >= mod_ - y) * mod_;
} else {
Assume(x < 2 * mod_ && y < 2 * mod_);
return x + y - (x >= 2 * mod_ - y) * (2 * mod_);
}
}
LIBCPPRIME_CONSTEXPR std::uint64_t sub(std::uint64_t x, std::uint64_t y) const noexcept {
if LIBCPPRIME_IF_CONSTEXPR (Strict) {
Assume(x < mod_ && y < mod_);
return x - y + (x < y) * mod_;
} else {
Assume(x < 2 * mod_ && y < 2 * mod_);
return x - y + (x < y) * (2 * mod_);
}
}
};
// clang-format off
LIBCPPRIME_CONSTEXPR const std::uint16_t Bases[256] = {
1216,1836,8885,4564,10978,5228,15613,13941,1553,173,3615,3144,10065,9259,233,2362,6244,6431,10863,5920,6408,6841,22124,2290,45597,6935,4835,7652,1051,445,5807,842,1534,22140,1282,1733,347,6311,14081,11157,186,703,9862,15490,1720,17816,10433,49185,2535,9158,2143,2840,664,29074,24924,1035,41482,1065,10189,8417,130,4551,5159,48886,
786,1938,1013,2139,7171,2143,16873,188,5555,42007,1045,3891,2853,23642,148,3585,3027,280,3101,9918,6452,2716,855,990,1925,13557,1063,6916,4965,4380,587,3214,1808,1036,6356,8191,6783,14424,6929,1002,840,422,44215,7753,5799,3415,231,2013,8895,2081,883,3855,5577,876,3574,1925,1192,865,7376,12254,5952,2516,20463,186,
5411,35353,50898,1084,2127,4305,115,7821,1265,16169,1705,1857,24938,220,3650,1057,482,1690,2718,4309,7496,1515,7972,3763,10954,2817,3430,1423,714,6734,328,2581,2580,10047,2797,155,5951,3817,54850,2173,1318,246,1807,2958,2697,337,4871,2439,736,37112,1226,527,7531,5418,7242,2421,16135,7015,8432,2605,5638,5161,11515,14949,
748,5003,9048,4679,1915,7652,9657,660,3054,15469,2910,775,14106,1749,136,2673,61814,5633,1244,2567,4989,1637,1273,11423,7974,7509,6061,531,6608,1088,1627,160,6416,11350,921,306,18117,1238,463,1722,996,3866,6576,6055,130,24080,7331,3922,8632,2706,24108,32374,4237,15302,287,2296,1220,20922,3350,2089,562,11745,163,11951 };
// clang-format on
LIBCPPRIME_CONSTEXPR bool IsPrime32(const std::uint32_t x) noexcept {
if (x < 151321) {
const std::uint32_t a = Divu128(272518712866683587ull % x, 10755835586592736005ull, x).low;
if (a == 0) return false;
if (x < 11881) return GCD32(a, x) == 1;
const std::uint32_t b = Divu128(827936745744686818ull % x, 10132550402535125089ull, x).low;
if (b == 0) return false;
if (x < 39601) return GCD32((a * b) % x, x) == 1;
const std::uint32_t c = Divu128(9647383993136055606ull % x, 17068348107132031867ull, x).low * a * b % x;
if (c == 0) return false;
if (x < 85849) return GCD32(c, x) == 1;
const std::uint32_t d = Divu128(5118528107581154032ull % x, 7251394891134766417ull, x).low;
if (d == 0) return false;
return GCD32(std::uint64_t(c) * d % x, x) == 1;
}
const std::uint32_t h = x * 0xad625b89;
std::uint32_t d = x - 1;
std::uint32_t pw = static_cast<std::uint32_t>(Bases[h >> 24]);
std::uint32_t s = CountrZero(d);
d >>= s;
std::uint32_t cur = pw;
if (d != 1) {
pw = std::uint64_t(pw) * pw % x;
d >>= 1;
while (d != 1) {
std::uint32_t tmp = std::uint64_t(pw) * pw % x;
if (d & 1) cur = std::uint64_t(cur) * pw % x;
pw = tmp;
d >>= 1;
}
cur = std::uint64_t(cur) * pw % x;
}
bool flag = cur == 1 || cur == x - 1;
if (x % 4 == 3) return flag;
if (flag) return true;
while (--s) {
cur = std::uint64_t(cur) * cur % x;
if (cur == x - 1) return true;
}
return false;
}
} // namespace internal
} // namespace cppr
#endif
namespace cppr {
namespace internal {
LIBCPPRIME_CONSTEXPR const std::uint32_t FlagTable10[32] = { 0xa08a28acu, 0x28208a20u, 0x2088288u, 0x800228a2u, 0x20a00a08u, 0x80282088u, 0x800800a2u, 0x8028228u, 0xa20a082u, 0x22880020u, 0x28020800u, 0x88208082u, 0x2022020u, 0x8828028u, 0x8008a202u, 0x20880880u, 0x20000a00u, 0xa082008u, 0x82820802u, 0x800a20u, 0x28208au, 0x20080822u, 0x20808020u, 0x2208088u, 0x20080022u, 0x28a00a00u, 0x8a200080u, 0x8a2000u, 0x808800u, 0x2082202u, 0x80820880u, 0x28220020u };
LIBCPPRIME_CONSTEXPR bool IsPrime10(const std::uint64_t n) noexcept {
return (FlagTable10[n / 32] >> (n % 32)) & 1;
}
LIBCPPRIME_CONSTEXPR std::uint64_t GetLucasBase(const std::uint64_t x) noexcept {
std::uint32_t tmp = x % 5;
if (tmp == 2 || tmp == 3) return 5;
tmp = x % 13;
if (tmp == 2 || tmp == 5 || tmp == 6 || tmp == 7 || tmp == 8 || tmp == 11) return 13;
tmp = x % 17;
if (tmp == 3 || tmp == 5 || tmp == 6 || tmp == 7 || tmp == 10 || tmp == 11 || tmp == 12 || tmp == 14) return 17;
tmp = x % 21;
if (tmp == 2 || tmp == 8 || tmp == 10 || tmp == 11 || tmp == 13 || tmp == 19) return 21;
tmp = x % 29;
if (tmp == 0) return 0;
if (tmp == 2 || tmp == 3 || tmp == 8 || tmp == 10 || tmp == 11 || tmp == 12 || tmp == 14 || tmp == 15 || tmp == 17 || tmp == 18 || tmp == 19 || tmp == 21 || tmp == 26 || tmp == 27) return 29;
if (0x02030213u >> (x & 31) & 1) {
std::int32_t k = 32 - (CountlZero(x - 1) >> 1);
std::uint64_t s = 1ull << k, t = (s + (x >> k)) >> 1;
while (t < s) {
s = t;
t = (s + x / s) >> 1;
}
if (s * s == x) return 0;
}
std::uint64_t Z = 33;
while (Z < x) {
std::uint64_t a = Z, n = x;
bool res = false;
while (a != 0) {
std::int32_t s = CountrZero(a);
a >>= s;
res ^= ((s & 1) && ((n & 0b111) == 3 || (n & 0b111) == 5));
res ^= ((a & 0b11) == 3 && (n & 0b11) == 3);
std::uint64_t tmp = n;
n = a;
a = tmp % n;
}
if (n == 1 && res) break;
Z += 4;
}
if (Z >= x) return 1;
return Z;
}
LIBCPPRIME_CONSTEXPR bool IsPrime64MillerRabin(const std::uint64_t x) noexcept {
MontgomeryModint64Impl<0> mint;
mint.set(x);
const std::int32_t S = CountrZero(x - 1);
const std::uint64_t D = (x - 1) >> S;
const auto one = mint.one(), mone = mint.neg(one);
auto test2 = [=](std::uint64_t base1, std::uint64_t base2) -> bool {
auto a = one, b = one;
auto c = mint.build(base1), d = mint.build(base2);
std::uint64_t ex = D;
while (ex != 1) {
auto e = mint.mul(c, c), f = mint.mul(d, d);
if (ex & 1) a = mint.mul(a, c), b = mint.mul(b, d);
c = e, d = f;
ex >>= 1;
}
a = mint.mul(a, c), b = mint.mul(b, d);
bool res1 = mint.same(a, one) || mint.same(a, mone);
bool res2 = mint.same(b, one) || mint.same(b, mone);
if (!(res1 && res2)) {
for (std::int32_t i = 0; i != S - 1; ++i) {
a = mint.mul(a, a), b = mint.mul(b, b);
res1 |= mint.same(a, mone), res2 |= mint.same(b, mone);
}
if (!res1 || !res2) return false;
}
return true;
};
auto test3 = [=](std::uint64_t base1, std::uint64_t base2, std::uint64_t base3) -> bool {
auto a = one, b = one, c = one;
auto d = mint.build(base1), e = mint.build(base2), f = mint.build(base3);
std::uint64_t ex = D;
while (ex != 1) {
const auto g = mint.mul(d, d), h = mint.mul(e, e), i = mint.mul(f, f);
if (ex & 1) a = mint.mul(a, d), b = mint.mul(b, e), c = mint.mul(c, f);
d = g, e = h, f = i;
ex >>= 1;
}
a = mint.mul(a, d), b = mint.mul(b, e), c = mint.mul(c, f);
bool res1 = mint.same(a, one) || mint.same(a, mone);
bool res2 = mint.same(b, one) || mint.same(b, mone);
bool res3 = mint.same(c, one) || mint.same(c, mone);
if (!(res1 && res2 && res3)) {
for (std::int32_t i = 0; i != S - 1; ++i) {
a = mint.mul(a, a), b = mint.mul(b, b), c = mint.mul(c, c);
res1 |= mint.same(a, mone), res2 |= mint.same(b, mone), res3 |= mint.same(c, mone);
}
if (!res1 || !res2 || !res3) return false;
}
return true;
};
auto test4 = [=](std::uint64_t base1, std::uint64_t base2, std::uint64_t base3, std::uint64_t base4) -> bool {
auto a = one, b = one, c = one, d = one;
auto e = mint.build(base1), f = mint.build(base2), g = mint.build(base3), h = mint.build(base4);
std::uint64_t ex = D;
while (ex != 1) {
auto i = mint.mul(e, e), j = mint.mul(f, f), k = mint.mul(g, g), l = mint.mul(h, h);
if (ex & 1) a = mint.mul(a, e), b = mint.mul(b, f), c = mint.mul(c, g), d = mint.mul(d, h);
e = i, f = j, g = k, h = l;
ex >>= 1;
}
a = mint.mul(a, e), b = mint.mul(b, f), c = mint.mul(c, g), d = mint.mul(d, h);
bool res1 = mint.same(a, one) || mint.same(a, mone);
bool res2 = mint.same(b, one) || mint.same(b, mone);
bool res3 = mint.same(c, one) || mint.same(c, mone);
bool res4 = mint.same(d, one) || mint.same(d, mone);
if (!(res1 && res2 && res3 && res4)) {
for (std::int32_t i = 0; i != S - 1; ++i) {
a = mint.mul(a, a), b = mint.mul(b, b), c = mint.mul(c, c), d = mint.mul(d, d);
res1 |= mint.same(a, mone), res2 |= mint.same(b, mone), res3 |= mint.same(c, mone), res4 |= mint.same(d, mone);
}
if (!res1 || !res2 || !res3 || !res4) return false;
}
return true;
};
// These bases were discovered by Steve Worley and Jim Sinclair.
if (x < 585226005592931977ull) {
if (x < 7999252175582851ull) {
if (x < 350269456337ull) return test3(4230279247111683200ull, 14694767155120705706ull, 16641139526367750375ull);
else if (x < 55245642489451ull) return test2(2ull, 141889084524735ull) && test2(1199124725622454117ull, 11096072698276303650ull);
else return test2(2ull, 4130806001517ull) && test3(149795463772692060ull, 186635894390467037ull, 3967304179347715805ull);
} else return test3(2ull, 123635709730000ull, 9233062284813009ull) && test3(43835965440333360ull, 761179012939631437ull, 1263739024124850375ull);
} else return test3(2ull, 325ull, 9375ull) && test4(28178ull, 450775ull, 9780504ull, 1795265022ull);
}
template<std::int32_t Strict> LIBCPPRIME_CONSTEXPR bool IsPrime64BailliePSW(const std::uint64_t x) noexcept {
MontgomeryModint64Impl<Strict> mint;
mint.set(x);
const std::int32_t S = CountrZero(x - 1);
const std::uint64_t D = (x - 1) >> S;
const auto one = mint.one(), mone = mint.neg(one);
auto miller_rabin_test = [&]() -> bool {
auto a = one, b = mint.raw(2);
std::uint64_t ex = D;
while (ex != 1) {
auto c = mint.mul(b, b);
if (ex & 1) a = mint.mul(a, b);
b = c;
ex >>= 1;
}
a = mint.mul(a, b);
bool flag = mint.same(a, one) || mint.same(a, mone);
if (x % 4 == 3) return flag;
if (flag) return true;
for (std::int32_t i = 0; i != S - 1; ++i) {
a = mint.mul(a, a);
if (mint.same(a, mone)) return true;
}
return false;
};
if (!miller_rabin_test()) return false;
std::uint64_t Z = GetLucasBase(x);
if (Z <= 1) return Z == 1;
const std::uint64_t Q = mint.raw(x - (Z - 1) / 4);
std::uint64_t u = one, v = one, Qn = Q;
std::uint64_t k = (x + 1) << CountlZero(x + 1);
Z = mint.raw(Z);
std::uint64_t t = (x >> 1) + 1;
for (k <<= 1; k; k <<= 1) {
u = mint.mul(u, v);
v = mint.sub(mint.mul(v, v), mint.add(Qn, Qn));
Qn = mint.mul(Qn, Qn);
if (k >> 63) {
std::uint64_t uu = mint.add(u, v);
uu = (uu >> 1) + ((uu & 1) ? t : 0);
v = mint.add(mint.mul(Z, u), v);
v = (v >> 1) + ((v & 1) ? t : 0);
u = uu;
Qn = mint.mul(Qn, Q);
}
}
if (mint.is_zero(u) || mint.is_zero(v)) return true;
std::uint64_t f = (x + 1) & ~x;
for (f >>= 1; f; f >>= 1) {
u = mint.mul(u, v);
v = mint.sub(mint.mul(v, v), mint.add(Qn, Qn));
if (mint.is_zero(v)) return true;
Qn = mint.mul(Qn, Qn);
}
return false;
}
} // namespace internal
LIBCPPRIME_CONSTEXPR bool IsPrimeNoTable(std::uint64_t n) noexcept {
if (n < 1024) return internal::IsPrime10(n);
else {
if ((n & 1) == 0 || 6148914691236517205u >= 12297829382473034411u * n || 3689348814741910323u >= 14757395258967641293u * n || 2635249153387078802u >= 7905747460161236407u * n || 1676976733973595601u >= 3353953467947191203u * n || 1418980313362273201u >= 5675921253449092805u * n || 1085102592571150095u >= 17361641481138401521u * n) return false;
if (n <= 0xffffffff) return internal::IsPrime32(n);
else if (n < (std::uint64_t(1) << 62)) return internal::IsPrime64MillerRabin(n);
else if (n < (std::uint64_t(1) << 63)) return internal::IsPrime64BailliePSW<1>(n);
else return internal::IsPrime64BailliePSW<2>(n);
}
}
} // namespace cppr
#endif
#include <cstring>
#include <unistd.h>
namespace itype {
using i8 = std::int8_t;
using u8 = std::uint8_t;
using i16 = std::int16_t;
using u16 = std::uint16_t;
using i32 = std::int32_t;
using u32 = std::uint32_t;
using i64 = std::int64_t;
using u64 = std::uint64_t;
} // namespace itype
namespace ctype {
using c8 = char;
}
template<itype::u32> constexpr auto InttoStr = [] {
struct {
ctype::c8 table[40004] = {};
} res;
for (itype::u32 i = 0; i != 10000; ++i) {
res.table[4 * i + 0] = (i / 1000 + '0');
res.table[4 * i + 1] = (i / 100 % 10 + '0');
res.table[4 * i + 2] = (i / 10 % 10 + '0');
res.table[4 * i + 3] = (i % 10 + '0');
}
return res;
}();
template<class Stream> constexpr itype::u64 Parseu64(Stream& stream) {
itype::u64 res = 0;
{
itype::u64 v;
std::memcpy(&v, stream.current(), 8);
if (!((v ^= 0x3030303030303030) & 0xf0f0f0f0f0f0f0f0)) {
stream.skip(8);
itype::u64 u;
std::memcpy(&u, stream.current(), 8);
if (!((u ^= 0x3030303030303030) & 0xf0f0f0f0f0f0f0f0)) {
v = (v * 10 + (v >> 8)) & 0x00ff00ff00ff00ff;
u = (u * 10 + (u >> 8)) & 0x00ff00ff00ff00ff;
v = (v * 100 + (v >> 16)) & 0x0000ffff0000ffff;
u = (u * 100 + (u >> 16)) & 0x0000ffff0000ffff;
v = (v * 10000 + (v >> 32)) & 0x00000000ffffffff;
u = (u * 10000 + (u >> 32)) & 0x00000000ffffffff;
res = v * 100000000 + u;
stream.skip(8);
} else {
v = (v * 10 + (v >> 8)) & 0x00ff00ff00ff00ff;
v = (v * 100 + (v >> 16)) & 0x0000ffff0000ffff;
v = (v * 10000 + (v >> 32)) & 0x00000000ffffffff;
res = v;
}
}
}
itype::u64 buf;
std::memcpy(&buf, stream.current(), 8);
{
itype::u32 v = buf;
if (!((v ^= 0x30303030) & 0xf0f0f0f0)) {
buf >>= 32;
v = (v * 10 + (v >> 8)) & 0x00ff00ff;
v = (v * 100 + (v >> 16)) & 0x0000ffff;
res = 10000 * res + v;
stream.skip(4);
}
}
{
itype::u16 v = buf;
if (!((v ^= 0x3030) & 0xf0f0)) {
buf >>= 16;
v = (v * 10 + (v >> 8)) & 0x00ff;
res = 100 * res + v;
stream.skip(2);
}
}
{
const ctype::c8 v = ctype::c8(buf) ^ 0x30;
const bool f = !(v & 0xf0);
res = f ? 10 * res + v : res;
stream.skip(f + 1);
}
return res;
}
template<class Stream> constexpr void Formatu64(Stream& stream, itype::u64 n) {
auto copy1 = [&](itype::u32 x) {
itype::u32 off = (x < 10) + (x < 100) + (x < 1000);
std::memcpy(stream.current(), InttoStr<0>.table + (4 * x + off), 4);
stream.skip(4 - off);
};
auto copy2 = [&](itype::u32 x) {
std::memcpy(stream.current(), InttoStr<0>.table + 4 * x, 4);
stream.skip(4);
};
if (n < 10000000000000000) {
if (n < 1000000000000) {
if (n < 100000000) {
if (n < 10000) copy1(n);
else {
copy1(n / 10000);
copy2(n % 10000);
}
} else {
copy1(n / 100000000);
copy2(n / 10000 % 10000);
copy2(n % 10000);
}
} else {
copy1(n / 1000000000000);
copy2(n / 100000000 % 10000);
copy2(n / 10000 % 10000);
copy2(n % 10000);
}
} else {
copy1(n / 10000000000000000);
copy2(n / 1000000000000 % 10000);
copy2(n / 100000000 % 10000);
copy2(n / 10000 % 10000);
copy2(n % 10000);
}
}
template<itype::u32 Bufsize = (1 << 17)> class BasicReader {
itype::i32 fd = 0;
ctype::c8 buf[Bufsize + 1] = {};
ctype::c8 *cur = buf, *eof = buf;
public:
BasicReader() {}
BasicReader(itype::i32 filehandle) : fd(filehandle) {}
BasicReader(const BasicReader& rhs) {
fd = rhs.fd;
std::memcpy(buf, rhs.buf, rhs.eof - rhs.cur);
cur = buf + (rhs.cur - rhs.buf);
eof = buf + (rhs.cur - rhs.eof);
}
BasicReader& operator=(const BasicReader& rhs) {
fd = rhs.fd;
std::memcpy(buf, rhs.buf, rhs.eof - rhs.cur);
cur = buf + (rhs.cur - rhs.buf);
eof = buf + (rhs.cur - rhs.eof);
return *this;
}
void reload() {
if (eof == buf + Bufsize || eof == cur || [&] {
auto p = cur;
while (*p >= '!') ++p;
return p;
}() == eof) [[likely]] {
itype::u32 rem = eof - cur;
std::memmove(buf, cur, rem);
*(eof = buf + rem + read(fd, buf + rem, Bufsize - rem)) = '\0';
cur = buf;
}
}
void reload(itype::u32 len) {
if (avail() < len) [[unlikely]]
reload();
}
itype::u32 avail() const { return eof - cur; }
const ctype::c8* current() const { return cur; }
void skip(itype::u32 n) { cur += n; }
};
BasicReader rd;
template<itype::u32 Bufsize = (1 << 17)> class BasicWriter {
itype::i32 fd = 1;
ctype::c8 buf[Bufsize + 1] = {};
ctype::c8 *cur = buf, *eof = buf + Bufsize;
public:
BasicWriter() {}
BasicWriter(itype::i32 filehandle) : fd(filehandle) {}
BasicWriter(const BasicWriter& rhs) {
fd = rhs.fd;
std::memcpy(buf, rhs.buf, rhs.cur - rhs.buf);
cur = buf + (rhs.cur - rhs.buf);
}
BasicWriter& operator=(const BasicWriter& rhs) {
fd = rhs.fd;
std::memcpy(buf, rhs.buf, rhs.cur - rhs.buf);
cur = buf + (rhs.cur - rhs.buf);
return *this;
}
void reload() {
[[maybe_unused]] itype::i32 tmp = write(fd, buf, cur - buf);
cur = buf;
}
void reload(itype::u32 len) {
if (eof - cur < len) [[unlikely]]
reload();
}
itype::u32 avail() const { return eof - cur; }
ctype::c8* current() { return cur; }
void skip(itype::u32 n) { cur += n; }
};
BasicWriter wt;
int main() {
rd.reload(20);
int n = Parseu64(rd);
for (int i = 0; i != n; ++i) {
rd.reload(20);
unsigned long long m = Parseu64(rd);
Formatu64(wt, m);
wt.reload(4);
*wt.current() = ' ';
wt.skip(1);
*wt.current() = '0' + cppr::IsPrimeNoTable(m);
wt.skip(1);
*wt.current() = '\n';
wt.skip(1);
}
wt.reload();
}