結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー 👑 MizarMizar
提出日時 2023-01-10 03:33:52
言語 Rust
(1.77.0 + proconio)
結果
AC  
実行時間 15 ms / 9,973 ms
コード長 47,510 bytes
コンパイル時間 15,999 ms
コンパイル使用メモリ 390,608 KB
実行使用メモリ 6,820 KB
最終ジャッジ日時 2024-12-21 05:34:22
合計ジャッジ時間 14,606 ms
ジャッジサーバーID
(参考情報)
judge2 / judge3
このコードへのチャレンジ
(要ログイン)

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 1 ms
6,820 KB
testcase_01 AC 1 ms
6,816 KB
testcase_02 AC 1 ms
6,820 KB
testcase_03 AC 1 ms
6,816 KB
testcase_04 AC 11 ms
6,816 KB
testcase_05 AC 12 ms
6,820 KB
testcase_06 AC 10 ms
6,816 KB
testcase_07 AC 10 ms
6,820 KB
testcase_08 AC 10 ms
6,820 KB
testcase_09 AC 15 ms
6,816 KB
権限があれば一括ダウンロードができます

ソースコード

diff #

// -*- coding:utf-8-unix -*-
extern crate core;
// 除算/剰余算インラインセンブリコード切替: inline assembly divrem code switch
macro_rules! cond_div {
    ($nonzero_div:expr, $asm_div:expr) => {
        //$nonzero_div // 1.51.0 or later
        $asm_div // 1.59.0 or later (feature = "asm")
    }
}

// モンゴメリ剰余乗算: Montgomery modular multiplication
pub 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;
}
pub 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;
}
pub trait MontModN<T, BitCountType=u32> {
    fn modn(&self, x: T) -> T;
}
pub 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()
}

macro_rules! impl_mont_val_pr { ($x:ty) => { impl MontVal<$x> for Mont<$x> {
    fn n(&self) -> $x { self.n }
    fn ni(&self) -> $x { self.ni }
    fn nh(&self) -> $x { self.nh }
    fn r(&self) -> $x { self.r }
    fn rn(&self) -> $x { self.rn }
    fn r2(&self) -> $x { self.r2 }
    fn d(&self) -> $x { self.d }
    fn k(&self) -> u32 { self.k }
}}}
impl_mont_val_pr!(u64);
impl_mont_val_pr!(u32);

macro_rules! impl_mont_ops { ($x:ty, $y:ty, $y2:ty, $w:literal) => { impl MontOps<$y> for $x {
    fn add(&self, a: $y, b: $y) -> $y {
        // == a + b (mod n)
        debug_assert!(a < self.n());
        debug_assert!(b < self.n());
        let (t, f) = a.overflowing_sub(self.n.wrapping_sub(b));
        t.wrapping_add((f as $y).wrapping_neg() & self.n())
    }
    fn sub(&self, a: $y, b: $y) -> $y {
        // == a - b (mod n)
        debug_assert!(a < self.n());
        debug_assert!(b < self.n());
        let (t, f) = a.overflowing_sub(b);
        t.wrapping_add((f as $y).wrapping_neg() & self.n())
    }
    fn div2(&self, ar: $y) -> $y {
        // == ar / 2 (mod n)
        debug_assert!(ar < self.n());
        let t = ar >> 1;
        t.wrapping_add((ar & 1).wrapping_neg() & self.nh())
    }
    fn mrmul(&self, ar: $y, br: $y) -> $y {
        // == (ar * br) / r (mod n)
        debug_assert!(ar < self.n());
        debug_assert!(br < self.n());
        let (n, ni) = (self.n(), self.ni());
        let t: $y2 = (ar as $y2) * (br as $y2);
        let (t, f) = ((t >> $w) as $y).overflowing_sub((((((t as $y).wrapping_mul(ni)) as $y2) * (n as $y2)) >> $w) as $y);
        t.wrapping_add((f as $y).wrapping_neg() & n)
    }
    fn mr(&self, ar: $y) -> $y {
        // == ar / r (mod n)
        debug_assert!(ar < self.n());
        let (n, ni) = (self.n(), self.ni());
        let (t, f) = (((((ar.wrapping_mul(ni)) as $y2) * (n as $y2)) >> $w) as $y).overflowing_neg();
        t.wrapping_add((f as $y).wrapping_neg() & n)
    }
    fn ar(&self, a: $y) -> $y {
        // == a * r (mod n)
        debug_assert!(a < self.n());
        self.mrmul(a, self.r2())
    }
    fn pow(&self, mut ar: $y, mut b: $y) -> $y {
        // == ((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); }
        }
    }
    fn powodd(&self, mut ar: $y, mut b: $y) -> $y {
        // == ((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); }
        }
    }
}}}
impl_mont_ops!(Mont<u64>, u64, u128, 64);
impl_mont_ops!(Mont<u32>, u32, u64, 32);

#[allow(unused)]
fn asm_div32(a: u32, b: u32) -> u32 { // a / b (32bit)
    debug_assert_ne!(b, 0);
    cond_div!(
        unsafe { a / std::num::NonZeroU32::new_unchecked(b) },
        if cfg!(target_arch = "x86_64") {
            let mut quot: u32;
            unsafe{core::arch::asm!(
                "div {0:e}",
                in(reg) b,
                inout("eax") a => quot,
                inout("edx") 0u32 => _,
                options(pure, nomem, nostack)
            )}
            debug_assert_eq!(a / b, quot);
            quot
        } else {
            unsafe { a / std::num::NonZeroU32::new_unchecked(b) }
        }
    )
}
#[allow(unused)]
fn asm_rem32(a: u32, b: u32) -> u32 { // a mod b (32bit)
    debug_assert_ne!(b, 0);
    cond_div!(
        unsafe { a % std::num::NonZeroU32::new_unchecked(b) },
        if cfg!(target_arch = "x86_64") {
            let mut rem: u32;
            unsafe{core::arch::asm!(
                "div {0:e}",
                in(reg) b,
                inout("eax") a => _,
                inout("edx") 0u32 => rem,
                options(pure, nomem, nostack)
            )}
            debug_assert_eq!(a % b, rem);
            rem
        } else {
            unsafe { a % std::num::NonZeroU32::new_unchecked(b) }
        }
    )
}
#[allow(unused)]
fn asm_div64(a: u64, b: u64) -> u64 { // a / b (64bit)
    debug_assert_ne!(b, 0);
    cond_div!(
        unsafe { a / std::num::NonZeroU64::new_unchecked(b) },
        if cfg!(target_arch = "x86_64") {
            let mut quot: u64;
            unsafe{core::arch::asm!(
                "div {0}",
                in(reg) b,
                inout("rax") a => quot,
                inout("rdx") 0u64 => _,
                options(pure, nomem, nostack)
            )}
            debug_assert_eq!(a / b, quot);
            quot
        } else {
            unsafe { a / std::num::NonZeroU64::new_unchecked(b) }
        }
    )
}
#[allow(unused)]
fn asm_rem64(a: u64, b: u64) -> u64 { // a mod b (64bit)
    debug_assert_ne!(b, 0);
    cond_div!(
        unsafe { a % std::num::NonZeroU64::new_unchecked(b) },
        if cfg!(target_arch = "x86_64") {
            let mut rem: u64;
            unsafe{core::arch::asm!(
                "div {0}",
                in(reg) b,
                inout("rax") a => _,
                inout("rdx") 0u64 => rem,
                options(pure, nomem, nostack)
            )}
            debug_assert_eq!(a % b, rem);
            rem
        } else {
            unsafe { a % std::num::NonZeroU64::new_unchecked(b) }
        }
    )
}
#[allow(unused)]
fn asm_r2_32(n: u32) -> (u32, u32) { // 2**32 mod n, 2**64 mod n (32bit)
    debug_assert!(n > 1 && (n & 1) == 1); // cause throw when n <= 1 (zero divide or quot overflow)
    cond_div!(
        ((1u32.wrapping_neg() % n).wrapping_add(1), ((1u64.wrapping_neg() % (n as u64)) as u32).wrapping_add(1)),
        if cfg!(target_arch = "x86_64") {
            let (mut r1, mut r2): (u32, u32);
            unsafe{core::arch::asm!(
                "xor eax, eax",
                "mov edx, 1",
                "div {0:e}",
                "mov {1:e}, edx",
                "xor eax, eax",
                "div {0:e}",
                in(reg) n,
                out(reg) r1,
                out("eax") _,
                out("edx") r2,
                options(pure, nomem, nostack)
            )}
            debug_assert_eq!((1u32.wrapping_neg() % n).wrapping_add(1), r1);
            debug_assert_eq!(((1u64.wrapping_neg() % (n as u64)) as u32).wrapping_add(1), r2);
            (r1, r2)
        } else { unsafe {(
            (1u32.wrapping_neg() % std::num::NonZeroU32::new_unchecked(n)).wrapping_add(1),
            ((1u64.wrapping_neg() % std::num::NonZeroU64::new_unchecked(n as u64)) as u32).wrapping_add(1)
        )}}
    )
}
#[allow(unused)]
fn asm_r2_64(n: u64) -> (u64, u64) { // 2**64 mod n, 2**128 mod n (64bit)
    debug_assert!(n > 1 && (n & 1) == 1); // cause throw when n <= 1 (zero divide or quot overflow)
    cond_div!(
        unsafe{(
            (1u64.wrapping_neg() % std::num::NonZeroU64::new_unchecked(n)).wrapping_add(1),
            ((1u128.wrapping_neg() % std::num::NonZeroU128::new_unchecked(n as u128)) as u64).wrapping_add(1)
        )},
        if cfg!(target_arch = "x86_64") {
            let (mut r1, mut r2): (u64, u64);
            unsafe{core::arch::asm!(
                "xor rax, rax",
                "mov rdx, 1",
                "div {0}",
                "mov {1}, rdx",
                "xor rax, rax",
                "div {0}",
                in(reg) n,
                out(reg) r1,
                out("rax") _,
                out("rdx") r2,
                options(pure, nomem, nostack)
            )}
            debug_assert_eq!((1u64.wrapping_neg() % n).wrapping_add(1), r1);
            debug_assert_eq!(((1u128.wrapping_neg() % (n as u128)) as u64).wrapping_add(1), r2);
            (r1, r2)
        } else { unsafe {(
            (1u64.wrapping_neg() % std::num::NonZeroU64::new_unchecked(n)).wrapping_add(1),
            ((1u128.wrapping_neg() % std::num::NonZeroU128::new_unchecked(n as u128)) as u64).wrapping_add(1)
        )}}
    )
}

macro_rules! impl_mont_modn { ($x:ty) => { impl MontModN<$x> for Mont<$x> {
    fn modn(&self, x: $x) -> $x { if x < self.n() { x } else {
        cond_div!(
            match stringify!($x) {
                "u64" => unsafe { (x as u64 % std::num::NonZeroU64::new_unchecked(self.n() as u64)) as $x },
                "u32" => unsafe { (x as u32 % std::num::NonZeroU32::new_unchecked(self.n() as u32)) as $x },
                _ => unimplemented!(),
            },
            match stringify!($x) {
                "u64" => asm_rem64(x as u64, self.n() as u64) as $x,
                "u32" => asm_rem32(x as u32, self.n() as u32) as $x,
                _ => unimplemented!(),
            }
        )
    }}
}}}
impl_mont_modn!(u64);
impl_mont_modn!(u32);

macro_rules! impl_mont_new { ($x:ty, $x2:ty) => { impl MontNew<$x> for Mont<$x> {
    fn new(n: $x) -> Self {
        debug_assert!(n > 1 && (n & 1) == 1);
        // // n is odd number, n = 2*k+1, n >= 1, n < 2**BITS, k is non-negative integer, k >= 0, k < 2**(BITS-1)
        // // As (k*(k+1)) is even, (k*(k+1)/2) is an integer.
        // ni0 := n; // = 2*k+1 = (1+(2**3)*((k*(k+1)/2)**1))/n
        /*
        // ni1 := ni0 * (2 - (n * ni0)); // = (1+(2**6)*((k*(k+1)/2)**2))/(2*k+1)
        // ni2 := ni0 * (2 - (n * ni1)); // = (1+(2**12)*((k*(k+1)/2)**4))/(2*k+1)
        // ni3 := ni0 * (2 - (n * ni2)); // = (1+(2**24)*((k*(k+1)/2)**8))/(2*k+1)
        // ni4 := ni0 * (2 - (n * ni3)); // = (1+(2**48)*((k*(k+1)/2)**16))/(2*k+1)
        // ni5 := ni0 * (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
        let mut ni = n;
        for _ in 0..5 {
            ni = ni.wrapping_mul((2 as $x).wrapping_sub(n.wrapping_mul(ni)));
        }
        // */
        /*
        // 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
        let mut ni = n;
        for _ in 0..3 {
            let t = n.wrapping_mul(ni);
            ni = t.wrapping_sub(3).wrapping_mul(t).wrapping_add(3).wrapping_mul(ni);
        }
        // */
        //*
        let n32 = n as u32;
        let ni = n32;
        let t = n32.wrapping_mul(ni);
        let ni = t.wrapping_sub(3).wrapping_mul(t).wrapping_add(3).wrapping_mul(ni);
        let t = n32.wrapping_mul(ni);
        let ni = t.wrapping_sub(3).wrapping_mul(t).wrapping_add(3).wrapping_mul(ni) as $x;
        let t = n.wrapping_mul(ni);
        let ni = t.wrapping_sub(3).wrapping_mul(t).wrapping_add(3).wrapping_mul(ni) as $x;
        // */
        debug_assert_eq!(n.wrapping_mul(ni), 1); // n * ni == 1 (mod 2**BITS)
        let nh = (n >> 1) + 1; // == (n + 1) / 2 == floor(n / 2) + 1  { n is odd }
        let (r, r2) = cond_div!( // == ( 2**BITS, (2**BITS)**2 ) (mod n)
            (
                (((1 as $x).wrapping_neg() % n).wrapping_add(1)),
                ((((1 as $x2).wrapping_neg() % (n as $x2)) as $x).wrapping_add(1)),
            ),
            match stringify!($x) {
                "u64" => { let (r1, r2) = asm_r2_64(n as u64); (r1 as $x, r2 as $x) },
                "u32" => { let (r1, r2) = asm_r2_32(n as u32); (r1 as $x, r2 as $x) },
                _ => unimplemented!(),
            }
        );
        let rn = n - r;
        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 }
    }
}}}
impl_mont_new!(u64, u128);
impl_mont_new!(u32, u64);

macro_rules! fn_isqrt_n { ($sqrt_newton:ident, $issq_newton:ident, $issq_mod32:ident, $issq_mod3960:ident, $issq_mod4095:ident, $issq:ident, $x:ty, $w:literal) => {
// 擬平方数判定: pseudo-square-number determination (mod32 7/32 0.2188) (true: Could be a square number, false: Not a square number)
fn $issq_mod32(x: $x) -> bool { // The detection rate for odd square number is similar to mod8 (only x mod 8 == 1)
    (0x02030213u32 >> ((x as u32) & 31)) & 1 == 1
}
#[allow(unused)]
// 擬平方数判定: pseudo-square-number determination (mod3960: 216/3960 0.0545) (true: Could be a square number, false: Not a square number)
fn $issq_mod3960(x: $x) -> bool {
    const SQTABLE_MOD3960: [u64; 62] = [0x2001002010213,0x1200001000020001,0x10020000010100,0x210000010,0x201000201,0x120000000010,0x200000000011002,0x101001000001000,0x110020100,0x2010000000001,0x1001000200,0x100000002110000,0x210010000,0x201000001000200,0x110000100010,0x2,0x200001000000001,0x100000110020002,0x10000010000010,0x1000000000000300,0x100120000010000,0x10002000002,0x1001200001000,0x2100000020000,0x3010000000010,0x100201,0x100020000120002,0x10001000000000,0x201000001021000,0x2000100,0x1002000002010000,0x1000020100001001,0x100000100000000,0x3000000001010,0x1000000201000001,0x20000100000,0x10000000001,0x200000200021100,0x10100100000100,0x1002000011000010,0x21000000000,0x20100120000,0x10000000011002,0x200000001001200,0x100000100000,0x2011002010001,0x1200000000000000,0x100020000,0x10000211000000,0x1000201000001,0x102000000000010,0x10010002001002,0x1000000000,0x100100000100,0x10002000200,0x1000101000000201,0x2030000,0x10000000010000,0x1000300000200,0x100000100100,0x2000010,0x1000];
    let p = (x % 3960) as usize;
    (SQTABLE_MOD3960[p >> 6] >> (p & 63)) & 1 == 1
}
// 擬平方数判定: pseudo-square-number determination (mod4095 336/4095 0.0821) (true: Could be a square number, false: Not a square number)
fn $issq_mod4095(x: $x) -> 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
}
// ニュートン法による整数平方根: integer square-root by Newton's method
fn $sqrt_newton(x: $x) -> $x { // == floor(sqrt(x))
    if x <= 1 { return x; }
    let k = $w - ((x - 1).leading_zeros() >> 1);
    let mut s = (1 as $x) << k; // s = 2**k
    let mut t = (s + (x >> k)) >> 1; // t = (s + x/s)/2
    // while loop count (= divide count) may be { u64: max 6 times, u32: max 5 times }
    // s > floor(sqrt(x)) -> floor(sqrt(x)) <= t < s
    // s == floor(sqrt(x)) -> s == floor(sqrt(x)) <= t <= floor(sqrt(x)) + 1
    while t < s {
        s = t;
        t = (s + match stringify!($x) {
            "u64" => asm_div64(x as u64, s as u64) as $x,
            "u32" => asm_div32(x as u32, s as u32) as $x,
            _ => unimplemented!(),
        }) >> 1;
    }
    s
}
// 平方数判定(ニュートン法): square-number determination (Newton's method) (true: Its a square number, false: Not a square number)
fn $issq_newton(x: $x) -> bool {
    let sqrt = $sqrt_newton(x); // == floor(sqrt(x))
    sqrt * sqrt == x
}
// 平方数判定(複合): square-number determination (combined) (true: Its a square number, false: Not a square number)
fn $issq(x: $x) -> bool {
    $issq_mod32(x) && // 擬平方数判定: pseudo-square-number determination (mod32 7/32 0.2188)
    //$issq_mod3960(x) && // 擬平方数判定: pseudo-square-number determination (mod3960: 216/3960 0.0545)
    $issq_mod4095(x) && // 擬平方数判定: pseudo-square-number determination (mod4095 336/4095 0.0821)
    $issq_newton(x) // 平方数判定(ニュートン法): square-number determination (Newton's method)
}
}}
fn_isqrt_n!(isqrt_64_newton, issq_u64_newton, issq_u64_mod32, issq_u64_mod3960, issq_u64_mod4095, issq_u64, u64, 32);
fn_isqrt_n!(isqrt_32_newton, issq_u32_newton, issq_u32_mod32, issq_u32_mod3960, issq_u32_mod4095, issq_u32, u32, 16);

macro_rules! fn_jacobi { ($jacobi:ident, $u:ty, $i:ty, $($x:literal),*) => {
// ヤコビ記号: Jacobi symbol
fn $jacobi(ai: $i, mut n: $u) -> i32 {
    debug_assert!((n & 1) == 1 && ai > <$i>::MIN);
    let (mut a, mut j) = (ai.abs() as $u, if ai < 0 && (n & 3) == 3 { -1 } else { 1 });
    while a != 0 {
        let ba = a.trailing_zeros();
        if (ba & 1) != 0 && (n.wrapping_add(2) & 5) == 5 { j = -j; }
        a >>= ba;
        if (a & n & 3) == 3 { j = -j; }
        n = match a {
            $($x => n % $x,)*
            _ => match stringify!($u) {
                "u64" => asm_rem64(n as u64, a as u64) as $u,
                "u32" => asm_rem32(n as u32, a as u32) as $u,
                _ => n % a,
            },
        };
        std::mem::swap(&mut n, &mut a);
        /*
        if a > (n >> 1) {
            a = n - a;
            if (n & 3) == 3 { j = -j; }
        }
        // */
    }
    if n == 1 { j } else { 0 }
}}}
fn_jacobi!(jacobi_64, u64, i64, 1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31);
fn_jacobi!(jacobi_32, u32, i32, 1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31);

macro_rules! trait_primetest { ($tr:ident, $jacobi:ident, $issq:ident, $u:ty, $i:ty, $w:literal) => {
pub trait $tr<T> : MontOps<$u> + MontVal<$u> + MontModN<$u> {
    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(r, r), d);
        if br == r || br == rn { return true; }
        for _ in 1..k {
            br = self.mrmul(br, br);
            if br == rn { return true; }
        }
        false
    }
    fn primetest_miller(&self, mut base: $u) -> 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
    }
    fn primetest_lucas_dprobe(&self) -> Option<$i> {
        let n = self.n();
        let mut d = 5 as $i;
        for i in 0u32.. {
            debug_assert!(i < 64);
            match $jacobi(d, n) {
                -1 => { return Some(d); },
                0 => if ((d.abs()) as $u) < n { return None; },
                _ => {},
            }
            if i == 32 && $issq(n) { return None; }
            if (i & 1) == 1 { d = 2 - d; } else { d = -(d + 2); }
        }
        None
    }
    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 d = match self.primetest_lucas_dprobe() {
            Some(d) => d,
            None => { return false; },
        };
        let qm = self.ar(self.modn(if d < 0 {((1 - d) as $u) / 4} else {n - ((d - 1) as $u) / 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 = self.ar(if d < 0 { n - self.modn((-d) as $u) } else { self.modn(d as $u) });
        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 >> ($w - 1)) != 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
    }
    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()
    }
}}}
trait_primetest!(PrimeTestU64Trait, jacobi_64, issq_u64, u64, i64, 64);
trait_primetest!(PrimeTestU32Trait, jacobi_32, issq_u32, u32, i32, 32);
impl PrimeTestU64Trait<u64> for Mont<u64> {}
impl PrimeTestU32Trait<u32> for Mont<u32> {}

// Deterministic variants of the Miller-Rabin primality test http://miller-rabin.appspot.com/
pub 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,
    ];
    fn primetest_miller_1bases(&self) -> bool {
        debug_assert!(self.n() < 341531);
        Self::MILLER_U64_1BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_2bases(&self) -> bool {
        debug_assert!(self.n() < 1050535501);
        Self::MILLER_U64_2BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_3bases(&self) -> bool {
        debug_assert!(self.n() < 350269456337);
        Self::MILLER_U64_3BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_4bases(&self) -> bool {
        debug_assert!(self.n() < 55245642489451);
        Self::MILLER_U64_4BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_5bases(&self) -> bool {
        debug_assert!(self.n() < 7999252175582851);
        Self::MILLER_U64_5BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_6bases(&self) -> bool {
        debug_assert!(self.n() < 585226005592931977);
        Self::MILLER_U64_6BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_7bases(&self) -> bool {
        Self::MILLER_U64_7BASES.iter().all(|&base| self.primetest_miller(base))
    }
    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<u64> for Mont<u64> {}

// Deterministic variants of the Miller-Rabin primality test http://miller-rabin.appspot.com/
pub trait MillerU32<T> : PrimeTestU32Trait<T> {
    const MILLER_U32_1BASES: [u32; 1] = [921211727];
    const MILLER_U32_2BASES: [u32; 2] = [1143370,2350307676];
    const MILLER_U32_3BASES: [u32; 3] = [2,7,61];
    fn primetest_miller_1bases(&self) -> bool {
        debug_assert!(self.n() < 49141);
        Self::MILLER_U32_1BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_2bases(&self) -> bool {
        debug_assert!(self.n() < 360018361);
        Self::MILLER_U32_2BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_3bases(&self) -> bool {
        Self::MILLER_U32_3BASES.iter().all(|&base| self.primetest_miller(base))
    }
    fn primetest_miller_mbases(&self) -> bool {
        // Deterministic variants of the Miller-Rabin primality test http://miller-rabin.appspot.com/
        match self.n() {
            0..=49140     => Self::MILLER_U32_1BASES.iter(),
            0..=360018360 => Self::MILLER_U32_2BASES.iter(),
            _             => Self::MILLER_U32_3BASES.iter(),
        }.all(|&base| self.primetest_miller(base))
    }
}
impl MillerU32<u32> for Mont<u32> {}

// 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()
}
// Baillie–PSW primarity test (u32)
pub fn primetest_u32_bpsw(n: u32) -> bool {
    if n == 2 { return true; }
    if n == 1 || (n & 1) == 0 { return false; }
    let mont = Mont::<u32>::new(n);
    mont.primetest_bpsw()
}
// 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()
}
// Miller-Rabin primarity test (u32)
pub fn primetest_u32_miller_mbases(n: u32) -> bool {
    if n == 2 { return true; }
    if n == 1 || (n & 1) == 0 { return false; }
    let mont = Mont::<u32>::new(n);
    mont.primetest_miller_mbases()
}

#[cfg(test)]
mod tests {
    use crate::*;
    enum MontTestType { U32, U64, }

    // 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]),
    ];
    // #SPSP-2 Miller-Rabin base 2 (up to 1e7)
    // 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 )
    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
    ];
    // #SLPSP Strong Lucas-Selfridge (up to 1e7)
    // 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 )
    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
    ];

    fn exec_eq(n: u64, f: bool, ty: &MontTestType) {
        if n < 3 || (n & 1) == 0 { return; }
        match ty {
            MontTestType::U64 => {
                let mont = Mont::<u64>::new(n);
                // Deterministic variants of the Miller-Rabin primality test
                // http://miller-rabin.appspot.com/
                if n < 341531 { assert_eq!(mont.primetest_miller_1bases(), f); }
                if n < 1050535501 { assert_eq!(mont.primetest_miller_2bases(), f); }
                if n < 350269456337 { assert_eq!(mont.primetest_miller_3bases(), f); }
                if n < 55245642489451 { assert_eq!(mont.primetest_miller_4bases(), f); }
                if n < 7999252175582851 { assert_eq!(mont.primetest_miller_5bases(), f); }
                if n < 585226005592931977 { assert_eq!(mont.primetest_miller_6bases(), f); }
                assert_eq!(mont.primetest_miller_7bases(), f);
                assert_eq!(mont.primetest_miller_mbases(), f);
                // Baillie–PSW primality test
                assert_eq!(mont.primetest_bpsw(), f);
            },
            MontTestType::U32 => {
                if n > 0xffff_ffff { return; }
                let mont = Mont::<u32>::new(n as u32);
                // Deterministic variants of the Miller-Rabin primality test
                if n < 49140     { assert_eq!(mont.primetest_miller_1bases(), f); }
                if n < 360018360 { assert_eq!(mont.primetest_miller_2bases(), f); }
                assert_eq!(mont.primetest_miller_3bases(), f);
                assert_eq!(mont.primetest_miller_mbases(), f);
                // Baillie–PSW primality test
                assert_eq!(mont.primetest_bpsw(), f);
            },
        }
    }
    fn gen_nbits(ty: &MontTestType) {
        // ten least k's for which (2**n)-k is prime
        // https://primes.utm.edu/lists/2small/0bit.html
        for (bit_ref, kvec) in PRIMES_LAST10.iter() {
            let bit = *bit_ref;
            if bit > 32 {
                match ty {
                    MontTestType::U32 => continue,
                    _ => {},
                }
            }
            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();
                exec_eq(n, f, ty);
            }
        }
    }
    #[test]
    fn nbits_u64() {
        let start_time = std::time::Instant::now();
        gen_nbits(&MontTestType::U64);
        eprint!("nbits_u64 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn nbits_u32() {
        let start_time = std::time::Instant::now();
        gen_nbits(&MontTestType::U32);
        eprint!("nbits_u32 {}us\n", start_time.elapsed().as_micros());
    }

    fn exec_base2(n: u64, ty: &MontTestType) -> (bool, bool) {
        match ty {
            MontTestType::U64 => {
                let mont = Mont::<u64>::new(n);
                let res_base2 = mont.primetest_base2();
                let res_mbases = mont.primetest_miller_mbases();
                (res_base2, res_mbases)
            },
            MontTestType::U32 => {
                let mont = Mont::<u32>::new(n as u32);
                let res_base2 = mont.primetest_base2();
                let res_mbases = mont.primetest_miller_mbases();
                (res_base2, res_mbases)
            },
        }
    }
    fn gen_base2_1e7(ty: &MontTestType) {
        // 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 result: Vec<u64> = (3..10_000_000).filter(|&n| {
            if (n & 1) == 0 { return false; }
            let (res_base2, res_mbases) = exec_base2(n, ty);
            assert!(!res_mbases || res_base2);
            res_mbases != res_base2
        }).collect();
        assert_eq!(Vec::from(SPRP_2), result);
    }
    #[test]
    fn base2_u64_1e7() {
        let start_time = std::time::Instant::now();
        gen_base2_1e7(&MontTestType::U64);
        eprint!("base2_u64_1e7 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn base2_u32_1e7() {
        let start_time = std::time::Instant::now();
        gen_base2_1e7(&MontTestType::U32);
        eprint!("base2_u32_1e7 {}us\n", start_time.elapsed().as_micros());
    }

    fn exec_lucas(n: u64, ty: &MontTestType) -> (bool, bool) {
        match ty {
            MontTestType::U64 => {
                let mont = Mont::<u64>::new(n);
                let res_lucas = mont.primetest_lucas();
                let res_mbases = mont.primetest_miller_mbases();
                (res_lucas, res_mbases)
            },
            MontTestType::U32 => {
                let mont = Mont::<u32>::new(n as u32);
                let res_lucas = mont.primetest_lucas();
                let res_mbases = mont.primetest_miller_mbases();
                (res_lucas, res_mbases)
            },
        }
    }
    fn gen_lucas_1e7(ty: &MontTestType) {
        // 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 result: Vec<u64> = (3..10_000_000).filter(|&n| {
            if (n & 1) == 0 { return false; }
            let (res_base2, res_mbases) = exec_lucas(n, ty);
            assert!(!res_mbases || res_base2);
            res_mbases != res_base2
        }).collect();
        assert_eq!(Vec::from(SLPRP), result);
    }
    #[test]
    fn lucas_u64_1e7() {
        let start_time = std::time::Instant::now();
        gen_lucas_1e7(&MontTestType::U64);
        eprint!("lucas_u64_1e7 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn lucas_u32_1e7() {
        let start_time = std::time::Instant::now();
        gen_lucas_1e7(&MontTestType::U32);
        eprint!("lucas_u32_1e7 {}us\n", start_time.elapsed().as_micros());
    }

    fn exec_bpsw(n: u64, ty: &MontTestType) {
        match ty {
            MontTestType::U64 => {
                let mont = Mont::<u64>::new(n);
                let res_bpsw = mont.primetest_bpsw();
                let res_mbases = mont.primetest_miller_mbases();
                assert_eq!(res_mbases, res_bpsw);
            },
            MontTestType::U32 => {
                let mont = Mont::<u32>::new(n as u32);
                let res_bpsw = mont.primetest_bpsw();
                let res_mbases = mont.primetest_miller_mbases();
                assert_eq!(res_mbases, res_bpsw);
            },
        }
    }
    fn gen_bpsw(minn: u64, maxn: u64, ty: &MontTestType) {
        // Baillie–PSW primarity test
        for n in minn..maxn {
            if n < 3 || (n & 1) == 0 { continue; }
            exec_bpsw(n, ty);
        }
    }
    #[test]
    fn bpsw_u64_1e7() { // 24bit
        let start_time = std::time::Instant::now();
        gen_bpsw(3, 10_000_000, &MontTestType::U64);
        eprint!("bpsw_u64_1e7 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn bpsw_u32_1e7() { // 24bit
        let start_time = std::time::Instant::now();
        gen_bpsw(3, 10_000_000, &MontTestType::U32);
        eprint!("bpsw_u32_1e7 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn bpsw_u64_4e9() { // 32bit
        let start_time = std::time::Instant::now();
        gen_bpsw(4_000_000_000, 4_010_000_000, &MontTestType::U64);
        eprint!("bpsw_u64_4e9 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn bpsw_u32_4e9() { // 32bit
        let start_time = std::time::Instant::now();
        gen_bpsw(4_000_000_000, 4_010_000_000, &MontTestType::U32);
        eprint!("bpsw_u32_4e9 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn bpsw_u64_1e16() { // 54bit
        let start_time = std::time::Instant::now();
        gen_bpsw(10_000_000_000_000_000, 10_000_000_010_000_000, &MontTestType::U64);
        eprint!("bpsw_u64_1e16 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn bpsw_u64_9e18() { // 63bit
        let start_time = std::time::Instant::now();
        gen_bpsw(9_000_000_000_000_000_000, 9_000_000_000_010_000_000, &MontTestType::U64);
        eprint!("bpsw_u64_9e18 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn bpsw_u64_10e18() { // 64bit
        let start_time = std::time::Instant::now();
        gen_bpsw(10_000_000_000_000_000_000, 10_000_000_000_010_000_000, &MontTestType::U64);
        eprint!("bpsw_u64_10e18 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn bpsw_u64_18e18() { // 64bit
        let start_time = std::time::Instant::now();
        gen_bpsw(18_000_000_000_000_000_000, 18_000_000_000_010_000_000, &MontTestType::U64);
        eprint!("bpsw_u64_18e18 {}us\n", start_time.elapsed().as_micros());
    }

    #[test]
    fn sq_u64() {
        let start_time = std::time::Instant::now();
        for i in 0xff00_0000u64..=0xffff_ffff {
            let n = i * i;
            assert!(issq_u64(n));
        }
        eprint!("sq_u64 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn sq_u64_p1() {
        let start_time = std::time::Instant::now();
        for i in 0xff00_0000u64..=0xffff_ffff {
            let n = i * i + 1;
            assert!(!issq_u64(n));
        }
        eprint!("sq_u64_p1 {}us\n", start_time.elapsed().as_micros());
    }
    #[test]
    fn sq_u64_m1() {
        let start_time = std::time::Instant::now();
        for i in 0xff00_0000u64..=0xffff_ffff {
            let n = i * i - 1;
            assert!(!issq_u64(n));
        }
        eprint!("sq_u64_m1 {}us\n", start_time.elapsed().as_micros());
    }

}

pub fn main() {
    use std::io::{stdin, stdout, Read, Write};
    let start_time = std::time::Instant::now();
    let ibuf = {
        let stdin = stdin();
        let mut handle = stdin.lock();
        let mut ibuf = Vec::<u8>::new();
        let _len = handle.read_to_end(&mut ibuf).unwrap();
        ibuf
    };
    let mut it = ibuf.iter().cloned().enumerate();
    let n = {
        let mut a = if let Some((_, c @ b'0'..=b'9')) = it.next() { (c & 0x0f) as usize } else { unreachable!() };
        while let Some((_, c @ b'0'..=b'9')) = it.next() { a = a * 10 + (c & 0x0f) as usize; }
        a
    };
    let xv: Vec<(u64, &[u8])> = (0..n).map(|_| {
        let (i, mut a) = loop { if let Some((i, c @ b'0'..=b'9')) = it.next() { break (i, (c & 0x0f) as u64) }};
        loop { match it.next() {
            Some((_, c @ b'0'..=b'9')) => { a = a * 10 + (c & 0x0f) as u64; },
            Some((j, _)) => break (a, &ibuf[i..j]),
            None => break (a, &ibuf[i..]),
        }}
    }).collect();
    let dur_input = start_time.elapsed();
    let xv: Vec<(&[u8], bool)> = xv.iter().cloned().map(|(x, tr)| {
        (tr, primetest_u64_bpsw(x))
    }).collect();
    let dur_calc = start_time.elapsed();
    {
        let mut obin = Vec::<u8>::with_capacity(n * 22);
        for (tr, res) in xv.iter().cloned() {
            obin.extend_from_slice(tr);
            obin.extend_from_slice(if res { b" 1\n" } else { b" 0\n" });
        }
        let mut out = stdout();
        out.write_all(obin.as_slice()).unwrap();
    }
    let dur_output = start_time.elapsed();
    eprint!(
        "{}us\n{}us\n{}us\n",
        dur_input.as_micros(),
        dur_calc.as_micros(),
        dur_output.as_micros(),
    );
}
0