Rustで多倍長整数と素数生成

暗号のお勉強をしているときにRustで何かしらプログラムを書こうと思ったので、Rustでの多倍長整数の取扱についての軽いメモ。以下のソースコードはNバイトの素数をMiller–Rabinの素数判定法を用いて確率的に生成するもの（まともにデバッグしていないので間違えているかも）。多倍長整数演算はnum_bigintを使用。

use num_bigint_dig::BigUint;
use rand::{RngCore, rngs::ThreadRng};
use competitive::prime::*;
use num_traits::identities::Zero;
use num_traits::One;
use num_traits::Pow;

pub fn gen_prime<const N: usize>(r: usize) -> BigUint {
    let primes =  primes(1_000_000);
    let mut rng = ThreadRng::default();
    let mut vv = [0u8; N];
    rng.fill_bytes(&mut vv);
    vv[0] |= 1 << 7;
    vv[N-1] |= 1;
    let mut ret = BigUint::from_bytes_be(&vv);
    loop {
        if primes.iter().copied().all(|p| &ret % p != BigUint::zero()) {
            if miller_rabin::<N>(&ret, r, &mut rng) {
                return ret
            }
        }
        rng.fill_bytes(&mut vv);
        vv[0] |= 1 << 7;
        vv[N-1] |= 1;
        ret = BigUint::from_bytes_be(&vv);
    }
} 

fn miller_rabin<const N: usize>(n: &BigUint, r: usize, rng: &mut ThreadRng) -> bool {
    let n1 = n.clone()-1u32;
    let s = {
        let mut s = 0;
        for i in 1.. {
            let bb: BigUint = 2u32.into();
            if &n1 % bb.pow(i) != BigUint::zero() {
                s = i-1;
                break;
            }
        }
        s
    };
    let d = &n1/2u32.pow(s);
    let mut vv = [0u8; N];
    rng.fill_bytes(&mut vv);
    let mut a = BigUint::from_bytes_be(&vv);
    while &a >= &n1 {
        rng.fill_bytes(&mut vv);
        a = BigUint::from_bytes_be(&vv);
    }
    for _ in 0..r {
        let mut y = a.modpow(&d,&n);
        if y == BigUint::one() || &y == &n1 {
            continue;
        } else {
            for _ in 1..s {
                y = &y*&y;
                if &y == &n1 {
                    continue;
                } else {
                    return false
                }
            }
        }
        rng.fill_bytes(&mut vv);
        a = BigUint::from_bytes_be(&vv);
        while &a >= &n1 {
            rng.fill_bytes(&mut vv);
            a = BigUint::from_bytes_be(&vv);
        }
    }
    true
}

#[test]
fn test() {
    for _ in 0..5 {
        println!("{}",gen_prime::<64>(3))
    }
}

テストの結果:

running 1 test
8384442391023910089559381911244353428872487399318179273396015886507100438953585665678365906489416840722213922620210237392111730172861086168993311923836111
7117174289562103609173135325238483107034344612555783222340008053564899447108584024808296706859414511376712644654207288058426605121566231239674354322869571
7563023309084596682599574034479919223499374791081086371863943325349323174848426538472921938956966836116894178088405681548677619092250925604872140829006499
6845191684205106300741531169084809632281172184322263496188751375542125818145681444532252331962842985727123170223256606669735938766457801060667384315729663
11140042094397453071834468857251093913894980235722953128972718807724718935294988778661504441980549669968722894429691286957748143827788668127013507759578059
test prime::test ... ok

軽い説明:
まず、1,000,000未満の素数を生成する。続いて、Nバイトの配列をランダムに生成する（ThreadRngはCryptoRngを実装している）。このとき、最上位・最下位のビットは1とする。続いて、その配列を多倍長符号なし整数に変化する。あとは先程生成した1,000,000未満の素数で割り切れないかつMiller–Rabinの素数判定法でacceptされたら終了、そうでなければNバイトの配列をランダムに生成してやり直す。

このプログラムの素数生成ではp-1法などで簡単に因数分解される場合を考慮していないので、間違ってもRSA等の素数生成でこのプログラムを使ってはならない（そもそもまともにデバッグをしていない）