結果
| 問題 |
No.8030 ミラー・ラビン素数判定法のテスト
|
| ユーザー |
👑  Mizar
|
| 提出日時 | 2022-09-02 16:24:00 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 33 ms / 9,973 ms |
| コード長 | 8,040 bytes |
| コンパイル時間 | 257 ms |
| コンパイル使用メモリ | 33,920 KB |
| 最終ジャッジ日時 | 2025-02-07 00:36:29 |
|
ジャッジサーバーID （参考情報） |
judge3 / judge2 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 10 |
コンパイルメッセージ
main.cpp: In function ‘int main(int, char**)’:
main.cpp:181:10: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
181 | scanf("%d", &n);
| ~~~~~^~~~~~~~~~
main.cpp:184:14: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
184 | scanf("%llu", &x);
| ~~~~~^~~~~~~~~~~~
ソースコード
//#pragma GCC target ("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,avx2,bmi2,lzcnt,tune=native")
//#pragma GCC target ("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native")
//#pragma GCC target ("sse4")
#pragma GCC optimize("O3")
//#pragma GCC optimize ("tree-vectorize")
//#pragma GCC optimize("unroll-loops")
#ifndef NDEBUG
#define NDEBUG
#endif
#include <cassert>
#include <ctime>
#include <cstdio>
#include <cstdbool>
#include <cstdint>
#include <initializer_list>
using u32 = uint32_t;
using u64 = uint64_t;
using u128 = __uint128_t;
using f64 = double;
class U64Mont {
private:
static u64 _ni(u64 n) noexcept { // n * ni == 1 (mod 2**64)
// // n is odd number, n = 2*k+1, n >= 1, n < 2**64, k is non-negative integer, k >= 0, k < 2**63
// ni0 := n; // = 2*k+1 = (1+(2**3)*((k*(k+1)/2)**1))/(2*k+1)
u64 ni = n;
// ni1 := ni0 * (2 - (n * ni0)); // = (1-(2**6)*((k*(k+1)/2)**2))/(2*k+1)
// ni2 := ni1 * (2 - (n * ni1)); // = (1-(2**12)*((k*(k+1)/2)**4))/(2*k+1)
// ni3 := ni2 * (2 - (n * ni2)); // = (1-(2**24)*((k*(k+1)/2)**8))/(2*k+1)
// ni4 := ni3 * (2 - (n * ni3)); // = (1-(2**48)*((k*(k+1)/2)**16))/(2*k+1)
// ni5 := ni4 * (2 - (n * ni4)); // = (1-(2**96)*((k*(k+1)/2)**32))/(2*k+1)
// // (n * ni5) mod 2**64 = ((2*k+1) * ni5) mod 2**64 = 1 mod 2**64
for (int i = 0; i < 5; ++i) {
ni = ni * (2 - n * ni);
}
assert(n * ni == 1); // n * ni == 1 (mod 2**64)
return ni;
}
static u64 _n1(u64 n) noexcept { // == n - 1
return n - 1;
}
static u64 _nh(u64 n) noexcept { // == (n + 1) / 2
return (n >> 1) + 1;
}
static u64 _r(u64 n) noexcept { // == 2**64 (mod n)
return (-n) % n;
}
static u64 _rn(u64 n) noexcept { // == -1 * (2**64) (mod n)
return n - _r(n);
}
static u64 _r2(u64 n) noexcept { // == 2**128 (mod n)
return (u64)((-((u128)n)) % ((u128)n));
}
static u32 _k(u64 n) noexcept { // == trailing_zeros(n - 1)
// https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html#Other-Builtins
return __builtin_ctzll(_n1(n));
}
static u64 _d(u64 n) noexcept { // == (n - 1) >> trailing_zeros(n - 1) // n == 2**k * d + 1
return _n1(n) >> _k(n);
}
public:
const u64 n; // == n
const u64 ni; // n * ni == 1 (mod 2**64)
const u64 n1; // == n - 1
const u64 nh; // == (n + 1) / 2
const u64 r; // == 2**64 (mod n)
const u64 rn; // == -1 * (2**64) (mod n)
const u64 r2; // == 2**128 (mod n)
const u64 d; // == (n - 1) >> trailing_zeros(n - 1) // n == 2**k * d + 1
const u32 k; // == trailing_zeros(n - 1)
U64Mont(u64 n) noexcept
: n(n), ni(_ni(n)), n1(_n1(n)), nh(_nh(n)), r(_r(n)), rn(_rn(n)), r2(_r2(n)), d(_d(n)), k(_k(n))
{ assert((n & 1) == 1); }
u64 add(u64 a, u64 b) const noexcept {
// add(a, b) == a + b (mod n)
assert(a < n);
assert(b < n);
unsigned long long t, u;
// https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html#Integer-Overflow-Builtins
bool fa = __builtin_uaddll_overflow(a, b, &t);
bool fs = __builtin_usubll_overflow(t, n, &u);
return (fa ? u : (fs ? t : u));
}
u64 sub(u64 a, u64 b) const noexcept {
// sub(a, b) == a - b (mod n)
assert(a < n);
assert(b < n);
unsigned long long t;
// https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html#Integer-Overflow-Builtins
bool f = __builtin_usubll_overflow(a, b, &t);
return t + (f ? n : 0);
}
u64 div2(u64 ar) const noexcept {
// div2(ar) == ar / 2 (mod n)
assert(ar < n);
return (ar >> 1) + ((ar & 1) == 0 ? 0 : nh);
}
u64 mrmul(u64 ar, u64 br) const noexcept {
// mrmul(ar, br) == (ar * br) / r (mod n)
// R == 2**64
// gcd(N, R) == 1
// N * ni mod R == 1
// 0 <= ar < N < R
// 0 <= br < N < R
// T := ar * br
// t := floor(T / R) - floor(((T * ni mod R) * N) / R)
// if t < 0 then return t + N else return t
assert(ar < n);
assert(br < n);
u128 t = ((u128)ar) * ((u128)br);
unsigned long long w;
// https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html#Integer-Overflow-Builtins
bool f = __builtin_usubll_overflow((unsigned long long)(t >> 64), (unsigned long long)((((u128)(((u64)t) * ni)) * ((u128)n)) >> 64), &w);
return w + (f ? n : 0);
}
u64 mr(u64 ar) const noexcept {
// mr(ar) == ar / r (mod n)
// R == 2**64
// gcd(N, R) == 1
// N * ni mod R == 1
// 0 <= ar < N < R
// t := floor(ar / R) - floor(((ar * ni mod R) * N) / R)
// if t < 0 then return t + N else return t
assert(ar < n);
u64 v = (u64)((((u128)(ar * ni)) * ((u128)n)) >> 64);
return v == 0 ? 0 : n - v;
}
u64 ar(u64 a) const noexcept {
// ar(a) == a * r (mod n)
assert(a < n);
return mrmul(a, r2);
}
u64 pow(u64 ar, u64 b) const noexcept {
// pow(ar, b) == ((ar / r) ** b) * r (mod n)
assert(ar < n);
if (b == 0) { return r; }
for (; (b & 1) == 0; b >>= 1) { ar = mrmul(ar, ar); }
u64 tr = ar;
for (b >>= 1; b != 0; b >>= 1) {
ar = mrmul(ar, ar);
if ((b & 1) != 0) { tr = mrmul(tr, ar); }
}
return tr;
}
};
U64Mont u64mont_new(u64 n) noexcept { return U64Mont(n); }
u64 u64mont_add(const U64Mont *const mont, u64 ar, u64 br) noexcept { return mont->add(ar, br); }
u64 u64mont_sub(const U64Mont *const mont, u64 ar, u64 br) noexcept { return mont->sub(ar, br); }
u64 u64mont_div2(const U64Mont *const mont, u64 ar) noexcept { return mont->div2(ar); }
u64 u64mont_mrmul(const U64Mont *const mont, u64 ar, u64 br) noexcept { return mont->mrmul(ar, br); }
u64 u64mont_mr(const U64Mont *const mont, u64 a) noexcept { return mont->mr(a); }
u64 u64mont_ar(const U64Mont *const mont, u64 a) noexcept { return mont->ar(a); }
u64 u64mont_pow(const U64Mont *const mont, u64 ar, u64 b) noexcept { return mont->pow(ar, b); }
int ctzll(unsigned long long v) noexcept { return __builtin_ctzll(v); }
int clzll(unsigned long long v) noexcept { return __builtin_clzll(v); }
bool miller_rabin(u64 n) noexcept {
// Deterministic variants of the Miller-Rabin primality test
// http://miller-rabin.appspot.com/
if (n == 2) { return true; }
if (n < 2 || (n & 1) == 0) { return false; }
U64Mont mont(n);
for (const u64 base : {2,325,9375,28178,450775,9780504,1795265022}) {
u64 a = base;
if (a >= n) { a %= n; if (a == 0) { continue; } }
u64 tr = mont.pow(mont.ar(a), mont.d);
if (tr == mont.r) { continue; }
for (u32 j = 1; tr != mont.rn; ++j) {
if (j >= mont.k) { return false; }
tr = mont.mrmul(tr, tr);
}
}
return true;
}
int main(int, char**) {
struct timespec time_monotonic_start, time_process_start, time_monotonic_end, time_process_end;
clock_gettime(CLOCK_MONOTONIC, &time_monotonic_start);
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time_process_start);
int n;
scanf("%d", &n);
for(int i = 0; i < n; ++i) {
unsigned long long x;
scanf("%llu", &x);
printf("%llu %d\n", x, miller_rabin((u64)x) ? 1 : 0);
}
clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &time_process_end);
clock_gettime(CLOCK_MONOTONIC, &time_monotonic_end);
f64 d_sec_monotonic =
(f64)(time_monotonic_end.tv_sec - time_monotonic_start.tv_sec) +
(f64)(time_monotonic_end.tv_nsec - time_monotonic_start.tv_nsec) / (1000 * 1000 * 1000);
f64 d_sec_process =
(f64)(time_process_end.tv_sec - time_process_start.tv_sec) +
(f64)(time_process_end.tv_nsec - time_process_start.tv_nsec) / (1000 * 1000 * 1000);
fprintf(stderr, "time_monotonic:%f, time_process:%f\n", d_sec_monotonic, d_sec_process);
}