１．はじめに
Ｘ(旧Twitter)で「素数生成プログラミングコンテスト」を開催すると言う下記のメッセージを見かけ、面白そうだったので参加してみることにしました。とは言っても事前に参加申請するわけではなく、いい案が思い付いたらＸでリプライするという運用なので気軽です。
以前も Twitter でエラトステネスの篩を使ってペラン擬素数を求めるプログラムの高速化について議論したことがあり、Z80の最適化と同様にこのような話題は好物なのです。
２．エラトステネスの篩
素数の列挙と言えば「エラトステネスの篩」ですが、パソコン関連の雑誌でもたまに見かけることがあります。国会図書館のサイトで検索してみたらコンピュータサイエンス誌のｂｉｔに何回か掲載されたことがあったようです。
ネット上にも沢山の情報がありますが、「いかたこのたこつぼ」というサイトの「
素数判定・列挙」のページに良くまとめられた情報がありました。
３．素数探索プログラム
はじめに素数を列挙していくプログラムを考えてみました。この時、素数か否かのチェックロジックは検討対象外とし、暫定処理として「それまでに抽出した素数で割り切れなかったら素数と判断する」と言うロジックにしました（後で判ったことですが素数判定ではこの方法が最強の手法の一つのようです）。
素数探索処理を高速化するためにはチェック対象を如何に絞っていくか肝になります。
この時の思考過程の概要を要約すると
- 最初に思い付くのは偶数のチェック対象から外すことです。
- 次に2と3の倍数をチェック対象から除外すればもっと速くなります。
- 5迄の倍数を除くには30の倍数＋1,7,11,13,17,19,23,29の確認で済みます。
- 30の値が出てくるのはそれまでの素数の最小公倍数（=素数を掛け合わせた値）が30だからです。即ち素数の倍数の情報は最小公倍数の周期性があります。
- であれば5よりもっと大きい素数までの倍数を列挙した配列を作り、その配列情報は最小公倍数の周期で無限大まで使えるはず。
- でも倍数情報をプログラムに埋め込むのは大変 ⇒ エラトステネスの篩を使えば自動生成できる。
ということで応募したものが下記になります。まぁ、ソース中の Until を N に変更すれば「エラトステネスの篩」そのものになるのですが、今回の手法はメモリの制限ぎりぎりまでエラトステネスの篩を使い、更にその先の素数を求めていく場合などには有効な手法と言えるかもしれません。
https://paiza.io/projects/35yPplZc8zelsL_aR7QbDg
４．「エラトステネスの篩」の高速化
次にエラトステネスの篩の高速化について考えてみました。エラトステネスの篩は良くできていて検出済みの素数の最大値（Ｎ)迄の倍数のフラグを更新（素数ではないというフラグに変更する）し、フラグ更新されなかった数を素数と判定するという実に直球的なロジックです。
少し考えると気が付きますが、ｋの倍数のフラグを更新する場合は
- 更新は k*k から始めればいい
- 偶数倍は素数では無いので倍数としては k*k に 2*k づつ足していった数のフラグをＮ以下の間更新すればいい
- 倍数のフラグ更新処理は k がルートＮ以下の時だけでよい
の条件で行えばいいことになります。
フラグ更新と言っても対応するフラグの配列を一定の値（例えば１）に変更するだけなので、高速化のために既に更新済みならスキップする等の条件を追加すると寧ろ遅くなるケースが殆どです。
そこで次の気付きを基に高速化してみました。
- 前述しましたが倍数でのフラグの更新状態はそれまで検出した素数の最小公倍数の周期で繰り返され、かつ素数が小さい程、更新密度が高い。そこで始めのうちは１周期分のフラグのみ更新し、後から後方に一括コピーする
- フラグ更新時に更新済みの3の倍数のスキップ処理をレジスタのみで処理可能な軽量な手法で実装
- ループ内の条件文を削減するためにループを3つのフェーズに分解
これらを実装したものが
https://paiza.io/projects/Ru-Azjr-DH13rFtulZWjFg です。ソースも貼っておきます。
その後、上記の篩の１周期は中心点に対して点対称であることに気が付き、半周期だけ求めた後にコピーするようにしてみましたが高速化への効果は殆ど認められませんでした。
|
エラトステネスの篩の高速化（Ｃ＋＋）
|
#include <chrono>
#include <stdio.h>
unsigned long *P, N = 10000000;
unsigned long get_prime(){
char *era = new char[N/2+1](); // 初期化有 篩は奇数部のみ
using prime_t = unsigned long;
using commul_t = struct { prime_t until; prime_t commul; };
commul_t com[] = { // 素数とそれらの最小公倍数のデータ
// { 7>>1, 210>>1 }, // 細かく刻んでも効果が出ないので一部コメントアウト
// { 11>>1, 2310>>1 },
{ 13>>1, 30030>>1 },
{ 17>>1, 510510>>1 },
{ 0>>1, N>>1 }
};
commul_t *pcom = com;
prime_t i, prime,pidx=0;
prime_t sqrtn=0; // Nの平方根（sqrtn）を求める
for ( prime_t n=N,tmp,bit=(prime_t)1<<(sizeof(prime_t)*8-2); bit; bit >>= 2 ) {
sqrtn >>= 1;
if ( n >= (tmp = (sqrtn << 1) | bit ) ) {
n -= tmp;
sqrtn |= bit;
}
}
prime_t N2=N>>1, sqrtn2=sqrtn>>1, until=com[sizeof(com)/sizeof(commul_t)-2].until;
char *puntil = &era[pcom->commul + pcom->until];
for ( i=1; i<=until; i++ ) {
if ( era[i] == 0 ) {
P[++pidx] = prime = (i<<1)+1;
prime_t mod3 = prime % 3, add3 = prime<<1; // 3の倍数のフラグを書換しないことで高速化
for ( char *pera=&era[(prime*prime)>>1]; pera<=puntil; ) {
*pera = 1;
if ( --mod3 == 0 ) { // prime=3の場合のみこの３の倍数スキップは機能しない
mod3 = 2;
pera += add3;
} else {
pera += prime;
}
}
if ( i == pcom->until ) { // フラグの後方コピー処理
char *spnt=&era[i+1],*dpnt=spnt+pcom->commul;
pcom++;
puntil = &era[pcom->commul + pcom->until];
while ( dpnt!=puntil ) {
*dpnt++ = *spnt++;
}
}
}
}
for ( ; i<=sqrtn2; i++ ) { // フラグの後方コピーが不要なループ
if ( era[i] == 0 ) {
P[++pidx] = prime = (i<<1)+1;
prime_t mod3 = prime % 3, add3 = prime<<1; // 3の倍数のフラグを書換しないことで高速化
for ( char *pera=&era[(prime*prime)>>1]; pera<=puntil; ) {
*pera = 1;
if ( --mod3 == 0 ) {
mod3 = 2;
pera += add3;
} else {
pera += prime;
}
}
}
}
for ( prime_t i=sqrtn2+1; i<N2; i++ ) { // フラグ更新が不要なループ
if ( era[i] == 0 ) {
P[++pidx] = (i<<1)+1;
}
}
delete[] era;
return pidx;
}
int main(){
unsigned long bound;
P = new unsigned long[999999]; P[0] = 2;
auto start = std::chrono::system_clock::now();
bound = get_prime();
auto end = std::chrono::system_clock::now();
printf("%lu （%lu番目の素数）です¥n", P[bound], bound + 1);
printf("Time: %u(ms)¥n", std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count());
delete[] P;
return 0;
}
|
|
５．「エラトステネスの篩」の亜種
エラトステネスの篩の亜種として
- 検出済みのそれぞれの素数の倍数をチェック対象の値以上のもの一つだけを昇順にリストで管理（倍数リスト）する
- チェック対象の値が倍数リストの先頭値（＝最小値）と同じであればチェック対象は合成数なので先頭の倍数を次の倍数にアップデート
- チェック対象が倍数リストの先頭と違う値の時は素数と判定し、倍数を倍数リストに追加
- チェック対象数に2を加え次のチェック対象としＮになるまで(2)からの処理を繰り返す
というものを考えてみました。
https://paiza.io/projects/tLZU5FwgeSiGpPbN9MpeoA
素数判定と抽出した素数の倍数を倍数リストへ追加する処理は軽いのですが倍数のアップデート処理が重く、処理時間が 800ms 程で上記の高速化したもの（15ms）よりかなり遅くなってしまいました。
リストも貼っておきます。
|
エラトステネスの篩の亜種（Ｃ＋＋）
|
#include <chrono>
#include <stdio.h>
unsigned long *P, N = 10000000;
unsigned long get_prime(){
using prime_t = unsigned long;
using primul_t = struct PRIMUL { PRIMUL *next; prime_t mul; prime_t add; };
primul_t *primul = new primul_t[ 999999 ]();
prime_t pidx=0;
primul_t *pp,*ppadd,*root,*tail,*last,**prev;
primul[0] = {nullptr, N+1, 0}; // 所謂番兵
tail = &primul[0];
P[++pidx] = 3;
primul[pidx] = {tail, 3*3, 3 << 1};
last = &primul[pidx];
prime_t sqrtn=0; // Nの平方根（sqrtn）を求める
for ( prime_t n=N,tmp,bit=(prime_t)1<<(sizeof(prime_t)*8-2); bit; bit >>= 2 ) {
sqrtn >>= 1;
if ( n >= (tmp = (sqrtn << 1) | bit ) ) {
n -= tmp;
sqrtn |= bit;
}
}
root = &primul[pidx];
for ( prime_t val=5; val<=N; val+=2 ) {
pp = root;
prev = &root;
if ( val == pp->mul ) { // リスト内の倍数の最小値と一致すれば合成数
prime_t mul;
mul = pp->mul += pp->add; // 倍数をアップデート
*prev = pp->next; // ppノードをリストから外す
ppadd = pp;
pp = pp->next;
while ( mul <= N ) {
while ( mul > pp->mul ) {
prev = &pp->next;
pp = pp->next;
}
if ( mul == pp->mul ) {
mul += ppadd->add; // 被っていれば次の倍数
ppadd->mul = mul;
} else {
ppadd->next = *prev; // ppの前に追加
*prev = ppadd;
if ( *prev == tail ) {
last = ppadd;
}
break;
}
}
} else { // 倍数の最小値と一致していなければ素数
P[++pidx] = val;
if ( val <= sqrtn ) {
ppadd = &primul[pidx];
*ppadd = { tail, val*val, val << 1};
last->next = ppadd; // リストの最後に追加
last = ppadd;
}
}
}
delete[] primul;
return pidx;
}
int main(){
unsigned long bound;
P = new unsigned long[999999]; P[0] = 2;
auto start = std::chrono::system_clock::now();
bound = get_prime();
auto end = std::chrono::system_clock::now();
printf("%lu （%lu番目の素数）です¥n", P[bound], bound + 1);
printf("Time: %u(ms)¥n", std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count());
delete[] P;
return 0;
}
|
|
６．その他
思い付きでパイプライン処理の効率低下の原因となる分岐命令の追加無しでスキップ処理を実装（結構めんどい）してみましたが高速化の効果は認められませんでした。
また、思考過程で５以上の素数の二乗に２を加えると３の倍数になることに気付き、ＡＩに確認してみたところそんなことは無いと言われググっても出てこないので新発見か？
と一瞬思いましたが冷静に考えたら直ぐに証明できたので文書化する程の事でも無い為、検索でヒットしなかったようですｗ
posted by skyriver at 00:33|
Comment(0)
|
その他
|
|
)
){kind=link}