結果
| 問題 | No.3030 ミラー・ラビン素数判定法のテスト |
| ユーザー |
👑  Mizar
|
| 提出日時 | 2022-09-05 12:28:55 |
| 言語 | Rust (1.77.0 + proconio) |
| 結果 |
AC
|
| 実行時間 | 14 ms / 9,973 ms |
| コード長 | 42,978 bytes |
| コンパイル時間 | 20,516 ms |
| コンパイル使用メモリ | 406,808 KB |
| 実行使用メモリ | 6,820 KB |
| 最終ジャッジ日時 | 2024-11-20 11:53:05 |
| 合計ジャッジ時間 | 14,962 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge5 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 1 ms
6,816 KB |
| testcase_01 | AC | 1 ms
6,820 KB |
| testcase_02 | AC | 1 ms
6,820 KB |
| testcase_03 | AC | 1 ms
6,816 KB |
| testcase_04 | AC | 10 ms
6,816 KB |
| testcase_05 | AC | 10 ms
6,820 KB |
| testcase_06 | AC | 9 ms
6,816 KB |
| testcase_07 | AC | 8 ms
6,820 KB |
| testcase_08 | AC | 9 ms
6,816 KB |
| testcase_09 | AC | 14 ms
6,816 KB |
ソースコード
// -*- coding:utf-8-unix -*-
use std::num::{NonZeroU128,NonZeroU64};
// 除数をNonZeroU型にした除算/剰余算の選択
macro_rules! cond_nzdiv {
($nzdiv:expr, $no_nzdiv:expr) => {
$nzdiv // 1.51.0以降 feature = "nonzero_div"
//$no_nzdiv // 1.50.0以前
}
}
// NonZeroU64/NonderoU128型を使用するか
const USE_NZU64: bool = true;
// モンゴメリ剰余乗算 (Montgomery modular multiplication)
trait MontVal<T, BitCountType=u32> {
fn n(&self) -> T;
fn ni(&self) -> T;
fn nh(&self) -> T;
fn r(&self) -> T;
fn rn(&self) -> T;
fn r2(&self) -> T;
fn d(&self) -> T;
fn k(&self) -> BitCountType;
}
trait MontOps<T, BitCountType=u32> {
fn add(&self, a: T, b: T) -> T;
fn sub(&self, a: T, b: T) -> T;
fn div2(&self, ar: T) -> T;
fn mrmul(&self, ar: T, br: T) -> T;
fn mr(&self, ar: T) -> T;
fn ar(&self, a: T) -> T;
fn pow(&self, ar: T, b: T) -> T;
fn powodd(&self, ar: T, b: T) -> T;
}
trait MontModN<T, BitCountType=u32> {
fn modn(&self, x: T) -> T;
}
trait MontNew<T, U=T, BitCountType=u32> {
fn new(n: T) -> Self;
}
struct Mont<T, BitCountType=u32> {
n: T, // n is odd, and n > 2
ni: T, // n * ni == 1 (mod 2**64)
nh: T, // == (n + 1) / 2
r: T, // == 2**64 (mod n)
rn: T, // == -(2**64) (mod n)
r2: T, // == 2**128 (mod n)
d: T, // == (n - 1) >> (n - 1).trailing_zeros()
k: BitCountType, // == (n - 1).trailing_zeros()
}
impl MontVal<u64> for Mont<NonZeroU64> {
#[inline] fn n(&self) -> u64 { self.n.get() }
#[inline] fn ni(&self) -> u64 { self.ni.get() }
#[inline] fn nh(&self) -> u64 { self.nh.get() }
#[inline] fn r(&self) -> u64 { self.r.get() }
#[inline] fn rn(&self) -> u64 { self.rn.get() }
#[inline] fn r2(&self) -> u64 { self.r2.get() }
#[inline] fn d(&self) -> u64 { self.d.get() }
#[inline] fn k(&self) -> u32 { self.k }
}
impl MontVal<u64> for Mont<u64> {
#[inline] fn n(&self) -> u64 { self.n }
#[inline] fn ni(&self) -> u64 { self.ni }
#[inline] fn nh(&self) -> u64 { self.nh }
#[inline] fn r(&self) -> u64 { self.r }
#[inline] fn rn(&self) -> u64 { self.rn }
#[inline] fn r2(&self) -> u64 { self.r2 }
#[inline] fn d(&self) -> u64 { self.d }
#[inline] fn k(&self) -> u32 { self.k }
}
macro_rules! mont_ops_u64 { ($x:ty) => { impl MontOps<u64> for $x {
#[inline]
fn add(&self, a: u64, b: u64) -> u64 {
// == 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 {
// == 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 {
// == 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 {
// == (ar * br) / r (mod n)
debug_assert!(ar < self.n());
debug_assert!(br < self.n());
let (n, ni) = (self.n(), self.ni());
let t: u128 = (ar as u128) * (br as u128);
let (t, f) = ((t >> 64) as u64).overflowing_sub((((((t as u64).wrapping_mul(ni)) as u128) * (n as u128)) >> 64) as u64);
if f { t.wrapping_add(n) } else { t }
}
#[inline]
fn mr(&self, ar: u64) -> u64 {
// == ar / r (mod n)
debug_assert!(ar < self.n());
let (n, ni) = (self.n(), self.ni());
let (t, f) = (((((ar.wrapping_mul(ni)) as u128) * (n as u128)) >> 64) as u64).overflowing_neg();
if f { t.wrapping_add(n) } else { t }
}
#[inline]
fn ar(&self, a: u64) -> u64 {
// == 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); // b is odd
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); }
}
}
}}}
mont_ops_u64!(Mont<NonZeroU64>);
mont_ops_u64!(Mont<u64>);
impl MontModN<u64> for Mont<NonZeroU64> {
#[inline] fn modn(&self, x: u64) -> u64 { if x < self.n() { x } else { cond_nzdiv!(x % self.n, x % self.n()) } }
}
impl MontModN<u64> for Mont<u64> {
#[inline] fn modn(&self, x: u64) -> u64 { if x < self.n() { x } else { x % self.n() } }
}
impl MontNew<NonZeroU64, u64> for Mont<NonZeroU64> {
#[inline]
fn new(n: NonZeroU64) -> Self {
debug_assert_eq!(n.get() & 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.get();
// 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 _ in 0..5 {
ni = ni.wrapping_mul(2u64.wrapping_sub(n.get().wrapping_mul(ni)));
}
debug_assert_eq!(n.get().wrapping_mul(ni), 1); // n * ni == 1 (mod 2**64)
let nh = (n.get() >> 1) + 1;
let r = cond_nzdiv!(
1u64.wrapping_neg() % n,
1u64.wrapping_neg() % n.get()
) + 1; // == 2**64 (mod n)
let rn = n.get() - r;
let r2 = cond_nzdiv!(
1u128.wrapping_neg() % NonZeroU128::from(n),
1u128.wrapping_neg() % NonZeroU128::from(n).get()
) as u64 + 1; // == 2**128 (mod n)
let mut d = n.get() - 1;
let k = d.trailing_zeros();
d >>= k;
let ni = NonZeroU64::new(ni).unwrap();
let nh = NonZeroU64::new(nh).unwrap();
let r = NonZeroU64::new(r).unwrap();
let rn = NonZeroU64::new(rn).unwrap();
let r2 = NonZeroU64::new(r2).unwrap();
let d = NonZeroU64::new(d).unwrap();
debug_assert_eq!(Self { n, ni, nh, r, rn, r2, d, k }.mr(r.get()), 1); // r / r == 1 (mod n)
debug_assert_eq!(Self { n, ni, nh, r, rn, r2, d, k }.mrmul(1, r2.get()), r.get()); // r2 / r == r (mod n)
Self { n, ni, nh, r, rn, r2, d, k }
}
}
impl MontNew<u64> for Mont<u64> {
#[inline]
fn new(n: u64) -> Self {
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 * (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 _ in 0..5 {
ni = ni.wrapping_mul(2u64.wrapping_sub(n.wrapping_mul(ni)));
}
debug_assert_eq!(n.wrapping_mul(ni), 1); // n * ni == 1 (mod 2**64)
let nh = (n >> 1) + 1;
let r = (1u64.wrapping_neg() % n) + 1; // == 2**64 (mod n)
let rn = n - r;
let r2 = ((1u128.wrapping_neg() % (n as u128)) as u64) + 1; // == 2**128 (mod n)
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 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
}
#[inline]
// 平方数判定 (true: 平方数である, false: 平方数ではない)
fn issq_u64_n(x: u64) -> bool {
let sqrt = isqrt64n(x);
sqrt * sqrt == x
}
#[inline]
// 擬平方数判定 mod64 12/64 0.1875 (true: 平方数かもしれない, false: 平方数ではない)
fn issq_u64_mod64(x: u64) -> bool {
(0x202021202030213u64 >> (x & 63)) & 1 == 1
}
#[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
}
// 平方数判定 (true: 平方数である, false: 平方数ではない)
#[inline]
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 {
if n == 0 { return if a == 1 || a == -1 { 1 } else { 0 }; }
let (mut a, mut j): (u64, i32) = if a >= 0 { (a as u64, 1) } else if (n & 3) == 3 { ((-a) as u64, -1) } else { ((-a) as u64, 1) };
while a != 0 {
let ba = a.trailing_zeros();
if (ba & 1) != 0 && ((n & 7) == 3 || (n & 7) == 5) { j = -j; }
a >>= a.trailing_zeros();
if (a & n & 3) == 3 { j = -j; }
std::mem::swap(&mut n, &mut a);
a %= n;
if a > (n >> 1) {
a = n - a;
if (n & 3) == 3 { j = -j; }
}
}
if n == 1 { j } else { 0 }
}
trait PrimeTestU64Trait<T> : MontOps<u64> + MontVal<u64> + MontModN<u64> {
#[inline]
fn primetest_base2(&self) -> bool {
// Mirrer-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) = (self.r(), self.rn(), self.d(), self.k());
let mut br = self.powodd(self.add(self.r(), self.r()), d);
if br == r || br == rn { return true; }
for _ in 1..k {
br = self.mrmul(br, br);
if br == rn { return true; }
}
false
}
#[inline]
fn primetest_miller(&self, mut base: u64) -> bool {
// Miller-Rabin primality test
let (n, r, rn, d, k) = (self.n(), self.r(), self.rn(), self.d(), self.k());
if base > n { base = self.modn(base); if base == 0 { return true; } }
let mut tr = self.pow(self.ar(base), d);
if tr == r || tr == rn { return true; }
for _ in 1..k { tr = self.mrmul(tr, tr); if tr == rn { return true; } }
false
}
#[inline]
fn primetest_lucas(&self) -> bool {
// Lucas primality test
// strong Lucas pseudoprimes ( https://oeis.org/A217255 ): 5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,...
let n = self.n();
let mut d: i64 = 5;
for i in 0u32.. {
debug_assert!(i < 64);
match jacobi(d, n) {
-1 => break,
0 => if ((d.abs()) as u64) < n { return false; },
_ => {},
}
if i == 8 && issq_u64(n) { return false; }
if (i & 1) == 1 { d = 2 - d; } else { d = -(d + 2); }
}
let qm = self.ar(self.modn(if d < 0 {((1 - d) as u64) / 4} else {n - ((d - 1) as u64) / 4}));
let mut k = (n + 1) << (n + 1).leading_zeros();
let mut um = self.r();
let mut vm = self.r();
let mut qn = qm;
let dm: u64 = self.ar(if d < 0 { n - self.modn((-d) as u64) } else { self.modn(d as u64) });
k <<= 1;
while k != 0 {
um = self.mrmul(um, vm);
vm = self.sub(self.mrmul(vm, vm), self.add(qn, qn));
qn = self.mrmul(qn, qn);
if (k >> 63) != 0 {
let uu = self.div2(self.add(um, vm));
vm = self.div2(self.add(self.mrmul(dm, um), vm));
um = uu;
qn = self.mrmul(qn, qm);
}
k <<= 1;
}
if um == 0 || vm == 0 {
return true;
}
let mut x = (n + 1) & (!n);
x >>= 1;
while x != 0 {
um = self.mrmul(um, vm);
vm = self.sub(self.mrmul(vm, vm), self.add(qn, qn));
if vm == 0 {
return true;
}
qn = self.mrmul(qn, qn);
x >>= 1;
}
false
}
#[inline]
fn primetest_bpsw(&self) -> bool {
// Mirrer-Rabin primality test (base 2)
// strong pseudoprimes to base 2 ( https://oeis.org/A001262 ): 2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,...
self.primetest_base2() &&
// Lucas primality test
// strong Lucas pseudoprimes ( https://oeis.org/A217255 ): 5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,...
self.primetest_lucas()
}
}
impl PrimeTestU64Trait<NonZeroU64> for Mont<NonZeroU64> {}
impl PrimeTestU64Trait<u64> for Mont<u64> {}
// Deterministic variants of the Miller-Rabin primality test http://miller-rabin.appspot.com/
trait MillerU64<T> : PrimeTestU64Trait<T> {
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]
fn primetest_miller_1bases(&self) -> bool {
debug_assert!(self.n() < 341531);
Self::MILLER_U64_1BASES.iter().all(|&base| self.primetest_miller(base))
}
#[inline]
fn primetest_miller_2bases(&self) -> bool {
debug_assert!(self.n() < 1050535501);
Self::MILLER_U64_2BASES.iter().all(|&base| self.primetest_miller(base))
}
#[inline]
fn primetest_miller_3bases(&self) -> bool {
debug_assert!(self.n() < 350269456337);
Self::MILLER_U64_3BASES.iter().all(|&base| self.primetest_miller(base))
}
#[inline]
fn primetest_miller_4bases(&self) -> bool {
debug_assert!(self.n() < 55245642489451);
Self::MILLER_U64_4BASES.iter().all(|&base| self.primetest_miller(base))
}
#[inline]
fn primetest_miller_5bases(&self) -> bool {
debug_assert!(self.n() < 7999252175582851);
Self::MILLER_U64_5BASES.iter().all(|&base| self.primetest_miller(base))
}
#[inline]
fn primetest_miller_6bases(&self) -> bool {
debug_assert!(self.n() < 585226005592931977);
Self::MILLER_U64_6BASES.iter().all(|&base| self.primetest_miller(base))
}
#[inline]
fn primetest_miller_7bases(&self) -> bool {
Self::MILLER_U64_7BASES.iter().all(|&base| self.primetest_miller(base))
}
#[inline]
fn primetest_miller_mbases(&self) -> bool {
// Deterministic variants of the Miller-Rabin primality test http://miller-rabin.appspot.com/
match self.n() {
0..=341530 => Self::MILLER_U64_1BASES.iter(),
0..=1050535500 => Self::MILLER_U64_2BASES.iter(),
0..=350269456336 => Self::MILLER_U64_3BASES.iter(),
0..=55245642489450 => Self::MILLER_U64_4BASES.iter(),
0..=7999252175582850 => Self::MILLER_U64_5BASES.iter(),
0..=585226005592931976 => Self::MILLER_U64_6BASES.iter(),
_ => Self::MILLER_U64_7BASES.iter(),
}.all(|&base| self.primetest_miller(base))
}
}
impl MillerU64<NonZeroU64> for Mont<NonZeroU64> {}
impl MillerU64<u64> for Mont<u64> {}
// Baillie–PSW primarity test (nonzerou64)
pub fn primetest_nzu64_bpsw(n: NonZeroU64) -> bool {
if n.get() == 2 { return true; }
if n.get() == 1 || (n.get() & 1) == 0 { return false; }
let mont = Mont::<NonZeroU64>::new(n);
mont.primetest_bpsw()
}
// Baillie–PSW primarity test (u64)
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);
mont.primetest_bpsw()
}
// Miller-Rabin primarity test (nonzerou64)
pub fn primetest_nzu64_miller_mbases(n: NonZeroU64) -> bool {
if n.get() == 2 { return true; }
if n.get() == 1 || (n.get() & 1) == 0 { return false; }
let mont = Mont::<NonZeroU64>::new(n);
mont.primetest_miller_mbases()
}
// Miller-Rabin primarity test (u64)
pub fn primetest_u64_miller_mbases(n: u64) -> bool {
if n == 2 { return true; }
if n == 1 || (n & 1) == 0 { return false; }
let mont = Mont::<u64>::new(n);
mont.primetest_miller_mbases()
}
#[cfg(test)]
mod tests {
use crate::*;
// ten least k's for which (2**n)-k is prime
// https://primes.utm.edu/lists/2small/0bit.html
const PRIMES_LAST10: [(u32,[u64;10]);57] = [
(8,[5,15,17,23,27,29,33,45,57,59]),
(9,[3,9,13,21,25,33,45,49,51,55]),
(10,[3,5,11,15,27,33,41,47,53,57]),
(11,[9,19,21,31,37,45,49,51,55,61]),
(12,[3,5,17,23,39,45,47,69,75,77]),
(13,[1,13,21,25,31,45,69,75,81,91]),
(14,[3,15,21,23,35,45,51,65,83,111]),
(15,[19,49,51,55,61,75,81,115,121,135]),
(16,[15,17,39,57,87,89,99,113,117,123]),
(17,[1,9,13,31,49,61,63,85,91,99]),
(18,[5,11,17,23,33,35,41,65,75,93]),
(19,[1,19,27,31,45,57,67,69,85,87]),
(20,[3,5,17,27,59,69,129,143,153,185]),
(21,[9,19,21,55,61,69,105,111,121,129]),
(22,[3,17,27,33,57,87,105,113,117,123]),
(23,[15,21,27,37,61,69,135,147,157,159]),
(24,[3,17,33,63,75,77,89,95,117,167]),
(25,[39,49,61,85,91,115,141,159,165,183]),
(26,[5,27,45,87,101,107,111,117,125,135]),
(27,[39,79,111,115,135,187,199,219,231,235]),
(28,[57,89,95,119,125,143,165,183,213,273]),
(29,[3,33,43,63,73,75,93,99,121,133]),
(30,[35,41,83,101,105,107,135,153,161,173]),
(31,[1,19,61,69,85,99,105,151,159,171]),
(32,[5,17,65,99,107,135,153,185,209,267]),
(33,[9,25,49,79,105,285,301,303,321,355]),
(34,[41,77,113,131,143,165,185,207,227,281]),
(35,[31,49,61,69,79,121,141,247,309,325]),
(36,[5,17,23,65,117,137,159,173,189,233]),
(37,[25,31,45,69,123,141,199,201,351,375]),
(38,[45,87,107,131,153,185,191,227,231,257]),
(39,[7,19,67,91,135,165,219,231,241,301]),
(40,[87,167,195,203,213,285,293,299,389,437]),
(41,[21,31,55,63,73,75,91,111,133,139]),
(42,[11,17,33,53,65,143,161,165,215,227]),
(43,[57,67,117,175,255,267,291,309,319,369]),
(44,[17,117,119,129,143,149,287,327,359,377]),
(45,[55,69,81,93,121,133,139,159,193,229]),
(46,[21,57,63,77,167,197,237,287,305,311]),
(47,[115,127,147,279,297,339,435,541,619,649]),
(48,[59,65,89,93,147,165,189,233,243,257]),
(49,[81,111,123,139,181,201,213,265,283,339]),
(50,[27,35,51,71,113,117,131,161,195,233]),
(51,[129,139,165,231,237,247,355,391,397,439]),
(52,[47,143,173,183,197,209,269,285,335,395]),
(53,[111,145,231,265,315,339,343,369,379,421]),
(54,[33,53,131,165,195,245,255,257,315,327]),
(55,[55,67,99,127,147,169,171,199,207,267]),
(56,[5,27,47,57,89,93,147,177,189,195]),
(57,[13,25,49,61,69,111,195,273,363,423]),
(58,[27,57,63,137,141,147,161,203,213,251]),
(59,[55,99,225,427,517,607,649,687,861,871]),
(60,[93,107,173,179,257,279,369,395,399,453]),
(61,[1,31,45,229,259,283,339,391,403,465]),
(62,[57,87,117,143,153,167,171,195,203,273]),
(63,[25,165,259,301,375,387,391,409,457,471]),
(64,[59,83,95,179,189,257,279,323,353,363]),
];
#[test]
fn test_nzu64_nbits() {
// ten least k's for which (2**n)-k is prime
// https://primes.utm.edu/lists/2small/0bit.html
let start_time = std::time::Instant::now();
for (bit_ref, kvec) in PRIMES_LAST10.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 n = NonZeroU64::new(n).unwrap();
let f = kvec.binary_search(&k).is_ok();
let mont = Mont::<NonZeroU64>::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.get() < 341531 { assert_eq!(mont.primetest_miller_1bases(), f); }
// n < 1050535501 : 2 bases Miller-Rabin primality test
if n.get() < 1050535501 { assert_eq!(mont.primetest_miller_2bases(), f); }
// n < 350269456337 : 3 bases Miller-Rabin primality test
if n.get() < 350269456337 { assert_eq!(mont.primetest_miller_3bases(), f); }
// n < 55245642489451 : 4 bases Miller-Rabin primality test
if n.get() < 55245642489451 { assert_eq!(mont.primetest_miller_4bases(), f); }
// n < 7999252175582851 : 5 bases Miller-Rabin primality test
if n.get() < 7999252175582851 { assert_eq!(mont.primetest_miller_5bases(), f); }
// n < 585226005592931977 : 6 bases Miller-Rabin primality test
if n.get() < 585226005592931977 { assert_eq!(mont.primetest_miller_6bases(), f); }
// n for all 64bit integer : 7 bases Miller-Rabin primality test
assert_eq!(mont.primetest_miller_7bases(), f);
// n for all 64bit integer : multi bases Miller-Rabin primality test
assert_eq!(mont.primetest_miller_mbases(), f);
// Baillie–PSW primality test
assert_eq!(mont.primetest_bpsw(), f);
// Lucas test
assert_eq!(mont.primetest_lucas(), mont.primetest_lucas());
}
}
eprint!("test_nzu64_nbits {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_u64_nbits() {
// ten least k's for which (2**n)-k is prime
// https://primes.utm.edu/lists/2small/0bit.html
let start_time = std::time::Instant::now();
for (bit_ref, kvec) in PRIMES_LAST10.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!(mont.primetest_miller_1bases(), f); }
// n < 1050535501 : 2 bases Miller-Rabin primality test
if n < 1050535501 { assert_eq!(mont.primetest_miller_2bases(), f); }
// n < 350269456337 : 3 bases Miller-Rabin primality test
if n < 350269456337 { assert_eq!(mont.primetest_miller_3bases(), f); }
// n < 55245642489451 : 4 bases Miller-Rabin primality test
if n < 55245642489451 { assert_eq!(mont.primetest_miller_4bases(), f); }
// n < 7999252175582851 : 5 bases Miller-Rabin primality test
if n < 7999252175582851 { assert_eq!(mont.primetest_miller_5bases(), f); }
// n < 585226005592931977 : 6 bases Miller-Rabin primality test
if n < 585226005592931977 { assert_eq!(mont.primetest_miller_6bases(), f); }
// n for all 64bit integer : 7 bases Miller-Rabin primality test
assert_eq!(mont.primetest_miller_7bases(), f);
// n for all 64bit integer : multi bases Miller-Rabin primality test
assert_eq!(mont.primetest_miller_mbases(), f);
// Baillie–PSW primality test
assert_eq!(mont.primetest_bpsw(), f);
// Lucas test
assert_eq!(mont.primetest_lucas(), mont.primetest_lucas());
}
}
eprint!("test_u64_nbits {}us\n", start_time.elapsed().as_micros());
}
// #SPSP-2 Miller-Rabin base 2 (up to 1e7)
const SPRP_2: [u64; 162] = [
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
];
#[test]
fn test_base2_nzu64_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::from(SPRP_2);
let start_time = std::time::Instant::now();
let result: Vec<u64> = (3..10_000_000).filter(|&n| {
if (n & 1) == 0 { return false; }
let nn = NonZeroU64::new(n).unwrap();
let mont = Mont::<NonZeroU64>::new(nn);
let res_3bases = mont.primetest_miller_3bases();
let res_base2 = mont.primetest_base2();
assert!(!res_3bases || res_base2);
res_3bases != res_base2
}).collect();
assert_eq!(assumed, result);
eprint!("test_base2_nzu64_1e7 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_base2_u64_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::from(SPRP_2);
let start_time = std::time::Instant::now();
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 = mont.primetest_miller_3bases();
let res_base2 = mont.primetest_base2();
res_3bases != res_base2
}).collect();
assert_eq!(assumed, result);
eprint!("test_base2_u64_1e7 {}us\n", start_time.elapsed().as_micros());
}
// #SLPSP Strong Lucas-Selfridge (up to 1e7)
const SLPRP: [u64; 178] = [
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
];
#[test]
fn test_lucas_nzu64_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::from(SLPRP);
let start_time = std::time::Instant::now();
let result: Vec<u64> = (3..10_000_000).filter(|&n| {
if (n & 1) == 0 { return false; }
let nn = NonZeroU64::new(n).unwrap();
let mont = Mont::<NonZeroU64>::new(nn);
let res_3bases = mont.primetest_miller_3bases();
let res_lucas = mont.primetest_lucas();
assert!(!res_3bases || res_lucas);
res_3bases != res_lucas
}).collect();
assert_eq!(assumed, result);
eprint!("test_lucas_nzu64_1e7 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_lucas_u64_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::from(SLPRP);
let start_time = std::time::Instant::now();
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 = mont.primetest_miller_3bases();
let res_lucas = mont.primetest_lucas();
assert!(!res_3bases || res_lucas);
res_3bases != res_lucas
}).collect();
assert_eq!(assumed, result);
eprint!("test_lucas_u64_1e7 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_nzu64_1e7() { // 24bit
let start_time = std::time::Instant::now();
for n in 3..10_000_000 {
if (n & 1) == 0 { continue; }
let nn = NonZeroU64::new(n).unwrap();
let mont = Mont::<NonZeroU64>::new(nn);
let res_3bases = mont.primetest_miller_3bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_3bases, res_bpsw);
}
eprint!("test_bpsw_nzu64_1e7 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_u64_1e7() { // 24bit
let start_time = std::time::Instant::now();
for n in 3..10_000_000 {
if (n & 1) == 0 { continue; }
let mont = Mont::<u64>::new(n);
let res_3bases = mont.primetest_miller_3bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_3bases, res_bpsw);
}
eprint!("test_bpsw_u64_1e7 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_nzu64_4e9() { // 32bit
let start_time = std::time::Instant::now();
for n in 4_000_000_000..4_010_000_000 {
if (n & 1) == 0 { continue; }
let nn = NonZeroU64::new(n).unwrap();
let mont = Mont::<NonZeroU64>::new(nn);
let res_3bases = mont.primetest_miller_3bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_3bases, res_bpsw);
}
eprint!("test_bpsw_nzu64_4e9 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_u64_4e9() { // 32bit
let start_time = std::time::Instant::now();
for n in 4_000_000_000..4_010_000_000 {
if (n & 1) == 0 { continue; }
let mont = Mont::<u64>::new(n);
let res_3bases = mont.primetest_miller_3bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_3bases, res_bpsw);
}
eprint!("test_bpsw_u64_4e9 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_nzu64_1e16() { // 54bit
let start_time = std::time::Instant::now();
for n in 10_000_000_000_000_000..10_000_000_010_000_000 {
if (n & 1) == 0 { continue; }
let nn = NonZeroU64::new(n).unwrap();
let mont = Mont::<NonZeroU64>::new(nn);
let res_7bases = mont.primetest_miller_7bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_7bases, res_bpsw);
}
eprint!("test_bpsw_nzu64_1e16 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_u64_1e16() { // 54bit
let start_time = std::time::Instant::now();
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 = mont.primetest_miller_7bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_7bases, res_bpsw);
}
eprint!("test_bpsw_u64_1e16 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_nzu64_9e18() { // 63bit
let start_time = std::time::Instant::now();
for n in 9_000_000_000_000_000_000..9_000_000_000_010_000_000 {
if (n & 1) == 0 { continue; }
let nn = NonZeroU64::new(n).unwrap();
let mont = Mont::<NonZeroU64>::new(nn);
let res_7bases = mont.primetest_miller_7bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_7bases, res_bpsw);
}
eprint!("test_bpsw_nzu64_9e18 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_u64_9e18() { // 63bit
let start_time = std::time::Instant::now();
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 = mont.primetest_miller_7bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_7bases, res_bpsw);
}
eprint!("test_bpsw_u64_9e18 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_nzu64_10e18() { // 64bit
let start_time = std::time::Instant::now();
for n in 10_000_000_000_000_000_000..10_000_000_000_010_000_000 {
if (n & 1) == 0 { continue; }
let nn = NonZeroU64::new(n).unwrap();
let mont = Mont::<NonZeroU64>::new(nn);
let res_7bases = mont.primetest_miller_7bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_7bases, res_bpsw);
}
eprint!("test_bpsw_nzu64_10e18 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_u64_10e18() { // 64bit
let start_time = std::time::Instant::now();
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 = mont.primetest_miller_7bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_7bases, res_bpsw);
}
eprint!("test_bpsw_u64_10e18 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_nzu64_18e18() { // 64bit
let start_time = std::time::Instant::now();
for n in 18_000_000_000_000_000_000..18_000_000_000_010_000_000 {
if (n & 1) == 0 { continue; }
let nn = NonZeroU64::new(n).unwrap();
let mont = Mont::<NonZeroU64>::new(nn);
let res_7bases = mont.primetest_miller_7bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_7bases, res_bpsw);
}
eprint!("test_bpsw_nzu64_18e18 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_bpsw_u64_18e18() { // 64bit
let start_time = std::time::Instant::now();
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 = mont.primetest_miller_7bases();
let res_bpsw = mont.primetest_bpsw();
assert_eq!(res_7bases, res_bpsw);
}
eprint!("test_bpsw_u64_18e18 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_issq_u64() {
let start_time = std::time::Instant::now();
for i in 0xffff_0000u64..=0xffff_ffff {
let n = i * i;
assert!(issq_u64(n));
}
eprint!("test_issq_u64 {}us\n", start_time.elapsed().as_micros());
}
#[test]
fn test_issq_u64_p1() {
let start_time = std::time::Instant::now();
for i in 0xffff_0000u64..=0xffff_ffff {
let n = i * i + 1;
assert!(!issq_u64(n));
}
eprint!("test_issq_u64_p1 {}us\n", start_time.elapsed().as_micros());
}
}
fn main() {
use std::io::{stdin, stdout, BufReader, BufRead, BufWriter, Write};
let start_time = std::time::Instant::now();
let out = stdout();
let mut writer = BufWriter::new(out.lock());
let inp = stdin();
let mut reader = BufReader::new(inp.lock());
let mut ibuf = Vec::<u8>::new();
reader.read_until(b'\n', &mut ibuf).unwrap();
let tr = &ibuf[0..=((0..(ibuf.len())).rfind(|&i| ibuf[i] > 0x20).unwrap())];
let n = tr.iter().fold(0, |a, &c| a * 10 + (c & 0x0f) as usize);
let mut obuf = Vec::<u8>::with_capacity(n * 22);
for _ in 0..n {
ibuf.clear();
reader.read_until(b'\n', &mut ibuf).unwrap();
let tr = &ibuf[0..=((0..(ibuf.len())).rfind(|&i| ibuf[i] > 0x20).unwrap())];
let x = tr.iter().fold(0, |a, &c| a * 10 + (c & 0x0f) as u64);
let res = if USE_NZU64 {
primetest_nzu64_bpsw(unsafe { NonZeroU64::new_unchecked(x) })
} else {
primetest_u64_bpsw(x)
};
obuf.extend_from_slice(tr);
obuf.extend_from_slice(if res { b" 1\n" } else { b" 0\n" });
}
writer.write_all(obuf.as_slice()).unwrap();
writer.flush().unwrap();
eprintln!("{}us", start_time.elapsed().as_micros());
}