連休の直前まで連休の存在を認識していない人間が計画的な連休を過ごせる筈もなく、部屋でアニメのBDを見たりギターを弾いたりプログラム的なネタで遊んだりして過ごしている。
そんな暇つぶしの一環として記事でも書こうかと久々にblogを見てみたら、前回書いた記事に面白過ぎるコメントがついているが、これはこれで面白すぎるのでそのままにしておく。それはそれとして以前twitterで流れて来たコードで「これアカンやつや」と思ったものについて、何がアカンのか解説する記事でも書いてみようかと思った。
というのも、今仕事を受けている現場のエンジニアに「これアカンですね」と紹介したところイマイチ共感を得られなかった、という事実から、ひょっとしたら世間にはこれの何がダメなのか理解できないエンジニアがゴロゴロしているのではないかと思ったからだ。問題のコードを書かれた方にとってはこんなところでネタにされるのは降って湧いた災難のような話かもしれないが、ダメなコードを見てみぬふり出来ぬのがエンジニアの性である。
で、問題のコードはこちら。
このコードが何をやっているかを端的に言えば、
ある値が素数であるかどうかを判定するために、「2以上その値未満の自然数全て」で「実際に割り切れるかどうか試し」ている。そして、そのテストを探索範囲にある全ての値について行っている。
ダメな理由は次の通り。というか、それなりにまともなプログラマならばコードを見た段階で「これアカンやつや」という感想を直感的に抱けるはず。
判定対象の値nが素数であるならば、(n-2) 回ループを回ることになる。nの値が億単位であれば、億単位の回数ループすることになる。現在のコンピュータでも億単位の回数をループさせるのはほとんどフリーズさせるに近い。
おそらく数学的な素養の無い方が書いたのだろうが、素人考えでもループの回数は約半分に減らせる。まともな解決策であれば、億単位以上の値でもループ回数を数万分の1に減らせる。ひとまず素人考え的に「偶数はテストしない」というところから見ていこう。
素人考え: 2以外の偶数は試す必要はない。
偶数で割り切れる自然数は、必ず2で割り切れる。つまり、「割って試す」値として、2以外の偶数を試す必要はない。これだけで、割って試す値の総数は約半分になる。最初に偶数であるか否かをテストし、奇数である場合のみ、3以上の奇数で試す。まずこれだけで判定対象は半分になる。
これは外側のループでも内側のループでも同様で、2以外の偶数はそもそも素数に登場せず、テストするまでもなく素数ではないため、最初から判定対象から外しておくべきだ。外側のループで偶数を除外すれば、内側のループでは2で割り切れるかどうかのテストも不要となり、内側のループでは3から開始される奇数のみでテストすればよくなる。
これにより、1億近い素数を判定にかけた場合でも5千万回ほどテストを省略できる。残る5千万回の奇数のうちに対象値を割り切れる値があれば、その時点でループは中断される。
しかし、実際に1億近い素数をテストにかけた場合、残る5千万回は全てループすることにはなってしまう。「半分も減った」のは事実だが、「まだ半分も残っている」のも事実だ。
…だが、この素人考えには重要なヒントがある。それを踏まえて「まともな方法」を解説する。
素人考えのその先:「エラトステネスの篩(ふるい)」
先ほどの素人考えで、「偶数で割り切れる自然数は、必ず2で割り切れる」と書いた。これを拡張し、以下のように考えてみる。
これを踏まえて素数の性質を再度確認する。
つまり、このプログラムの命題において、「割って試す値」としては、少なくとも「そこまでに登場した素数」以外を試す必要はない。以下要点。
さらにもう一つ。
つまり、n が素数か否かは√n 以下の素数で割ってみればわかる、ということになる。
件のプログラムは受け取った探索範囲の中から素数を探す、という機能要件が与えられているが、最初に探索範囲の最大値の平方根を上限として素数リストを作ってしまえば用が足りる、ということになる。また、素数リストの各値で割る際、その素数の2乗が判定対象値を超えた値では判定の必要が無い、ということでもある。
これを踏まえると判定を行うループの回数を劇的に減らすことができる。1億前後の値を判定する際も1万以下の素数との間で判定すればよく、1万以下の素数は1229個しかない。
一つの値につき1億回ループすることに比べれば、1/80000以下のループ回数で済んでしまう。
件のプログラムは動作サンプルにもあるとおり、1000未満の値であればまあそこそこ実用に足りるだろうが、億以上の素数でも同様に使えるか? を問うた時点で破綻する。そこが「これアカンやつや」の所以である。
この方法は古代ギリシャ人エラトステネス(B.C.275-B.C.194)が考案したとされる「エラトステネスの篩(ふるい)」と呼ばれる方法に基づいている。
21世紀のエンジニアが紀元前のギリシャ人に頭の出来を鼻で笑われるようでは先人に申し訳が立たない。先に挙げたコードに違和感を感じなかった方はぜひ覚えていただきたい。
そして実際のコード
理論だけ紹介して実装が無いと「口先ばっかりやないか」と言われそうなので、ちゃんとコードも実装してみてある。件のコードの置いてあるサイトに書かれているライセンスを見ると、改変した場合でもオリジナルコードの出所を書け、と書いてあるので、「すみませんが俺はそのコード要りません」という思いを込めて、結局フルスクラッチで書く、という方向で。
探索範囲は scanf() ではなくコマンドライン引数として与えるようになっていて、値をとるあたりの実装はえらく手抜きではあるが、探索範囲に100000000〜100020000 を指定した際の1081個の素数を見つける速度などを比べてみて欲しい。
…多分あちらのコードでは開始するなりフリーズしたかのように見えると思うけど。最初の素数100000007に出くわした時点で1億5回分のループを回る羽目になるので。
#include <stdio.h> #include <stdlib.h> #include <math.h> typedef struct PRIME { struct PRIME * next; int num; } PRIME; typedef struct PRIME_LIST { PRIME * begin; PRIME * end; } PRIME_LIST; static PRIME_LIST prime_list; void init(void) { prime_list.begin = NULL; prime_list.end = NULL; } /* 現在の素数履歴で、与えられた値を割り切れるものがあるかを判定 */ int isPrime(int num) { PRIME * prime; for(prime = prime_list.begin; /* 履歴先頭から判定 */ prime && prime->num * prime->num <= num; /* 末端に達するか、素数の値が num の平方根以上になったら終了 */ prime = prime->next) { /* 次の素数へ */ if((num % prime->num) == 0) return 0; /* 割ってみて余りが 0 ならば割り切れたので素数ではない。 */ } return 1; /* numを割り切れる素数が一つもなければ、numは素数 */ } /* 探索範囲の上限値を与え、 その上限値までの素数を判定するのに必要な素数履歴を作る。 判定に必要な値は必ず上限値の平方根以下の自然数となる。 */ int makePrimeList(int area_max) { int find_max = (int)sqrt((double)area_max); /* 2 以外は奇数のみを判定しリスト化する */ for(int i = 2; i <= find_max; i+=1+(i&1)) { /* その時点までに見つかっている素数の中で、 iを割り切れるものがあるかを判定 */ if(isPrime(i)) { /* 素数なので、素数リストに項目を追加する。*/ PRIME * prime; if(NULL == (prime = (PRIME *)malloc(sizeof(PRIME)))) { /* 素数リストの項目を作るのに失敗したらエラー */ return -1; } prime->next = NULL; prime->num = i; if(prime_list.end) { prime_list.end->next = prime; } else { prime_list.begin = prime; } prime_list.end = prime; } } return 0; } /* うん。まあ main() 関数は切り離せたほうがいいよね。 */ int Eratosthenes(int min, int max) { /* 与えられた探索範囲を判定する */ int prime_count = 0; if(!(min&1) && (min != 2)) min++; for(int n = min; n <= max; n += 1+(n&1)) { /* 2以外の偶数はそもそもテストの必要さえない。*/ if(isPrime(n)) { if(!(prime_count % 10)) printf("\n"); printf("%6d, ", n); prime_count++; } } return prime_count; } int main(int argc, char **argv) { /* 探索範囲をコマンド引数から取得する。わりと手抜き。 */ int min = atoi(argv[1]); int max = atoi(argv[2]); /* 上限と下限が入れ替わっているとアレなので、順番を整える。 */ if(min > max) { /* 論理演算を用いた値の交換。 なぜこれで値が入れ替わるのかわからない人は勉強すること。*/ min ^= max; max ^= min; min ^= max; } /* 素数履歴初期化 */ init(); /* √max 以下の素数履歴作成。*/ if(makePrimeList(max)) { return EXIT_FAILURE; /* 履歴を作れなかった */ } /* エラトステネスの篩により、探索範囲の素数を探す */ int prime_count = Eratosthenes(min, max); printf("\n\n見つけた素数: %d個\n", prime_count); return EXIT_SUCCESS; }
