結果
| 問題 | No.3030 ミラー・ラビン素数判定法のテスト |
| ユーザー |
👑  Mizar
|
| 提出日時 | 2022-09-20 11:46:08 |
| 言語 | Rust (1.77.0 + proconio) |
| 結果 |
AC
|
| 実行時間 | 16 ms / 9,973 ms |
| コード長 | 54,357 bytes |
| コンパイル時間 | 13,635 ms |
| コンパイル使用メモリ | 399,480 KB |
| 実行使用メモリ | 6,820 KB |
| 最終ジャッジ日時 | 2024-12-22 03:28:53 |
| 合計ジャッジ時間 | 14,629 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge2 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 1 ms
6,816 KB |
| testcase_01 | AC | 1 ms
6,816 KB |
| testcase_02 | AC | 1 ms
6,816 KB |
| testcase_03 | AC | 1 ms
6,816 KB |
| testcase_04 | AC | 12 ms
6,820 KB |
| testcase_05 | AC | 12 ms
6,820 KB |
| testcase_06 | AC | 10 ms
6,820 KB |
| testcase_07 | AC | 10 ms
6,816 KB |
| testcase_08 | AC | 10 ms
6,820 KB |
| testcase_09 | AC | 16 ms
6,816 KB |
ソースコード
// -*- coding:utf-8-unix -*-
// モンゴメリ剰余乗算 (Montgomery modular multiplication)
pub struct Mont<T, BitCounts=u32> {
n: T, // n is odd, and n > 2
ni: T, // n * ni == 1 (mod 2**64 {@T=u64})
nh: T, // == (n + 1) / 2
r: T, // == 2**64 {@T=u64} (mod n)
rn: T, // == -1 * (2**64 {T=u64}) (mod n)
r2: T, // == 2**128 {@T=u64} (mod n)
d: T, // == (n - 1) >> (n - 1).trailing_zeros()
k: BitCounts, // == (n - 1).trailing_zeros()
}
pub trait MontTrait<T> {
// create constant structures for montgomery modular arithmetic
fn new(n: T) -> Self;
// addmod(a, b) == a + b (mod n)
fn add(&self, a: T, b: T) -> T;
// submod(a, b) == a - b (mod n)
fn sub(&self, a: T, b: T) -> T;
// div2mod(ar) == ar / 2 (mod n)
fn div2(&self, ar: T) -> T;
// mrmul(ar, br) == (ar * br) / r (mod n)
// R == 2**64 {@T=u64}
// 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
fn mrmul(&self, ar: T, br: T) -> T;
// mr(ar) == ar / r (mod n)
// R == 2**64 {@T=u64}
// 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
fn mr(&self, ar: T) -> T;
// ar(a) == a * r (mod n)
fn ar(&self, a: T) -> T;
// powir(ar, b) == ((ar / r) ** b) * r (mod n)
fn pow(&self, ar: T, b: T) -> T;
// powir(ar, b) == ((ar / r) ** b) * r (mod n)
fn powodd(&self, ar: T, b: T) -> T;
}
impl MontTrait<u64> for Mont<u64> {
#[inline]
fn new(n: u64) -> Self {
// create constant structures for montgomery modular arithmetic
debug_assert_eq!(n & 1, 1);
// // 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)
let mut ni = n;
// ni1 := ni0 * ((n * ni0 - 3) * (n * ni0) + 3); // = (1+(2**9)*((k*(k+1)/2)**3))/(2*k+1)
// ni2 := ni1 * ((n * ni1 - 3) * (n * ni1) + 3); // = (1+(2**27)*((k*(k+1)/2)**9))/(2*k+1)
// ni3 := ni2 * ((n * ni2 - 3) * (n * ni2) + 3); // = (1+(2**81)*((k*(k+1)/2)**27))/(2*k+1)
// // (n * ni3) mod 2**64 = ((2*k+1) * ni3) mod 2**64 = 1 mod 2**64
for _ in 0..3 {
let t = n.wrapping_mul(ni);
ni = t.wrapping_sub(3).wrapping_mul(t).wrapping_add(3).wrapping_mul(ni);
}
debug_assert_eq!(n.wrapping_mul(ni), 1); // n * ni == 1 (mod 2**64)
let nh = (n >> 1) + 1; // == (n + 1) / 2
let r: u64 = n.wrapping_neg() % n; // == 2**64 (mod n)
let rn: u64 = n - r; // == -(2**64) (mod n)
let r2: u64 = ((n as u128).wrapping_neg() % (n as u128)) as u64; // == 2**128 (mod n)
// n == 2**k * d + 1
let mut d = n - 1;
let k = d.trailing_zeros();
d >>= k;
debug_assert_eq!(Self { n, ni, nh, r, rn, r2, d, k }.mr(r), 1); // r / r == 1 (mod n)
debug_assert_eq!(Self { n, ni, nh, r, rn, r2, d, k }.mrmul(1, r2), r); // r2 / r == r (mod n)
Self { n, ni, nh, r, rn, r2, d, k }
}
#[inline]
fn add(&self, a: u64, b: u64) -> u64 {
// addmod(a, b) == a + b (mod n)
debug_assert!(a < self.n);
debug_assert!(b < self.n);
let (t, fa) = a.overflowing_add(b);
let (u, fs) = t.overflowing_sub(self.n);
if fa { u } else { if fs { t } else { u } }
}
#[inline]
fn sub(&self, a: u64, b: u64) -> u64 {
// submod(a, b) == a - b (mod n)
debug_assert!(a < self.n);
debug_assert!(b < self.n);
let (t, f) = a.overflowing_sub(b);
if f { t.wrapping_add(self.n) } else { t }
}
#[inline]
fn div2(&self, ar: u64) -> u64 {
// div2mod(ar) == ar / 2 (mod n)
debug_assert!(ar < self.n);
let t = ar >> 1;
if (ar & 1) == 0 { t } else { t + self.nh }
}
#[inline]
fn mrmul(&self, ar: u64, br: u64) -> u64 {
// 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
debug_assert!(ar < self.n);
debug_assert!(br < self.n);
let t: u128 = (ar as u128) * (br as u128);
let (t, f) = ((t >> 64) as u64).overflowing_sub((((((t as u64).wrapping_mul(self.ni)) as u128) * (self.n as u128)) >> 64) as u64);
if f { t.wrapping_add(self.n) } else { t }
}
#[inline]
fn mr(&self, ar: u64) -> u64 {
// 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
debug_assert!(ar < self.n);
let (t, f) = (((((ar.wrapping_mul(self.ni)) as u128) * (self.n as u128)) >> 64) as u64).overflowing_neg();
if f { t.wrapping_add(self.n) } else { t }
}
#[inline]
fn ar(&self, a: u64) -> u64 {
// ar(a) == a * r (mod n)
debug_assert!(a < self.n);
self.mrmul(a, self.r2)
}
#[inline]
fn pow(&self, mut ar: u64, mut b: u64) -> u64 {
// == ((ar / r) ** b) * r (mod n)
debug_assert!(ar < self.n);
let mut t = if (b & 1) == 0 { self.r } else { ar };
loop {
b >>= 1;
if b == 0 { return t; }
ar = self.mrmul(ar, ar);
if (b & 1) != 0 { t = self.mrmul(t, ar); }
}
}
#[inline]
fn powodd(&self, mut ar: u64, mut b: u64) -> u64 {
// == ((ar / r) ** b) * r (mod n)
debug_assert!(ar < self.n);
debug_assert_eq!(b & 1, 1);
let mut t = ar;
loop {
b >>= 1;
if b == 0 { return t; }
ar = self.mrmul(ar, ar);
if (b & 1) != 0 { t = self.mrmul(t, ar); }
}
}
}
impl MontTrait<u32> for Mont<u32> {
#[inline]
fn new(n: u32) -> Self {
// create constant structures for montgomery modular arithmetic
debug_assert_eq!(n & 1, 1);
// // 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**2)*((k*(k+1))**1))/(2*k+1)
let mut ni = n;
// ni1 := ni0 * (2 - (n * ni0)); // = (1-(2**4)*((k*(k+1))**2))/(2*k+1)
// ni2 := ni1 * (2 - (n * ni1)); // = (1-(2**8)*((k*(k+1))**4))/(2*k+1)
// ni3 := ni2 * (2 - (n * ni2)); // = (1-(2**16)*((k*(k+1))**8))/(2*k+1)
// ni4 := ni3 * (2 - (n * ni3)); // = (1-(2**32)*((k*(k+1))**16))/(2*k+1)
// // (n * ni4) mod 2**32 = ((2*k+1) * ni5) mod 2**32 = 1 mod 2**32
for _ in 0..4 {
ni = ni.wrapping_mul(2u32.wrapping_sub(n.wrapping_mul(ni)));
}
debug_assert_eq!(n.wrapping_mul(ni), 1); // n * ni == 1 (mod 2**32)
let nh = (n >> 1) + 1; // == (n + 1) / 2
let r: u32 = n.wrapping_neg() % n; // == 2**32 (mod n)
let rn: u32 = n - r; // == -(2**32) (mod n)
let r2: u32 = ((n as u64).wrapping_neg() % (n as u64)) as u32; // == 2**64 (mod n)
// n == 2**k * d + 1
let mut d = n - 1;
let k = d.trailing_zeros();
d >>= k;
debug_assert_eq!(Self { n, ni, nh, r, rn, r2, d, k }.mr(r), 1); // r / r == 1 (mod n)
debug_assert_eq!(Self { n, ni, nh, r, rn, r2, d, k }.mrmul(1, r2), r); // r2 / r == r (mod n)
Self { n, ni, nh, r, rn, r2, d, k }
}
#[inline]
fn add(&self, a: u32, b: u32) -> u32 {
// addmod(a, b) == a + b (mod n)
debug_assert!(a < self.n);
debug_assert!(b < self.n);
let (t, fa) = a.overflowing_add(b);
let (u, fs) = t.overflowing_sub(self.n);
if fa || !fs { u } else { t }
}
#[inline]
fn sub(&self, a: u32, b: u32) -> u32 {
// submod(a, b) == a - b (mod n)
debug_assert!(a < self.n);
debug_assert!(b < self.n);
let (t, f) = a.overflowing_sub(b);
if f { t.wrapping_add(self.n) } else { t }
}
#[inline]
fn div2(&self, ar: u32) -> u32 {
// div2mod(ar) == ar / 2 (mod n)
debug_assert!(ar < self.n);
if (ar & 1) == 0 {
ar >> 1
} else {
(ar >> 1) + self.nh
}
}
#[inline]
fn mrmul(&self, ar: u32, br: u32) -> u32 {
// mrmul(ar, br) == (ar * br) / r (mod n)
// R == 2**32
// 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
debug_assert!(ar < self.n);
debug_assert!(br < self.n);
let t: u64 = (ar as u64) * (br as u64);
let (t, f) = ((t >> 32) as u32).overflowing_sub((((((t as u32).wrapping_mul(self.ni)) as u64) * (self.n as u64)) >> 32) as u32);
if f { t.wrapping_add(self.n) } else { t }
}
#[inline]
fn mr(&self, ar: u32) -> u32 {
// mr(ar) == ar / r (mod n)
// R == 2**32
// 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
debug_assert!(ar < self.n);
let (t, f) = (((((ar.wrapping_mul(self.ni)) as u64) * (self.n as u64)) >> 32) as u32).overflowing_neg();
if f { t.wrapping_add(self.n) } else { t }
}
#[inline]
fn ar(&self, a: u32) -> u32 {
// ar(a) == a * r (mod n)
debug_assert!(a < self.n);
self.mrmul(a, self.r2)
}
#[inline]
fn pow(&self, mut ar: u32, mut b: u32) -> u32 {
// powir(ar, b) == ((ar / r) ** b) * r (mod n)
debug_assert!(ar < self.n);
let mut t = if (b & 1) == 0 { self.r } else { ar };
loop {
b >>= 1;
if b == 0 { return t; }
ar = self.mrmul(ar, ar);
if (b & 1) != 0 { t = self.mrmul(t, ar); }
}
}
#[inline]
fn powodd(&self, mut ar: u32, mut b: u32) -> u32 {
// powir(ar, b) == ((ar / r) ** b) * r (mod n)
debug_assert!(ar < self.n);
debug_assert_eq!(b & 1, 1);
let mut t = ar;
loop {
b >>= 1;
if b == 0 { return t; }
ar = self.mrmul(ar, ar);
if (b & 1) != 0 { t = self.mrmul(t, ar); }
}
}
}
#[inline]
fn primetest_u64_miller_sub(mont: &Mont<u64>, mut base: u64) -> bool {
// Miller-Rabin primality test
let (n, r, rn, d, k) = (mont.n, mont.r, mont.rn, mont.d, mont.k);
if base > n { base %= n; if base == 0 { return true; } }
let mut tr = mont.powodd(mont.ar(base), d);
if tr == r || tr == rn { return true; }
for _ in 1..k { tr = mont.mrmul(tr, tr); if tr == rn { return true; } }
false
}
#[inline]
pub fn primetest_u32_miller_sub(mont: &Mont<u32>, mut base: u32) -> bool {
// Miller-Rabin primality test
let (n, r, rn, d, k) = (mont.n, mont.r, mont.rn, mont.d, mont.k);
if base > n { base %= n; if base == 0 { return true; } }
let mut tr = mont.powodd(mont.ar(base), d);
if tr == r || tr == rn { return true; }
for _ in 1..k { tr = mont.mrmul(tr, tr); if tr == rn { return true; } }
false
}
#[inline]
fn primetest_u64_base2_sub(mont: &Mont<u64>) -> bool {
// Miller-Rabin primality test (base 2)
// strong pseudoprimes to base 2 ( https://oeis.org/A001262 ): 2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,...
let (r, rn, d, k) = (mont.r, mont.rn, mont.d, mont.k);
let mut tr = mont.powodd(mont.add(r, r), d);
if tr == r || tr == rn { return true; }
for _ in 1..k { tr = mont.mrmul(tr, tr); if tr == rn { return true; } }
false
}
#[inline]
pub fn primetest_u32_base2_sub(mont: &Mont<u32>) -> bool {
// Miller-Rabin primality test (base 2)
// strong pseudoprimes to base 2 ( https://oeis.org/A001262 ): 2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,...
let (r, rn, d, k) = (mont.r, mont.rn, mont.d, mont.k);
let mut tr = mont.powodd(mont.add(r, r), d);
if tr == r || tr == rn { return true; }
for _ in 1..k { tr = mont.mrmul(tr, tr); if tr == rn { return true; } }
false
}
// 64bit整数平方根(固定ループ回数) -> (floor(sqrt(iv)), remain)
#[inline]
pub fn isqrt64f(iv: u64) -> (u64, u64) { isqrt64i(iv, 0) }
// 64bit整数平方根(可変ループ回数) -> (floor(sqrt(iv)), remain)
#[inline]
pub fn isqrt64d(iv: u64) -> (u64, u64) { isqrt64i(iv, iv.leading_zeros()) }
// 64bit整数平方根(lz:ケチるループ回数*2+(0~1)、内部実装) -> (floor(sqrt(iv)), remain)
#[inline]
fn isqrt64i(iv: u64, lz: u32) -> (u64, u64) {
let n = (64 >> 1) - (lz >> 1);
let s = (lz >> 1) << 1;
let t = n << 1;
let (mut a, mut b, c, d, e) = (
iv as u128,
0x0000_0000_0000_0000_4000_0000_0000_0000 >> s,
0xffff_ffff_ffff_fffe_0000_0000_0000_0000 >> s,
0x0000_0000_0000_0001_0000_0000_0000_0000 >> s,
0x0000_0000_0000_0000_ffff_ffff_ffff_ffff >> s,
);
for _ in 0..n {
if a >= b {
a -= b;
b = ((b + b) & c) + (b & e) + d;
} else {
b = ((b + b) & c) + (b & e);
}
a <<= 2;
}
((b >> t) as u64, (a >> t) as u64)
}
#[allow(unused)]
#[inline]
// ニュートン法による整数平方根
fn isqrt64n(x: u64) -> u64 {
if x <= 1 { return x; }
let k = 32 - ((x - 1).leading_zeros() >> 1);
let mut s = 1u64 << k; // s = 2**k
let mut t = (s + (x >> k)) >> 1; // t = (s + x/s)/2
// whileループ回数=除算回数は最大6回を想定
// s > floor(sqrt(x)) :: floor(sqrt(x)) <= t < s
// s == floor(sqrt(x)) :: floor(sqrt(x)) <= t <= floor(sqrt(x)) + 1
while t < s {
s = t;
t = (s + (x / s)) >> 1;
}
s
}
#[allow(unused)]
#[inline]
// 平方数判定 (true: 平方数である, false: 平方数ではない)
fn issq_u64_n(x: u64) -> bool {
let sqrt = isqrt64n(x);
sqrt * sqrt == x
}
#[allow(unused)]
#[inline]
// 擬平方数判定 mod64 12/64 0.1875 (true: 平方数かもしれない, false: 平方数ではない)
fn issq_u64_mod64(x: u64) -> bool {
(0x202021202030213u64 >> (x & 63)) & 1 == 1
}
#[allow(unused)]
#[inline]
// 擬平方数判定 mod4095 336/4095 0.0821 (true: 平方数かもしれない, false: 平方数ではない)
fn issq_u64_mod4095(x: u64) -> bool {
const SQTABLE_MOD4095: [u64; 64] = [0x2001002010213,0x4200001008028001,0x20000010004,0x80200082010,0x1800008200044029,0x120080000010,0x2200000080410400,0x8100041000200800,0x800004000020100,0x402000400082201,0x9004000040,0x800002000880,0x18002000012000,0x801208,0x26100000804010,0x80000080000002,0x108040040101045,0x20c00004000102,0x400000100c0010,0x1300000040208,0x804000020010000,0x1008402002400080,0x201001000200040,0x4402000000806000,0x10402000000,0x1040008001200801,0x4080000000020400,0x10083080000002,0x8220140000040000,0x800084020100000,0x80010400010000,0x1200020108008060,0x180000000,0x400002400000018,0x4241000200,0x100800000000,0x10201008400483,0xc008000208201000,0x800420000100,0x2010002000410,0x28041000000,0x4010080000024,0x400480010010080,0x200040028000008,0x100810084020,0x20c0401000080000,0x1000240000220000,0x4000020800,0x410000000480000,0x8004008000804201,0x806020000104000,0x2080002000211000,0x1001008001000,0x20000010024000,0x480200002040000,0x48200044008000,0x100000000010080,0x80090400042,0x41040200800200,0x4000020100110,0x2000400082200010,0x1008200000000040,0x2004800002,0x2002010000080];
let p = (x % 4095) as usize;
(SQTABLE_MOD4095[p >> 6] >> (p & 63)) & 1 == 1
}
#[allow(unused)]
#[inline]
// 擬平方数判定 mod4096 684/4096 0.1670 (true: 平方数かもしれない, false: 平方数ではない)
fn issq_u64_mod4096(x: u64) -> bool {
const SQTABLE_MOD4096: [u64; 64] = [0x202021202030213,0x202021202020213,0x202021202030212,0x202021202020212,0x202021202030213,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020213,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030213,0x202021202020213,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020213,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020213,0x202021202030212,0x202021202020212,0x202021202030213,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020213,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020213,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020213,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212,0x202021202030212,0x202021202020212];
let p = (x % 4096) as usize;
(SQTABLE_MOD4096[p >> 6] >> (p & 63)) & 1 == 1
}
#[allow(unused)]
#[inline]
// 擬平方数判定 mod7429 1080/7429 0.1454 (true: 平方数かもしれない, false: 平方数ではない)
fn issq_u64_mod7429(x: u64) -> bool {
const SQTABLE_MOD7429: [u64; 117] = [0x82801806010213,0xac8013020a22001,0x1800020200011400,0x6000840210200010,0x812400220008001,0x14024230100018,0x4a40088100118080,0x402098101100b000,0x4010100000200c4,0x20a00400002080b,0x200010804004d0,0xa20110080a040400,0x6408048458000,0x9080008100224,0x108001000000601a,0x510a00800040682,0x9000c04000002801,0x100404020c0006,0x100402040800010,0x26801440040102d0,0xc0802080211400,0x88530000a2014008,0x20221001640098,0x842000004010100,0x1000420030000280,0x4608220400092001,0x40081040180a8001,0x980101001200460,0x1211080020044c00,0x20010001c080a400,0x3021182800050020,0x500020840140008,0x2408800451000922,0x8008008000422010,0xf0204a200,0x8000416040400009,0x404140002a01011,0x4042248006010,0x440040801000082,0x4007820081010,0x40000b0013410a60,0x5020a02001000084,0x222000440098880,0x20000000010040a,0x80001c080688004,0x220181088080100,0x54080a8041042,0x1018003000064021,0x1000020004c02008,0x6100160009040521,0x182c00010408200,0xa080400008114,0x8c0449002902104,0x408100002010080,0x8000500042200008,0x8402000482009003,0x4154200a000110b0,0x600200002000a00,0x2000002100000200,0x5000090004004,0xa001301002001c,0xa000400008c580,0x3210440028885002,0x200018100410022,0x34100608004001,0x4800080042000e0,0x5500088001020c02,0x90032000a442a000,0x200c04000c0110,0xc0c09000001108,0x808410420a00300,0x800440200016008,0x8048000200058002,0x100604010600008,0x306002202108840,0x6040480140002880,0x4400801051010080,0x10800ae0405,0x3010083400024,0x4830291010000000,0x201010a08008400,0x300400800c504209,0x1440080880700020,0x880001b022220,0x1802000a0001000,0x10100850000a2003,0x8088040422420040,0x120050042b810451,0x84000058011808,0x402000501100102,0x805000900061000c,0x40200004100080,0x2c2010c00220c,0x2200104c0804841,0x400a202800840800,0x200014c001880030,0x2800001030800600,0x170002a0004460,0x2120081600000204,0x1004840204080,0x410060004000481,0x8204484409021,0x3c081200022100,0x8980011002221004,0x10800a4000020001,0x85340082001010,0x8001042003000008,0x5e0022800441002,0xc000d0001908800,0x8000020100100004,0x21a006000140c0,0x20a104118280848,0x702118080204000,0x6010240084004100,0x202900040000002,0x45001880050205,0x0];
let p = (x % 7429) as usize;
(SQTABLE_MOD7429[p >> 6] >> (p & 63)) & 1 == 1
}
#[allow(unused)]
#[inline]
// 擬平方数判定 mod33263 4560/33263 0.1371 (true: 平方数かもしれない, false: 平方数ではない)
fn issq_u64_mod33263(x: u64) -> bool {
const SQTABLE_MOD33263: [u64; 520] = [0xc002005212010293,0x201881000024089,0xe20001210410,0x9100028608000014,0x1000a8250001441,0x6220020020c00050,0x120c10080c812004,0x402201000620001,0x80490600038208,0xc022400108122001,0x5118801008000020,0x20200a000000,0x200288800210300,0x300000008600,0x404002d010080118,0x100d44810000004e,0x40701000108a001,0x890180400000a002,0x4a00141290000014,0x4184000206,0xa0a0020085090020,0x20b0000013004000,0x8405004080080,0x2112021800100002,0x8d0408440204104,0x40201120900a0003,0x1064000021800,0x201821000041100,0x288040080400c01,0x40801500072,0x38000040010088,0x400c838000040410,0x5002201680021000,0x48000a0c08000053,0x120120f00080118,0x488040018020040,0x2050008000808001,0x104204010010,0x700002004020801,0x1104800080900010,0x8008420008304000,0xe188080088000,0x1040800422016608,0x30100810c811080,0x2040000240200040,0x641c1000002080,0x8014042000188220,0x4000006040028002,0x2000803902008844,0x400200440100225,0x292100014041008,0x4080025024480000,0x1000011068004008,0xc24402000800022,0x13800002840002,0x20040000c00804,0x8000200100231040,0x889060004089800,0xc20001011412200,0x108e08001004480,0x80240001041912,0x4080408c4,0x410001c812040003,0x1002420001121,0x1490a00008208001,0x2400000102401400,0x2800080140000442,0x4220012020000111,0x109c800004a00094,0x20620000a224000,0x1004080881400200,0x1048300300044401,0x2000000008010100,0x1041000008020040,0x40101000e210890,0x104187000286680,0x24021080222110,0x100801200400020c,0x20060000c1080009,0x2020800002002040,0x8000a04106000,0x18020d0120001089,0xe40010800004a0,0x1030000012120000,0x816800400ac80020,0xa201801902000008,0x80048440080004,0x8a30200050510200,0xa002000700e0000,0x80418000043000,0x120608001018800,0x8000080020041842,0x420020a00600050,0x8204880010204,0x2002201000011008,0xc890400d0001c010,0xa100008504000500,0x2180000200100800,0x940400052608100,0x10e811290004,0x242080400035,0x1081000000801014,0x404b400380720400,0x5000c0008388000,0x5142400006000,0x4000022c090,0x2040011003008000,0x804024480000001,0x1020100080448,0x4800040000621c0,0x3000003040010846,0x400401444000001,0xa00110080000803,0x1060820089880a00,0x10400200028,0x80010100890c0,0x24a000403912012,0x2820008141c94000,0x800212044621000,0x400000068480,0x2000a0018001220,0x4128002000408001,0x200820000240,0x402000000311600,0x8800210800090600,0x420000c224020108,0x1080081118004028,0x104200000401102,0x1080080000c40,0x50000080a2400010,0x5000c204010114,0x418304008e480050,0x40a4000023018011,0x8000104000804482,0x2000541000009100,0x94010000a040a40,0x440004102211002,0x40020000003400,0x1025000401191,0x42100100800,0x400000cc1000042,0x21102062008812,0x48050000806824,0x200010412203008,0x10800200c0400083,0x28200410000080,0x628080018040082,0x90000045870002,0x8022c00600010800,0x8a04890000204243,0x2205020021080010,0x68040000c040204,0x8128000008489,0x98000804252,0x2022000608100800,0x80020113000000a0,0x40080400030020,0x1000082010018040,0x9000000222404400,0x600100a30a000504,0x4102408842200440,0x80100024190080,0x1040010000c0000,0x100448000000120c,0x1020004000408006,0x14068000a22c1000,0x1160012800040,0x4004c4804025310,0x21102b1000005040,0x834001c00200001,0x1000400000080106,0x14044010e0200,0x14000120208c0,0x400110890000005,0x8802044430002003,0x600801108c00483,0x888000041620000,0x100000408108060,0x200020442000032,0xa014000001a08008,0x21c200090222000,0x2004624000000080,0x81190200030400,0x2a0044005102500,0x9080081000940800,0x2022448211000,0xc014880184600,0x3040a000a0010100,0x1080022088000018,0x210400824008241c,0x4040080003100000,0x902008040800000,0x1800102201200180,0xa0001088440044,0x1064000002191201,0x7000000080c400,0x4000c81020100101,0x502242000012040,0x8010400804800180,0x40010242080004,0x10c4000803820,0x5452000280928,0x848100048008a042,0x820041a324100f0,0x408600040004000,0x4890050005041800,0xc020021001010000,0x10000c412004220,0x300000400089020,0x198100000260000,0x2042008100210008,0x8801200004020002,0x282004000100c10,0x282800210440020,0x300020400508000,0x28312000104900,0x240884000004041c,0x4000201292000100,0x70000002c0000040,0x80280000014011,0x8004180008006502,0x20000820d1020006,0x20012046000200,0xc44804004380000,0x221000001080000,0x88400220180211,0x120200a0006085,0x1004488022000000,0x1803100000880800,0x80400d0000821000,0x444110000020e0,0x1600100830a1,0x8018041220000022,0x2001408109000009,0x20c42040012000,0x710880108000,0xcb00014232040002,0x4000128800020,0x1510204010600200,0x406004102c09000,0x1100900810030,0x8022448300000100,0xe018090084404000,0x88c42001400c000,0x108011000704,0x344400008100,0xc0208000000300,0x200200080820000c,0x6200000000011,0x380f000844180,0x4800402191205000,0x210000094840000a,0x10201c2100003,0x2050028806848400,0x408000a00102100,0x202888000041,0x10e0080c00820001,0x460100220160002,0x1010080088040800,0x1832412010808,0x8601010804000882,0x40201041802200,0x8000880009,0x1c03a200000400,0x4002418401020050,0x8900000a40000800,0x2400100600800029,0x90140080042204,0x44020100810281,0x810004040024020,0x820040210c100008,0x382801400340020,0x8882600050418100,0x204210010580000,0x10080c4020044000,0x1008210011304,0x4180000086402440,0x201c1020100810,0x120c4040200200,0x4800021080004110,0x2000003004040000,0x9000a00581000022,0x2200120800081000,0x4489080000020140,0x21020128a0000001,0x84400a824000104,0x613002000040815,0x49c00083421200,0x20100200000,0x81400e000a821,0x440040002802000,0x610200008824000,0x84622241802000,0x20120c00020c90,0x8204000220200462,0xa001002409008041,0xd04010010000,0x448511000100a09,0x100000a00,0x868005002c000221,0x2118011010004088,0xa4e480402501002,0x20b100801110000,0xa040240700400108,0x10000180010200,0x10064400041a0010,0x2000120c018000,0x80a000108500,0x200c3000003300,0x20810840440288,0x88042241021c0006,0x3088000845048,0x400402000001001,0x1500080120003002,0xc00088082800044,0x724000c0008110,0x18000009080000,0x4003802240000,0x1008000c80080004,0x52204220100082,0x5184880680c0000,0x12000070000,0x8000a08808021810,0x24000205400,0x4001000010601021,0x210040000208040,0x4114100008089,0xa001100840028830,0x6600010618102008,0x5080080680040080,0xc0020000010000,0x28204011304100,0x8040086400048000,0x1021006010418,0x84440000068000a,0x25200000392000,0x400401c4000,0x1a121100b004,0x28804003400048,0x9022000101101000,0x20020c0000001048,0x801200040280,0x10020820008110c1,0x9810003820004423,0x200000222000104,0x1000c0002001004,0x200410810040100,0xa80008034082004,0x2041800000060,0x121400001c90009,0x4010000200e42402,0x5000109080110820,0x440020000004240,0x8111000008009110,0x9440000020a200,0x10d000800480,0x114006000a8a22,0x4480408408020300,0x2100000110000084,0x84020041c000020,0x1202104000200804,0x6440840180400621,0x2110004018000008,0xc00010210000a,0xc1100007202010,0x16000400a0c100,0x68041820c0002202,0x1088000800040880,0x4a0800005050020,0x2a0008005108040,0x2003000001210,0x4000928001008c8,0x4001080000500,0x12000022e0840005,0x108a0081005000,0x2204020402003008,0x8809008040001804,0x842000060c00041,0x810808008010080,0x4034042004020010,0x5120800a00038000,0xc000000228000401,0x4504000000281121,0x1100b00030850200,0x4021c102020000,0x80494040804260,0x2004c010002929,0x9014000704200800,0x100000000000,0x20004800018044,0x600824008000c038,0x4090000116000200,0x7000040184080402,0x181201021100080,0x9000804080c8080,0x200010010122000a,0x840800000400906,0x801a00100301003,0x2283040011145400,0x4008a0000400000,0x262080116989,0x90c0008801804020,0x412010200100000,0x24082040200,0x41a024000701c8,0x8480000030000462,0x8220000440000002,0x401000058000403,0x800412004600400,0x1000800009110010,0x700000208244031,0x10c000221089019,0x1008404020c22004,0x200e108000111800,0x800310400820048,0x458108020200,0x1140800008000a0,0x840000000420020,0x20101a000108f44,0x440000a400042301,0x1100000010010080,0x484000024410201c,0x820d0020002800,0x4044000880222,0x1800004000200000,0x8200011024000c0,0x8042020144210022,0x20b20400010cc008,0x4008002004406201,0x3000240800004080,0x80401001020120,0x4000500002008900,0x60010001041000,0x862c00001800,0x808029000000487,0x4242000240000021,0xc1088011900040a,0x24220040002112,0xc000028001020c04,0x4240010812200000,0x4102a09800150,0x180200800060041,0x4c000020130000,0x8011004000828002,0x4402400008000801,0x104810000004808,0x340600008230100,0x12000284204,0x188080002419,0x120002800000060c,0x2007008948082000,0x141220000080,0x80014140800a220,0x52001c0110003280,0x106080504c400002,0x2001200023100208,0x2220800021c1040,0x82800804140000,0x80020014000010e,0x4080040004042,0x1803400011088800,0x1880000443000,0x2000600020108c0,0x8008048800001062,0x3a040000410050,0x808001110054001,0x4000a08004800100,0x621c01009030,0x50008846a004000,0x6500002119001001,0x8b4400001,0x41810c000b00080,0x4081000200024820,0x100010401100018,0x4018000180104080,0xc100042208020,0x1006210010004000,0x284002000310,0x1c0222200880480,0x44040000222400,0x1001004108800,0x8003106448000282,0x800507000022082,0x20400010c8000a00,0x90800004a424020d,0x2000820001080400,0x8004020524200,0x902000108901080,0x8040001445a20020,0x870c00200048004,0x80300800000a3200,0x408000a0000,0x201019004008404,0x4220240402112041,0x20200020c40000,0x30240440010,0x2019800040024,0xc200001021,0xc410000520209040,0x1600100820400004,0x448028100021000,0x6810010004008050,0x480002035002010,0x4100001084004024,0x1c0400400129000,0x930000060c004,0x64001080a2050001,0x41222018000400,0x1008108200012,0xc0000203204000,0x8100c4400280420,0x20406210000810a,0x2080110841048088,0x3000000002000000,0x2000100800802000,0x80080200040081,0x2400040028010082,0x10b0009441020000,0x422402012000002,0x22000];
let p = (x % 33263) as usize;
(SQTABLE_MOD33263[p >> 6] >> (p & 63)) & 1 == 1
}
#[allow(unused)]
#[inline]
// 平方数判定 (true: 平方数である, false: 平方数ではない)
fn issq_u64(x: u64) -> bool {
issq_u64_mod64(x) &&
issq_u64_mod4095(x) &&
issq_u64_n(x)
}
// Jacobi symbol: ヤコビ記号
#[inline]
fn jacobi(a: i64, mut n: u64) -> i32 {
let (mut a, mut j): (u64, i32) = if a >= 0 { (a as u64, 1) } else if (n & 3) == 3 { (a.wrapping_neg() as u64, -1) } else { (a.wrapping_neg() as u64, 1) };
while a > 0 {
let ba = a.trailing_zeros();
a >>= ba;
if ((n & 7) == 3 || (n & 7) == 5) && (ba & 1) != 0 { j = -j; }
if (a & n & 3) == 3 { j = -j; }
let t = n; n = a; a = t; a %= n;
if a > (n >> 1) {
a = n - a;
if (n & 3) == 3 { j = -j; }
}
}
if n == 1 { j } else { 0 }
}
fn primetest_u64_lucas_sub(mont: &Mont<u64>) -> bool {
// Lucas primality test
// strong Lucas pseudoprimes ( https://oeis.org/A217255 ): 5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,...
let n = mont.n;
let mut d: i64 = 5;
for i in 0.. { // d = 5, -7, 11, -13, ... から jacobi(d, n) == -1 となる値を探す
assert!(i < 64); // n が平方数でなければ jacobi(d, n) == -1 を満たす値はすぐに見つかるはず
// 5 <= d.abs() < n かつ jacobi(d, n) == 0 ならば gcd(d, n) != 1 、つまり n は素数ではない
// n == u64::MAX の場合、 gcd(5, u64::MAX) == 5 より jacobi(5, u64::MAX) == 0 とるので
// (u64::MAX + 1) のオーバーフローは気にしなくて良い
match jacobi(d, n) {
-1 => break,
0 => if ((d.abs()) as u64) < n { return false; },
_ => {},
}
// n が平方数の場合、 jacobi(d, n) == -1 となる値は永遠に見つからないのである程度ループしたら平方数チェックを行う
//if i == 32 && isqrt64f(n).1 == 0 { return false; }
if i == 8 && issq_u64_n(n) { return false; }
if (i & 1) == 0 { d = -(d + 2); } else { d = 2 - d; }
}
let n1p = n.wrapping_add(1);
let qm = mont.ar(if d < 0 {(((1 - d) as u64) >> 2) % n} else {n - (((d - 1) as u64) >> 2) % n});
let mut k = n1p << n1p.leading_zeros();
let mut um = mont.r;
let mut vm = mont.r;
let mut qn = qm;
// n > 1, jacobi(d,n) == -1 より gcd(d,n) == 1, d % n != 0
let dm: u64 = mont.ar(if d < 0 { n - (((-d) as u64) % n) } else { (d as u64) % n });
k <<= 1;
while k > 0 {
um = mont.mrmul(um, vm);
vm = mont.sub(mont.mrmul(vm, vm), mont.add(qn, qn));
qn = mont.mrmul(qn, qn);
if (k >> 63) != 0 {
let mut uu = mont.add(um, vm);
uu = mont.div2(uu);
vm = mont.add(mont.mrmul(dm, um), vm);
vm = mont.div2(vm);
um = uu;
qn = mont.mrmul(qn, qm);
}
k <<= 1;
}
if um == 0 || vm == 0 {
return true;
}
let mut x = n1p & (!n);
x >>= 1;
while x > 0 {
um = mont.mrmul(um, vm);
vm = mont.sub(mont.mrmul(vm, vm), mont.add(qn, qn));
if vm == 0 {
return true;
}
qn = mont.mrmul(qn, qn);
x >>= 1;
}
false
}
// Baillie–PSW primarity test
#[inline]
fn primetest_u64_bpsw_sub(mont: &Mont<u64>) -> bool {
// Miller-Rabin primality test (base 2)
// strong pseudoprimes to base 2 ( https://oeis.org/A001262 ): 2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,...
primetest_u64_base2_sub(mont) &&
// Lucas primality test
// strong Lucas pseudoprimes ( https://oeis.org/A217255 ): 5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,...
primetest_u64_lucas_sub(mont)
}
// Baillie–PSW primarity test
pub fn primetest_u64_bpsw(n: u64) -> bool {
if n == 2 { return true; }
if n == 1 || (n & 1) == 0 { return false; }
let mont = Mont::<u64>::new(n);
primetest_u64_bpsw_sub(&mont)
}
// Deterministic variants of the Miller-Rabin primality test
// http://miller-rabin.appspot.com/
const MILLER_U64_1BASES: [u64; 1] = [
9345883071009581737,
];
const MILLER_U64_2BASES: [u64; 2] = [
336781006125,
9639812373923155,
];
const MILLER_U64_3BASES: [u64; 3] = [
4230279247111683200,
14694767155120705706,
16641139526367750375,
];
const MILLER_U64_4BASES: [u64; 4] = [
2,
141889084524735,
1199124725622454117,
11096072698276303650,
];
const MILLER_U64_5BASES: [u64; 5] = [
2,
4130806001517,
149795463772692060,
186635894390467037,
3967304179347715805,
];
const MILLER_U64_6BASES: [u64; 6] = [
2,
123635709730000,
9233062284813009,
43835965440333360,
761179012939631437,
1263739024124850375,
];
const MILLER_U64_7BASES: [u64; 7] = [
2,
325,
9375,
28178,
450775,
9780504,
1795265022,
];
#[inline]
pub fn primetest_u64_miller_1base_sub(mont: &Mont<u64>) -> bool {
assert!(mont.n < 341531);
MILLER_U64_1BASES.iter().all(|&base| primetest_u64_miller_sub(&mont, base))
}
#[inline]
pub fn primetest_u64_miller_2bases_sub(mont: &Mont<u64>) -> bool {
assert!(mont.n < 1050535501);
MILLER_U64_2BASES.iter().all(|&base| primetest_u64_miller_sub(&mont, base))
}
#[inline]
pub fn primetest_u64_miller_3bases_sub(mont: &Mont<u64>) -> bool {
assert!(mont.n < 350269456337);
MILLER_U64_3BASES.iter().all(|&base| primetest_u64_miller_sub(&mont, base))
}
#[inline]
pub fn primetest_u64_miller_4bases_sub(mont: &Mont<u64>) -> bool {
assert!(mont.n < 55245642489451);
MILLER_U64_4BASES.iter().all(|&base| primetest_u64_miller_sub(&mont, base))
}
#[inline]
pub fn primetest_u64_miller_5bases_sub(mont: &Mont<u64>) -> bool {
assert!(mont.n < 7999252175582851);
MILLER_U64_5BASES.iter().all(|&base| primetest_u64_miller_sub(&mont, base))
}
#[inline]
pub fn primetest_u64_miller_6bases_sub(mont: &Mont<u64>) -> bool {
assert!(mont.n < 585226005592931977);
MILLER_U64_6BASES.iter().all(|&base| primetest_u64_miller_sub(&mont, base))
}
#[inline]
pub fn primetest_u64_miller_7bases_sub(mont: &Mont<u64>) -> bool {
MILLER_U64_7BASES.iter().all(|&base| primetest_u64_miller_sub(&mont, base))
}
#[inline]
pub fn primetest_u64_miller_mbases_sub(mont: &Mont<u64>) -> bool {
match mont.n {
0..=341530 => MILLER_U64_1BASES.iter(),
0..=1050535500 => MILLER_U64_2BASES.iter(),
0..=350269456336 => MILLER_U64_3BASES.iter(),
0..=55245642489450 => MILLER_U64_4BASES.iter(),
0..=7999252175582850 => MILLER_U64_5BASES.iter(),
0..=585226005592931976 => MILLER_U64_6BASES.iter(),
_ => MILLER_U64_7BASES.iter(),
}.all(|&base| primetest_u64_miller_sub(&mont, base))
}
pub fn primetest_u64_miller_mbases(n: u64) -> bool {
// Deterministic variants of the Miller-Rabin primality test
// http://miller-rabin.appspot.com/
if n == 2 { return true; }
if n == 1 || (n & 1) == 0 { return false; }
let mont = Mont::<u64>::new(n);
primetest_u64_miller_mbases_sub(&mont)
}
#[inline]
pub fn primetest_u32_miller_mbases_sub(mont: &Mont<u32>) -> bool {
match mont.n {
// Deterministic variants of the Miller-Rabin primality test
// http://miller-rabin.appspot.com/
0..=49140 => [921211727].iter(),
0..=360018360 => [1143370,2350307676].iter(),
_ => [2,7,61].iter(),
}.all(|&base| primetest_u32_miller_sub(&mont, base))
}
pub fn primetest_u32_miller_mbases(n: u32) -> bool {
// Deterministic variants of the Miller-Rabin primality test
// http://miller-rabin.appspot.com/
if n == 2 { return true; }
if n == 1 || (n & 1) == 0 { return false; }
let mont = Mont::<u32>::new(n);
primetest_u32_miller_mbases_sub(&mont)
}
fn main() {
use std::io::{BufRead,Write};
let start_time = std::time::Instant::now();
let out = std::io::stdout();
let mut out = std::io::BufWriter::new(out.lock());
macro_rules! puts {($($format:tt)*) => (let _ = write!(out,$($format)*););}
let input = std::io::stdin();
let mut lines = std::io::BufReader::new(input.lock()).lines();
let n: usize = lines.next().unwrap().unwrap().parse().unwrap();
/*
let x: Vec<u64> = lines.take(n).map(|l| l.unwrap().parse().unwrap()).collect();
let elapsed1 = start_time.elapsed();
let res: Vec<bool> = x.iter().map(|&v| prime_test_bpsw(v)).collect();
let elapsed2 = start_time.elapsed();
for i in 0..n {
puts!("{} {}\n", x[i], if res[i] { "1" } else { "0" });
}
let elapsed3 = start_time.elapsed();
out.flush().unwrap();
let elapsed4 = start_time.elapsed();
eprint!(
" input: {}us\ncompute: {}us\n output: {}us\n wflush: {}us\n",
elapsed1.as_micros(),
elapsed2.as_micros(),
elapsed3.as_micros(),
elapsed4.as_micros(),
);
*/
for _ in 0..n {
let x: u64 = lines.next().unwrap().unwrap().parse().unwrap();
puts!("{} {}\n", x, if primetest_u64_bpsw(x) { "1" } else { "0" });
}
eprint!("{}us\n", start_time.elapsed().as_micros());
}
#[cfg(test)]
mod tests {
use crate::*;
#[test]
fn test_nbits() {
// ten least k's for which (2**n)-k is prime
// https://primes.utm.edu/lists/2small/0bit.html
let primes: Vec<(u32,Vec<u64>)> = vec![
(8,vec![5,15,17,23,27,29,33,45,57,59]),
(9,vec![3,9,13,21,25,33,45,49,51,55]),
(10,vec![3,5,11,15,27,33,41,47,53,57]),
(11,vec![9,19,21,31,37,45,49,51,55,61]),
(12,vec![3,5,17,23,39,45,47,69,75,77]),
(13,vec![1,13,21,25,31,45,69,75,81,91]),
(14,vec![3,15,21,23,35,45,51,65,83,111]),
(15,vec![19,49,51,55,61,75,81,115,121,135]),
(16,vec![15,17,39,57,87,89,99,113,117,123]),
(17,vec![1,9,13,31,49,61,63,85,91,99]),
(18,vec![5,11,17,23,33,35,41,65,75,93]),
(19,vec![1,19,27,31,45,57,67,69,85,87]),
(20,vec![3,5,17,27,59,69,129,143,153,185]),
(21,vec![9,19,21,55,61,69,105,111,121,129]),
(22,vec![3,17,27,33,57,87,105,113,117,123]),
(23,vec![15,21,27,37,61,69,135,147,157,159]),
(24,vec![3,17,33,63,75,77,89,95,117,167]),
(25,vec![39,49,61,85,91,115,141,159,165,183]),
(26,vec![5,27,45,87,101,107,111,117,125,135]),
(27,vec![39,79,111,115,135,187,199,219,231,235]),
(28,vec![57,89,95,119,125,143,165,183,213,273]),
(29,vec![3,33,43,63,73,75,93,99,121,133]),
(30,vec![35,41,83,101,105,107,135,153,161,173]),
(31,vec![1,19,61,69,85,99,105,151,159,171]),
(32,vec![5,17,65,99,107,135,153,185,209,267]),
(33,vec![9,25,49,79,105,285,301,303,321,355]),
(34,vec![41,77,113,131,143,165,185,207,227,281]),
(35,vec![31,49,61,69,79,121,141,247,309,325]),
(36,vec![5,17,23,65,117,137,159,173,189,233]),
(37,vec![25,31,45,69,123,141,199,201,351,375]),
(38,vec![45,87,107,131,153,185,191,227,231,257]),
(39,vec![7,19,67,91,135,165,219,231,241,301]),
(40,vec![87,167,195,203,213,285,293,299,389,437]),
(41,vec![21,31,55,63,73,75,91,111,133,139]),
(42,vec![11,17,33,53,65,143,161,165,215,227]),
(43,vec![57,67,117,175,255,267,291,309,319,369]),
(44,vec![17,117,119,129,143,149,287,327,359,377]),
(45,vec![55,69,81,93,121,133,139,159,193,229]),
(46,vec![21,57,63,77,167,197,237,287,305,311]),
(47,vec![115,127,147,279,297,339,435,541,619,649]),
(48,vec![59,65,89,93,147,165,189,233,243,257]),
(49,vec![81,111,123,139,181,201,213,265,283,339]),
(50,vec![27,35,51,71,113,117,131,161,195,233]),
(51,vec![129,139,165,231,237,247,355,391,397,439]),
(52,vec![47,143,173,183,197,209,269,285,335,395]),
(53,vec![111,145,231,265,315,339,343,369,379,421]),
(54,vec![33,53,131,165,195,245,255,257,315,327]),
(55,vec![55,67,99,127,147,169,171,199,207,267]),
(56,vec![5,27,47,57,89,93,147,177,189,195]),
(57,vec![13,25,49,61,69,111,195,273,363,423]),
(58,vec![27,57,63,137,141,147,161,203,213,251]),
(59,vec![55,99,225,427,517,607,649,687,861,871]),
(60,vec![93,107,173,179,257,279,369,395,399,453]),
(61,vec![1,31,45,229,259,283,339,391,403,465]),
(62,vec![57,87,117,143,153,167,171,195,203,273]),
(63,vec![25,165,259,301,375,387,391,409,457,471]),
(64,vec![59,83,95,179,189,257,279,323,353,363]),
];
for (bit_ref, kvec) in primes.iter() {
let bit = *bit_ref;
let &lastk = kvec.iter().last().unwrap();
for k in 1..=lastk {
if (k & 1) == 0 { continue; }
let n = (if bit < 64 { 1u64 << bit } else { 0u64 }).wrapping_sub(k);
let f = kvec.binary_search(&k).is_ok();
let mont = Mont::<u64>::new(n);
// Deterministic variants of the Miller-Rabin primality test
// http://miller-rabin.appspot.com/
// n < 341531 : 1 base Miller-Rabin primality test
if n < 341531 { assert_eq!(primetest_u64_miller_1base_sub(&mont), f); }
// n < 1050535501 : 2 bases Miller-Rabin primality test
if n < 1050535501 { assert_eq!(primetest_u64_miller_2bases_sub(&mont), f); }
// n < 350269456337 : 3 bases Miller-Rabin primality test
if n < 350269456337 { assert_eq!(primetest_u64_miller_3bases_sub(&mont), f); }
// n < 55245642489451 : 4 bases Miller-Rabin primality test
if n < 55245642489451 { assert_eq!(primetest_u64_miller_4bases_sub(&mont), f); }
// n < 7999252175582851 : 5 bases Miller-Rabin primality test
if n < 7999252175582851 { assert_eq!(primetest_u64_miller_5bases_sub(&mont), f); }
// n < 585226005592931977 : 6 bases Miller-Rabin primality test
if n < 585226005592931977 { assert_eq!(primetest_u64_miller_6bases_sub(&mont), f); }
// n for all 64bit integer : 7 bases Miller-Rabin primality test
assert_eq!(primetest_u64_miller_7bases_sub(&mont), f);
// n for all 64bit integer : multi bases Miller-Rabin primality test
assert_eq!(primetest_u64_miller_mbases_sub(&mont), f);
// Baillie–PSW primality test
assert_eq!(primetest_u64_bpsw_sub(&mont), f);
// Lucas test
assert_eq!(primetest_u64_lucas_sub(&mont), primetest_u64_lucas_sub(&mont));
// 32bit multi bases test
if n < 4294967296 {
let mont32 = Mont::<u32>::new(n as u32);
assert_eq!(primetest_u32_miller_mbases_sub(&mont32), f);
}
}
}
}
#[test]
fn test_base2_1e7() {
// Miller-Rabin primality test (base 2)
// strong pseudoprimes to base 2 ( https://oeis.org/A001262 ): 2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,...
// composite 2-SPRP list up to 2**64 ( http://miller-rabin.appspot.com/#links )
// Pseudoprime Statistics, Tables, and Data ( http://ntheory.org/pseudoprimes.html )
// Miller-Rabin base 2 data (up to 1e15) ( http://ntheory.org/data/spsps.txt )
let assumed: Vec<u64> = vec![ // #SPSP-2 Miller-Rabin base 2 (up to 1e7)
2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,65281,74665,80581,85489,88357,90751,104653,
130561,196093,220729,233017,252601,253241,256999,271951,280601,314821,357761,390937,458989,476971,
486737,489997,514447,580337,635401,647089,741751,800605,818201,838861,873181,877099,916327,976873,
983401,1004653,1016801,1023121,1082401,1145257,1194649,1207361,1251949,1252697,1302451,1325843,1357441,
1373653,1397419,1441091,1493857,1507963,1509709,1530787,1678541,1730977,1811573,1876393,1907851,1909001,
1969417,1987021,2004403,2081713,2181961,2205967,2264369,2269093,2284453,2304167,2387797,2419385,2510569,
2746477,2748023,2757241,2811271,2909197,2953711,2976487,3090091,3116107,3125281,3375041,3400013,3429037,
3539101,3567481,3581761,3605429,3898129,4181921,4188889,4335241,4360621,4469471,4502485,4513841,4682833,
4835209,4863127,5016191,5044033,5049001,5173169,5173601,5256091,5310721,5444489,5489641,5590621,5599765,
5672041,5681809,5919187,6140161,6226193,6233977,6334351,6368689,6386993,6787327,6836233,6952037,7177105,
7306261,7306561,7462001,7674967,7759937,7820201,7883731,8036033,8095447,8384513,8388607,8534233,8725753,
8727391,9006401,9056501,9069229,9073513,9371251,9564169,9567673,9588151,9729301,9774181,9863461,9995671
];
let result: Vec<u64> = (3..10_000_000).filter(|&n| {
if (n & 1) == 0 { return false; }
let mont = Mont::<u64>::new(n);
let res_3bases = primetest_u64_miller_3bases_sub(&mont);
let res_base2 = primetest_u64_base2_sub(&mont);
assert!(!res_3bases || res_base2);
res_3bases != res_base2
}).collect();
assert_eq!(assumed, result);
}
#[test]
fn test_base2_1e7u32() {
// Miller-Rabin primality test (base 2)
// strong pseudoprimes to base 2 ( https://oeis.org/A001262 ): 2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,...
// composite 2-SPRP list up to 2**64 ( http://miller-rabin.appspot.com/#links )
// Pseudoprime Statistics, Tables, and Data ( http://ntheory.org/pseudoprimes.html )
// Miller-Rabin base 2 data (up to 1e15) ( http://ntheory.org/data/spsps.txt )
let assumed: Vec<u32> = vec![ // #SPSP-2 Miller-Rabin base 2 (up to 1e7)
2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,65281,74665,80581,85489,88357,90751,104653,
130561,196093,220729,233017,252601,253241,256999,271951,280601,314821,357761,390937,458989,476971,
486737,489997,514447,580337,635401,647089,741751,800605,818201,838861,873181,877099,916327,976873,
983401,1004653,1016801,1023121,1082401,1145257,1194649,1207361,1251949,1252697,1302451,1325843,1357441,
1373653,1397419,1441091,1493857,1507963,1509709,1530787,1678541,1730977,1811573,1876393,1907851,1909001,
1969417,1987021,2004403,2081713,2181961,2205967,2264369,2269093,2284453,2304167,2387797,2419385,2510569,
2746477,2748023,2757241,2811271,2909197,2953711,2976487,3090091,3116107,3125281,3375041,3400013,3429037,
3539101,3567481,3581761,3605429,3898129,4181921,4188889,4335241,4360621,4469471,4502485,4513841,4682833,
4835209,4863127,5016191,5044033,5049001,5173169,5173601,5256091,5310721,5444489,5489641,5590621,5599765,
5672041,5681809,5919187,6140161,6226193,6233977,6334351,6368689,6386993,6787327,6836233,6952037,7177105,
7306261,7306561,7462001,7674967,7759937,7820201,7883731,8036033,8095447,8384513,8388607,8534233,8725753,
8727391,9006401,9056501,9069229,9073513,9371251,9564169,9567673,9588151,9729301,9774181,9863461,9995671
];
let result: Vec<u32> = (3..10_000_000).filter(|&n| {
if (n & 1) == 0 { return false; }
let mont = Mont::<u32>::new(n);
let res_mbases = primetest_u32_miller_mbases_sub(&mont);
let res_base2 = primetest_u32_base2_sub(&mont);
assert!(!res_mbases || res_base2);
res_mbases != res_base2
}).collect();
assert_eq!(assumed, result);
}
#[test]
fn test_lucas_1e7() {
// Lucas primality test
// strong Lucas pseudoprimes ( https://oeis.org/A217255 ): 5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,...
// Pseudoprime Statistics, Tables, and Data ( http://ntheory.org/pseudoprimes.html )
// Strong Lucas-Selfridge data (up to 1e15) ( http://ntheory.org/data/slpsps-baillie.txt )
let assumed: Vec<u64> = vec![ // #SLPSP Strong Lucas-Selfridge (up to 1e7)
5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,75077,97439,100127,113573,115639,130139,155819,
158399,161027,162133,176399,176471,189419,192509,197801,224369,230691,231703,243629,253259,268349,288919,
313499,324899,353219,366799,391169,430127,436409,455519,487199,510479,572669,611399,622169,635627,636199,
701999,794611,835999,839159,851927,871859,875879,887879,895439,950821,960859,1033997,1106327,1241099,
1256293,1308119,1311389,1388903,1422319,1501439,1697183,1711469,1777159,1981559,2003579,2263127,2435423,
2461211,2518889,2566409,2624399,2662277,2666711,2690759,2738969,2782079,2828699,2942081,2952071,3109049,
3165119,3175883,3179609,3204599,3373649,3399527,3410531,3441239,3452147,3479111,3498879,3579599,3684251,
3694079,3700559,3706169,3735521,3774377,3776219,3785699,3802499,3813011,3865319,3892529,3900797,3903791,
4067279,4109363,4226777,4309631,4322399,4368869,4403027,4563719,4828277,4870847,5133281,5208377,5299139,
5396999,5450201,5479109,5514479,5720219,5762629,5807759,5879411,5942627,6001379,6003923,6296291,6562891,
6641189,6668099,6784721,6784861,6863291,6893531,6965639,7017949,7163441,7199399,7241639,7353917,7453619,
7621499,8112899,8159759,8221121,8234159,8361989,8372849,8518127,8530559,8555009,8574551,8581219,8711699,
8817899,8990279,9049319,9335969,9401893,9485951,9587411,9713027,9793313,9800981,9827711,9922337,9965069
];
let result: Vec<u64> = (3..10_000_000).filter(|&n| {
if (n & 1) == 0 { return false; }
let mont = Mont::<u64>::new(n);
let res_3bases = primetest_u64_miller_3bases_sub(&mont);
let res_lucas = primetest_u64_lucas_sub(&mont);
assert!(!res_3bases || res_lucas);
res_3bases != res_lucas
}).collect();
assert_eq!(assumed, result);
}
#[test]
fn test_bpsw_1e7() { // 24bit
for n in 3..10_000_000 {
if (n & 1) == 0 { continue; }
let mont = Mont::<u64>::new(n);
let mont32 = Mont::<u32>::new(n as u32);
let res_3bases = primetest_u64_miller_3bases_sub(&mont);
let res_bpsw = primetest_u64_bpsw_sub(&mont);
let res_u32mbases = primetest_u32_miller_mbases_sub(&mont32);
assert_eq!(res_3bases, res_bpsw);
assert_eq!(res_3bases, res_u32mbases);
}
}
#[test]
fn test_bpsw_4e9() { // 32bit
for n in 4_000_000_000..4_010_000_000 {
if (n & 1) == 0 { continue; }
let mont = Mont::<u64>::new(n);
let mont32 = Mont::<u32>::new(n as u32);
let res_3bases = primetest_u64_miller_3bases_sub(&mont);
let res_bpsw = primetest_u64_bpsw_sub(&mont);
let res_u32mbases = primetest_u32_miller_mbases_sub(&mont32);
assert_eq!(res_3bases, res_bpsw);
assert_eq!(res_3bases, res_u32mbases);
}
}
#[test]
fn test_bpsw_1e16() { // 54bit
for n in 10_000_000_000_000_000..10_000_000_010_000_000 {
if (n & 1) == 0 { continue; }
let mont = Mont::<u64>::new(n);
let res_7bases = primetest_u64_miller_7bases_sub(&mont);
let res_bpsw = primetest_u64_bpsw_sub(&mont);
assert_eq!(res_7bases, res_bpsw);
}
}
#[test]
fn test_bpsw_9e18() { // 63bit
for n in 9_000_000_000_000_000_000..9_000_000_000_010_000_000 {
if (n & 1) == 0 { continue; }
let mont = Mont::<u64>::new(n);
let res_7bases = primetest_u64_miller_7bases_sub(&mont);
let res_bpsw = primetest_u64_bpsw_sub(&mont);
assert_eq!(res_7bases, res_bpsw);
}
}
#[test]
fn test_bpsw_10e18() { // 64bit
for n in 10_000_000_000_000_000_000..10_000_000_000_010_000_000 {
if (n & 1) == 0 { continue; }
let mont = Mont::<u64>::new(n);
let res_7bases = primetest_u64_miller_7bases_sub(&mont);
let res_bpsw = primetest_u64_bpsw_sub(&mont);
assert_eq!(res_7bases, res_bpsw);
}
}
#[test]
fn test_bpsw_18e18() { // 64bit
for n in 18_000_000_000_000_000_000..18_000_000_000_010_000_000 {
if (n & 1) == 0 { continue; }
let mont = Mont::<u64>::new(n);
let res_7bases = primetest_u64_miller_7bases_sub(&mont);
let res_bpsw = primetest_u64_bpsw_sub(&mont);
assert_eq!(res_7bases, res_bpsw);
}
}
}