結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー urectanc
提出日時 2025-09-03 19:42:26
言語 Rust
(1.83.0 + proconio)
結果
AC  
実行時間 28 ms / 9,973 ms
コード長 3,517 bytes
コンパイル時間 11,983 ms
コンパイル使用メモリ 398,316 KB
実行使用メモリ 7,720 KB
最終ジャッジ日時 2025-09-03 19:42:41
合計ジャッジ時間 13,960 ms
ジャッジサーバーID
（参考情報） 		judge4 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 10
権限があれば一括ダウンロードができます

ソースコード

diff # 

use proconio::input;
use std::io::Write;

fn main() {
    input! {
        n: usize,
        x: [u64; n],
    }
    let mut stdout = std::io::BufWriter::new(std::io::stdout().lock());

    for &x in &x {
        let ans = math::is_prime(x);
        writeln!(stdout, "{x} {}", if ans { '1' } else { '0' }).unwrap();
    }
}

#[allow(dead_code)]
mod math {
    pub fn is_prime(n: u64) -> bool {
        const WITNESS_3: [u64; 3] = [2, 7, 61];
        const WITNESS_7: [u64; 7] = [2, 325, 9375, 28178, 450775, 9780504, 1795265022];
        const THRESHOLD: u64 = 4_759_123_141;

        if n == 1 || n % 2 == 0 {
            return n == 2;
        }

        let witness = if n < THRESHOLD {
            &WITNESS_3[..]
        } else {
            &WITNESS_7[..]
        };

        let montgomery = Montgomery::new(n);
        let (one, minus_one) = (montgomery.encode(1), montgomery.encode(n - 1));
        let d = n >> (n - 1).trailing_zeros();

        for &a in witness {
            if n <= a {
                continue;
            }

            let mut d = d;
            let mut y = montgomery.pow(montgomery.encode(a), d);
            while d != n - 1 && y != one && y != minus_one {
                y = montgomery.mul(y, y);
                d <<= 1;
            }
            if y != minus_one && d & 1 == 0 {
                return false;
            }
        }

        true
    }

    pub struct Montgomery {
        n: u64,
        n_inv: u64,
        r2: u64,
    }

    impl Montgomery {
        pub fn new(modulus: u64) -> Self {
            assert_eq!(modulus % 2, 1);
            let n = modulus;
            let mut n_inv = 1u64;
            for _ in 0..6 {
                n_inv = n_inv.wrapping_mul(2u64.wrapping_sub(n.wrapping_mul(n_inv)));
            }
            let r2 = ((n as u128).wrapping_neg() % (n as u128)) as u64;
            Self { n, n_inv, r2 }
        }

        pub fn encode(&self, x: u64) -> u64 {
            self.mul(x, self.r2)
        }

        pub fn reduce(&self, x: u64) -> u64 {
            self.mul(x, 1)
        }

        pub fn add(&self, a: u64, b: u64) -> u64 {
            let (sum, carry) = a.overflowing_add(b);
            let (sub, borrow) = sum.overflowing_sub(self.n);
            if !carry && borrow {
                sum
            } else {
                sub
            }
        }

        pub fn mul(&self, a: u64, b: u64) -> u64 {
            let c = a as u128 * b as u128;
            let d = self.n_inv.wrapping_mul(c as u64);
            let e = ((d as u128 * self.n as u128) >> 64) as u64;
            let (sub, borrow) = ((c >> 64) as u64).overflowing_sub(e);
            if borrow {
                sub.wrapping_add(self.n)
            } else {
                sub
            }
        }

        fn mul_add_unnormalized(&self, a: u64, b: u64, c: u64) -> u64 {
            let d = a as u128 * b as u128;
            let e = ((d >> 64) as u64).wrapping_add(self.n).wrapping_add(c);
            let f = self.n_inv.wrapping_mul(d as u64);
            let g = ((f as u128 * self.n as u128) >> 64) as u64;
            e.wrapping_sub(g)
        }

        pub fn pow(&self, lhs: u64, mut e: u64) -> u64 {
            let one = self.encode(1);
            let mut res = one;
            let mut pow = lhs;
            while e > 0 {
                if e & 1 == 1 {
                    res = self.mul(res, pow);
                }
                pow = self.mul(pow, pow);
                e >>= 1;
            }
            res
        }
    }
}
0