この記事はD言語 Advent Calendar 2023の6日目です。

私のいる某分散SNSでいろんな言語で素数を列挙するのが流行ってるようなので
この流れに乗ってD言語で素数を求めるプログラムをいくつか書いてみます。
速度は後半にまとめて記載しています。

免責事項

現状有姿で提供され無保証です。 （PCが吹き飛んでも弁償できません。）
アルゴリズムは得意じゃないので速度は期待しないでください。
深夜テンションで書いたので雑実装です。
配列じゃなく、単に標準出力に出せればOKとします。

ナイーブな実装

prime1.d

import std;
void main(string[] args){
    auto n = args.length >= 2 ? args[1].to!size_t : 100;
    foreach(i; iota(2, n+1)){
        if(!iota(2, i-1).filter!(a=>i%a==0).array.length){
            writeln(i);
        }
    }
}

2からnまでの各数について、
2から自身-1の値で割ったときに割り切れる数が出ない値を表示します。
Dっぽく、iotaとfilterでレンジを活用していますがメモリ効率は極悪です。

$ dmd prime1.d
$ ./prime1 10
2
3
5
7

雑なエラトステネスの篩版

prime2.d

import std;
void main(string[] args){
    auto n = args.length >= 2 ? args[1].to!size_t : 100;
    auto r = n^^0.5;
    T[] f(T)(T[] x, T p){
        writeln(p);
        return (p>r)?x[1..$]:f(x.filter!(a=>(a%p!=0)).array, x[1]);
    }
    f(iota(2, n+1).array, 2).each!writeln;
}

大分速くなりましたが、再帰使ってるので篩の性能が活かせていません。

で、ここまでやって思い出したのですが・・・

エラトステネスの篩版（公式版）

はい、公式のサンプルにありました。
記事書き始めてから思い出したので仕方ないね。

prime3.d

/* Eratosthenes Sieve prime number calculation. */

import std.conv;
import std.stdio;
import std.range;

void main(string[] args)
{
    immutable max = (1 < args.length)
        ? args[1].to!size_t
        : 0x4000;
//    size_t count = 1; // we have 2.
    writeln(2);
    // flags[i] = isPrime(2 * i + 3)
    auto flags = new bool[(max - 1) / 2];
    flags[] = true;

    foreach (i; 0..flags.length)
    {
        if (!flags[i])
            continue;
        auto prime = i + i + 3;
        foreach (k; iota(i + prime, flags.length, prime))
            flags[k] = false;

        writeln(prime);
//        count++;
    }
//    writefln("%d primes", count);
}

けっこう昔からあったようです。
オリジナルの実装だと、素数の数をカウントするだけなのでそこだけ修正しています。

速度比較

$ time ./prime1 10000 | wc
   1229    1229    5948

real    0m0.268s
user    0m0.267s
sys     0m0.000s
$ time ./prime1 100000 | wc     $ time ./prime2 100000 | wc
   9592    9592   56126            9592    9592   56126

real    0m26.107s               real    0m0.027s
user    0m26.094s               user    0m0.029s
sys     0m0.000s                sys     0m0.000s
                                $ time ./prime2 1000000 | wc
                                  78498   78498  538468

                                real    0m0.477s
                                user    0m0.473s
                                sys     0m0.007s
                                $ time ./prime2 10000000 | wc   $ time ./prime3 10000000 | wc
                                 664579  664579 5227116          664579  664579 5227116

                                real    0m9.998s                real    0m0.137s
                                user    0m9.740s                user    0m0.145s
                                sys     0m0.253s                sys     0m0.004s
                                                                $ time ./prime3 100000000 | wc
                                                                5761455 5761455 51099000

                                                                real    0m1.714s
                                                                user    0m1.787s
                                                                sys     0m0.033s
                                                                $ time ./prime3 1000000000 | wc
                                                                50847534 50847534 501959790

                                                                real    0m18.823s
                                                                user    0m19.526s
                                                                sys     0m0.305s
                                                                $ time ./prime3 10000000000 | wc
                                                                455052511 455052511 4948214537

                                                                real    3m16.855s
                                                                user    3m24.470s
                                                                sys     0m2.815s

裏作業しながら測ってるのであくまで参考までに。

やはり変に横着してやるより公式のコードは早いですね。
prime3の100億で5GBくらいメモリ食ってました。

おまけ

D言語使いたるもの、コンパイル時に素数を出したいですね。（これがやりたかっただけ）

primect.d

import std;
void main(){
    immutable n = 10;
    foreach(i; aliasSeqOf!(iota(2, n+1))){
        static if(!iota(2, i-1).filter!(a=>i%a==0).array.length){
            pragma(msg, i);
        }
    }
}

$ dmd -c -o- primect.d 2>&1
2
3
5
7

さすがD、実行する必要すらないですね。
ナイーブ実装はこれを書くために用意したものだったり。（篩版は諦めた）

ちなみに迂闊に n = 10000; とかやってしまうとメモリが吹き飛んで死にます。
OOM-KillerさんがいなければOSごと即死だった・・・

感想

コンパイル時D言語くんの記事を書いたのが7年前だということに戦慄しています。

なお今回の最大の収穫は
今までずっと エラ「スト」テネス だと思ってたのですが、
正しくは エラ「トス」テネス だったと気づけたことでしたorz