素数を列挙する3 - dachs

キーボードが小さくなって誤入力倍増キャンペーン開催中．

##単純なエラトステネスの篩

use Benchmark;

$x=primefactory1();
$count=5000;
$i=0;

timethese($count, {'p1' => '$x->();', 
                  }
         );


sub primefactory1{
    my @primes=(2,3);
    my $state=-1;
    return sub{
        $state++;
        return $primes[$state] if $state==0 or $state==1;
        my $prime_cand=$primes[-1]+1;### (1)
        SIEVE:
        foreach my $d (2 .. int(sqrt($prime_cand))+1){
            #$i++; 
            if ($prime_cand % $d ==0){
                $prime_cand++;###(1)
                goto SIEVE; 
            }
        }
        @primes=(@primes,$prime_cand);
        return $primes[-1];
    }
}

これだと結果は

Benchmark: timing 5000 iterations of p1...
        p1:  4 wallclock secs ( 4.08 usr +  0.00 sys =  4.08 CPU) @ 1226.09/s (n=5000)

となります．そして，$iをいかして篩本体のループの回数を数えると1020148回でした．次に，(1)の部分を +2 にしてみます．素数の隣は素数ではないので一個とばしても問題ありません．結果は（見えてるのですが）

Benchmark: timing 5000 iterations of p1...
        p1:  4 wallclock secs ( 4.00 usr +  0.00 sys =  4.00 CPU) @ 1250.00/s (n=5000)

実際は何回もやってるのですが，ほとんど同じです．ループ回数を増やしたり統計処理をすべきかもしれませんが，やっていません．篩のループ回数は995844で多少減ってるのですが，減ってる部分はもともと大きなコストがなかった部分なので影響はそれほどではないということでしょう．

なお，篩の中でのsqrt(x)をx/2に変えるとここの計算量は減りますが，その分ループが大量に増えるので（1015150から57537466）論外です．また，(1)を代入だけにして直後に++する（+1よりも++演算子の方が速い）のも遅くなることも確認しました（代入＋「++」よりも「+1の結果を代入」の方が速いようです）．

となると本質は最初から自明だった気もしますが，篩の本体のループ（の中でコストの大きい部分）をいかに減らすかです．とにかくこれでは平方根以下の自然数を全部チェックする感じなので，そこをなんとかしようとするわけです．実際は「すでに分かっている素数」だけで十分です．

foreachを次のようにしてみます．

        foreach my $d ( @primes ){

結果，まず篩のループ回数は12609711と当然ながら増えて，

Benchmark: timing 5000 iterations of p1...
        p1:  5 wallclock secs ( 5.55 usr +  0.00 sys =  5.55 CPU) @ 901.55/s (n=5000)

時間もかかります．不要な大きい素数まで考慮してるせいもあるはずです．そこで，foreachにかける配列を先にgrepして

 SIEVE:
 my @reduced_primes=grep { $_<= int(sqrt($prime_cand))+1} @primes; 
 foreach my $d ( @reduced_primes ){

ということを考えると，これはgrepのコストが大きすぎて論外です．ループ回数自体は279478と激減ですが，

Benchmark: timing 5000 iterations of p1...
        p1: 32 wallclock secs (31.70 usr +  0.00 sys = 31.70 CPU) @ 157.71/s (n=5000)

という有様．grepしなくても，そもそもforeachで一個一個取り出しているのですから，grepみたいなことをすでにやってます．$dをチェックすることにします．

sub primefactory1{
    my @primes=(2,3);
    my $state=-1;
    return sub{
        $state++;
        return $primes[$state] if $state==0 or $state==1;
        my $prime_cand=$primes[-1]+2;##+1;
        my $th;
        SIEVE:
        $th = int(sqrt($prime_cand))+1;
        foreach my $d ( @primes ){
            last if $d>$th;
            if ( ($prime_cand % $d ==0)  ){
                #$prime_cand++;
                $prime_cand = $prime_cand + 2;
                goto SIEVE; 
            }
        }
        @primes=(@primes,$prime_cand);
        return $primes[-1];
    }
}

結果は，ループ回数284475で少し増えて（ループに入ってすぐ出る分が増えて）

Benchmark: timing 5000 iterations of p1...
        p1:  4 wallclock secs ( 3.92 usr +  0.00 sys =  3.92 CPU) @ 1274.86/s (n=5000)

となりました．微妙ですが，これが一番速いようです．

この「単純なエラトステネスの篩」の方針では，理屈ではこれ以上の効率化はないように思います．となると，Haskellのときのご指摘のような「mod 6で1か5」を探索する手法を考えてみようとなります．