Rustで素数を手軽に求める(primalライブラリ)

この記事は24 days of Rustの内容をなぞったものです。原文はMIT license.

Day 2 - primal | 24 days of Rust

primal素数判定、素数列挙、素因数分解、そしてn番目の素数の上界・下界推定を行うことができるライブラリである。アルゴリズムはエラトステネスのふるいだが、色々と最適化しているらしい。作者のパソコン(i7-3517U)で10^10未満の素数を数え上げるのに2.8秒程度だったそうな。

使い方

1000番目の素数を求めたいとする。

まずふるいを生成するのだが、このとき上界を決めてやる必要がある。適当に大きい数を与えてやるのがある意味楽だが、ここでは同crateのestimate_nth_primeが使える。これはn番目の素数の上界と下界をいい精度で予測してくれるらしい*1。使ってみよう。

extern crate primal;

fn main() {
    let n = 1000;
    let (lower, upper) = primal::estimate_nth_prime(n as u64);
    // => (7816, 8344)
}

上界からふるいを生成する。Sieve::new(upper)してやればよい。

    // fn main に以下を追加
    let sieve = Sieve::new(upper as usize);

1000番目の素数を求めるには以下のように書く。

    let p = sieve.primes_from(0).nth(n - 1).unwrap();
    println!("estimated lower bound: {}", lower);
    println!("estimated upper bound: {}", upper);
    println!("actual {}th prime: {}", n, p);
    // => estimated lower bound: 7816
    // => estimated upper bound: 8344
    // => actual 1000th prime: 7919

primes_from()素数イテレーターを返すメソッドである。イテレーターに対してRustは有用なメソッドを使うことができnthもその一つだ。nthは最初のn-1個のイテレーターを無視し、n番目の値を返す。これは0オリジンなのでn番目の要素を得るにはn-1を指定している。

素因数分解

factorメソッドにより素因数分解を行える。

    let x = 10u64.pow(8) + 1;
    // primal::estimate_nth_prime(x)による推定値
    let sieve = Sieve::new( 2038575494 ); 
    let factors = sieve.factor(x as usize);
    println!("{:?}", factors);
    // => Ok([(17, 1), (5882353, 1)])

返り値がResult型になっていることに注意してもらいたい。factor()の型は以下。

type Factors = Vec<(usize, usize)>;
fn factor(&self, n: usize) -> Result<Factors, (usize, Factors)>

因数分解がうまく行った場合にOkにくるまれたFactorsが返ってくる(これは素因数とその現れる回数のペア、の配列Vecになっている)。

Errはふるいの大きさが足りない場合に返ってくる。(usize, Factors)のうち、usizeは分解しきれなかった値、Factorsはそれまでに素因数分解した結果が入っている。

ちなみに10^8 + 1程度の数を素因数分解しようとすると、デバッグモードでコンパイルしたかリリースモードでしたかによって差が出てくる。timeコマンドでかかった時間を測ってみるといいかもしれない*2

*1:これがなぜうまく動くのかは、私にはよくわからない。ドキュメントを見ると論文へのリンクがある

*2:この記事はWindowsでやっているが、elapsedtimeコマンドというのを知った。僕の結果を見せるとそんなにいいマシン使ってないのがバレるので載せません(謎の見栄)。