// -*- coding:utf-8-unix -*- // モンゴメリ剰余乗算 (Montgomery modular multiplication) pub trait UMontTrait { 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) -> u32; // 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 // 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 // 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; } pub struct U64Mont { n: u64, // n is odd, and n > 2 ni: u64, // n * ni == 1 (mod 2**64) nh: u64, // == (n + 1) / 2 r: u64, // == 2**64 (mod n) rn: u64, // == -(2**64) (mod n) r2: u64, // == 2**128 (mod n) d: u64, // == (n - 1) >> (n - 1).trailing_zeros() k: u32, // == (n - 1).trailing_zeros() } impl UMontTrait for U64Mont { #[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 } #[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**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) // ni5 := ni4 * (2 - (n * ni4)); // = (1-(2**64)*((k*(k+1))**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; // == (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 || !fs { u } else { t } } #[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()); if (ar & 1) == 0 { ar >> 1 } else { (ar >> 1) + 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 { // 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 }; b >>= 1; while b != 0 { ar = self.mrmul(ar, ar); if (b & 1) != 0 { t = self.mrmul(t, ar); } b >>= 1; } t } } // 64bit整数平方根(固定ループ回数) -> (floor(sqrt(iv)), remain) #[inline] 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) + d + (b & e); } else { b = ((b + b) & c) + (b & e); } a <<= 2; } ((b >> t) as u64, (a >> t) as u64) } // 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) as u64, -1) } else { ((-a) as u64, 1) }; while a > 0 { let ba = a.trailing_zeros(); a >>= ba; if ((n & 7) == 3 || (n & 7) == 5) && (ba & 1) == 1 { j = -j; } if (a & n & 3) == 3 { j = -j; } let t = a; a = n; n = t; a %= n; if a > (n >> 1) { a = n - a; if (n & 3) == 3 { j = -j; } } } if n == 1 { j } else { 0 } } #[inline] fn prime_test_miller_sub(mont: &U64Mont, mut base: u64) -> 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 (n, r, d, k) = (mont.n(), mont.r(), mont.d(), mont.k()); if base >= n { base %= n; if base == 0 { return true; } } let mut br = mont.pow(mont.ar(2), d); if br == r || br == mont.rn { return true; } for _ in 1..k { br = mont.mrmul(br, br); if br == mont.rn { return true; } } false } #[inline] fn prime_test_miller2_sub(mont: &U64Mont) -> 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, d, k) = (mont.r(), mont.d(), mont.k()); let mut br = mont.pow(mont.ar(2), d); if br == r || br == mont.rn { return true; } for _ in 1..k { br = mont.mrmul(br, br); if br == mont.rn { return true; } } false } fn prime_test_lucas_sub(mont: &U64Mont) -> 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..64 { if jacobi(d, n) == -1 { break; } if i == 32 && isqrt64f(n).1 == 0 { return false; } if (i & 1) == 1 { d = 2 - d; } else { d = -(d + 2); } } let qm = mont.ar(if d < 0 {((1 - d) as u64) / 4 % n} else {n - ((d - 1) as u64) / 4 % n}); let mut k = (n + 1) << (n + 1).leading_zeros(); let mut um = mont.r(); let mut vm = mont.r(); let mut qn = qm; let dm: u64 = mont.ar(if d < 0 { let nd = ((-d) as u64) % n; if nd == 0 { 0 } else { n - nd } } 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 = (n + 1) & (!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 prime_test_bpsw_sub(mont: &U64Mont) -> 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,... prime_test_miller2_sub(mont) && // Lucas primality test // strong Lucas pseudoprimes ( https://oeis.org/A217255 ): 5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,... prime_test_lucas_sub(mont) } // Baillie–PSW primarity test pub fn prime_test_bpsw(n: u64) -> bool { if n == 2 { return true; } if n == 1 || (n & 1) == 0 { return false; } let mont = U64Mont::new(n); prime_test_bpsw_sub(&mont) } #[inline] fn prime_test_miller_1base(mont: &U64Mont) -> bool { // Deterministic variants of the Miller-Rabin primality test // http://miller-rabin.appspot.com/ assert!(mont.n() <= 341531); prime_test_miller_sub(mont, 9345883071009581737) } #[inline] fn prime_test_miller_2bases(mont: &U64Mont) -> bool { // Deterministic variants of the Miller-Rabin primality test // http://miller-rabin.appspot.com/ assert!(mont.n() <= 1050535501); prime_test_miller_sub(mont, 336781006125) && prime_test_miller_sub(mont, 9639812373923155) } #[inline] fn prime_test_miller_3bases(mont: &U64Mont) -> bool { // Deterministic variants of the Miller-Rabin primality test // http://miller-rabin.appspot.com/ assert!(mont.n() <= 350269456337); prime_test_miller_sub(mont, 4230279247111683200) && prime_test_miller_sub(mont, 14694767155120705706) && prime_test_miller_sub(mont, 16641139526367750375) } #[inline] fn prime_test_miller_4bases(mont: &U64Mont) -> bool { // Deterministic variants of the Miller-Rabin primality test // http://miller-rabin.appspot.com/ assert!(mont.n() <= 55245642489451); prime_test_miller_sub(mont, 2) && prime_test_miller_sub(mont, 141889084524735) && prime_test_miller_sub(mont, 1199124725622454117) && prime_test_miller_sub(mont, 11096072698276303650) } #[inline] fn prime_test_miller_5bases(mont: &U64Mont) -> bool { // Deterministic variants of the Miller-Rabin primality test // http://miller-rabin.appspot.com/ assert!(mont.n() <= 7999252175582851); prime_test_miller_sub(mont, 2) && prime_test_miller_sub(mont, 4130806001517) && prime_test_miller_sub(mont, 149795463772692060) && prime_test_miller_sub(mont, 186635894390467037) && prime_test_miller_sub(mont, 3967304179347715805) } #[inline] fn prime_test_miller_6bases(mont: &U64Mont) -> bool { // Deterministic variants of the Miller-Rabin primality test // http://miller-rabin.appspot.com/ assert!(mont.n() <= 585226005592931977); prime_test_miller_sub(mont, 2) && prime_test_miller_sub(mont, 123635709730000) && prime_test_miller_sub(mont, 9233062284813009) && prime_test_miller_sub(mont, 43835965440333360) && prime_test_miller_sub(mont, 761179012939631437) && prime_test_miller_sub(mont, 1263739024124850375) } #[inline] fn prime_test_miller_7bases(mont: &U64Mont) -> bool { // Deterministic variants of the Miller-Rabin primality test // http://miller-rabin.appspot.com/ prime_test_miller_sub(mont, 2) && prime_test_miller_sub(mont, 325) && prime_test_miller_sub(mont, 9375) && prime_test_miller_sub(mont, 28178) && prime_test_miller_sub(mont, 450775) && prime_test_miller_sub(mont, 9780504) && prime_test_miller_sub(mont, 1795265022) } pub fn prime_test_64_miller(n: u64) -> bool { if n == 2 { return true; } if n == 1 || (n & 1) == 0 { return false; } let mont = U64Mont::new(n); match n { // Deterministic variants of the Miller-Rabin primality test // http://miller-rabin.appspot.com/ 0..=341531 => prime_test_miller_1base(&mont), 0..=1050535501 => prime_test_miller_2bases(&mont), 0..=350269456337 => prime_test_miller_3bases(&mont), 0..=55245642489451 => prime_test_miller_4bases(&mont), 0..=7999252175582851 => prime_test_miller_5bases(&mont), 0..=585226005592931977 => prime_test_miller_6bases(&mont), _ => prime_test_miller_7bases(&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 = lines.take(n).map(|l| l.unwrap().parse().unwrap()).collect(); let elapsed1 = start_time.elapsed(); let res: Vec = 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 prime_test_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)> = 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 = U64Mont::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!(prime_test_miller_1base(&mont), f); } // n <= 1050535501 : 2 bases Miller-Rabin primality test if n <= 1050535501 { assert_eq!(prime_test_miller_2bases(&mont), f); } // n <= 350269456337 : 3 bases Miller-Rabin primality test if n <= 350269456337 { assert_eq!(prime_test_miller_3bases(&mont), f); } // n <= 55245642489451 : 4 bases Miller-Rabin primality test if n <= 55245642489451 { assert_eq!(prime_test_miller_4bases(&mont), f); } // n <= 7999252175582851 : 5 bases Miller-Rabin primality test if n <= 7999252175582851 { assert_eq!(prime_test_miller_5bases(&mont), f); } // n <= 585226005592931977 : 6 bases Miller-Rabin primality test if n <= 585226005592931977 { assert_eq!(prime_test_miller_6bases(&mont), f); } // n for all 64bit integer : 7 bases Miller-Rabin primality test assert_eq!(prime_test_miller_7bases(&mont), f); // Baillie–PSW primality test assert_eq!(prime_test_bpsw_sub(&mont), f); } } } #[test] fn test_base2_1e7() { // 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,... // 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 = 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 = (3..10_000_000).filter(|&n| { if n & 1 == 0 { return false; } let mont = U64Mont::new(n); let res_3bases = prime_test_miller_3bases(&mont); let res_base2 = prime_test_miller2_sub(&mont); assert!(!res_3bases || res_base2); res_3bases != 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 = 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 = (3..10_000_000).filter(|&n| { if n & 1 == 0 { return false; } let mont = U64Mont::new(n); let res_3bases = prime_test_miller_3bases(&mont); let res_lucas = prime_test_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 = U64Mont::new(n); let res_3bases = prime_test_miller_3bases(&mont); let res_bpsw = prime_test_bpsw_sub(&mont); assert_eq!(res_3bases, res_bpsw); } } #[test] fn test_bpsw_4e9() { // 32bit for n in 4_000_000_000..4_010_000_000 { if (n & 1) == 0 { continue; } let mont = U64Mont::new(n); let res_3bases = prime_test_miller_3bases(&mont); let res_bpsw = prime_test_bpsw_sub(&mont); assert_eq!(res_3bases, res_bpsw); } } #[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 = U64Mont::new(n); let res_7bases = prime_test_miller_7bases(&mont); let res_bpsw = prime_test_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 = U64Mont::new(n); let res_7bases = prime_test_miller_7bases(&mont); let res_bpsw = prime_test_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 = U64Mont::new(n); let res_7bases = prime_test_miller_7bases(&mont); let res_bpsw = prime_test_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 = U64Mont::new(n); let res_7bases = prime_test_miller_7bases(&mont); let res_bpsw = prime_test_bpsw_sub(&mont); assert_eq!(res_7bases, res_bpsw); } } }