結果

問題 No.8030 ミラー・ラビン素数判定法のテスト
ユーザー sakikuroe
提出日時 2024-12-09 22:40:16
言語 Rust
(1.83.0 + proconio)
結果
AC  
実行時間 59 ms / 9,973 ms
コード長 4,806 bytes
コンパイル時間 11,600 ms
コンパイル使用メモリ 402,284 KB
実行使用メモリ 6,820 KB
最終ジャッジ日時 2024-12-09 22:40:31
合計ジャッジ時間 13,187 ms
ジャッジサーバーID
（参考情報） 		judge4 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 10
権限があれば一括ダウンロードができます

ソースコード

diff # 

pub struct MontgomeryU64 {
    // r = 2 ** 64
    modulo: u64,       // odd and (m < 2 ** 63)
    modulo_prime: u64, // modulo * modulo_prime = -1 mod r
    r: u64,            // r mod m
    r2: u64,           // r^{2} mod m
    r3: u64,           // r^{3} mod m
}

impl MontgomeryU64 {
    /// Returns
    ///     (x, y, g)
    ///         s.t. a * x + b * y = g
    ///              g = gcd(a, b)
    fn extended_gcd(a: i128, b: i128) -> (i128, i128, u128) {
        let (mut xs, mut ys, mut s) = (1, 0, a);
        let (mut xt, mut yt, mut t) = (0, 1, b);

        while s % t != 0 {
            let q = s / t;
            let (xu, yu, u) = (xs - q * xt, ys - q * yt, s - q * t);
            (xs, ys, xt, yt) = (xt, yt, xu, yu);
            (s, t) = (t, u);
        }

        if t < 0 {
            (-xt, -yt, (-t) as u128)
        } else {
            (xt, yt, t as u128)
        }
    }

    pub fn new(modulo: u64) -> Self {
        assert!(modulo % 2 == 1);
        assert!(modulo < 1 << 63);

        let modulo_prime = {
            let (_, b, _) = Self::extended_gcd(1_i128 << 64, modulo as i128);
            if b <= 0 {
                (-b) as u64
            } else {
                (-b + (1_i128 << 64)) as u64
            }
        };
        let r = ((1_u128 << 64) % modulo as u128) as u64;
        let r2 = ((r as u128 * r as u128) % modulo as u128) as u64;
        let r3 = ((r2 as u128 * r as u128) % modulo as u128) as u64;

        Self {
            modulo,
            modulo_prime,
            r,
            r2,
            r3,
        }
    }

    /// Returns:
    ///     a * r mod N
    pub fn to_montgomery_form(&self, a: u64) -> u64 {
        self.mul(a, self.r2)
    }

    /// Montgomery reduction
    ///
    /// Returns:
    ///     t * r^{-1} mod N
    pub fn reduction(&self, t: u128) -> u64 {
        let t = {
            ((t + (((t as u64).wrapping_mul(self.modulo_prime)) as u128 * self.modulo as u128))
                >> 64) as u64
        };

        if t < self.modulo {
            t
        } else {
            t - self.modulo
        }
    }

    pub fn modulo(&self) -> u64 {
        self.modulo
    }

    pub fn add(&self, ar: u64, br: u64) -> u64 {
        let t = ar + br;
        if t < self.modulo {
            t
        } else {
            t - self.modulo
        }
    }

    pub fn sub(&self, ar: u64, br: u64) -> u64 {
        let (t, f) = ar.overflowing_sub(br);
        if !f {
            t
        } else {
            t.wrapping_add(self.modulo)
        }
    }

    pub fn mul(&self, ar: u64, br: u64) -> u64 {
        self.reduction(ar as u128 * br as u128)
    }

    pub fn inv(&self, ar: u64) -> Option<u64> {
        let (x, _, g) = Self::extended_gcd(ar as i128, self.modulo as i128);
        if g != 1 {
            None
        } else {
            Some(if x < 0 {
                self.mul((x as i64 + self.modulo as i64) as u64, self.r3)
            } else {
                self.mul(x as u64, self.r3)
            })
        }
    }

    pub fn pow(&self, ar: u64, mut n: usize) -> u64 {
        let mut res = self.r;
        let mut x = ar;

        while n > 0 {
            if n % 2 == 1 {
                res = self.mul(res, x);
            }
            x = self.mul(x, x);
            n /= 2;
        }

        res
    }
}

/// Returns:
///     if n is prime number:
///         true
///     else:
///         false
///
/// Algorithm:
///     Miller-Rabin
///
/// References:
///     - [Deterministic variants of the Miller-Rabin primality test. Miller-Rabin SPRP bases records](https://miller-rabin.appspot.com/)
///     - [64bit数の素数判定](https://zenn.dev/mizar/articles/791698ea860581)
pub fn is_prime(n: u64) -> bool {
    if n == 0 || n == 1 {
        return false;
    }

    if n == 2 {
        return true;
    }

    if n % 2 == 0 {
        return false;
    }

    let s = (n - 1).trailing_zeros();
    let d = (n - 1) >> s;
    let mont = MontgomeryU64::new(n);
    let mont_1 = mont.to_montgomery_form(1);
    let mont_n_1 = mont.to_montgomery_form(n - 1);

    let maybe_prime = |a| {
        let a = a % n;

        if a == 0 {
            return true;
        }

        let mut ad = mont.pow(mont.to_montgomery_form(a), d as usize);

        if ad == mont_1 || ad == mont_n_1 {
            return true;
        }

        for _ in 1..s {
            ad = mont.pow(ad, 2);
            if ad == mont_n_1 {
                return true;
            }
        }

        false
    };

    [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
        .into_iter()
        .all(maybe_prime)
}

use proconio::input;

fn main() {
    input! {
        q: u64,
    }

    for _ in 0..q {
        input! {
            n: u64,
        }

        println!("{} {}", n, if is_prime(n) { 1 } else { 0 });
    }
}
0