Ｎが合成数の時、Ｎは√Ｎ以下の約数を必ず持つことの証明

C - エラトステネスの篩

// エラトステネスの篩

#include <stdio.h>
#include <math.h>

#define NUM 50  // NUM以下の素数を求める

int main(void)
{
  unsigned i, j;
  unsigned sq_num = (int)sqrt((double)NUM);
  unsigned prime[NUM];

  // 1が立っているものが素数
  // 0が立っているものが合成数（素数ではない）

  for (i = 0; i < NUM; i++)
    prime[i] = 1;  // 全ての数に1を立てる
  prime[0] = 0;    // 1は素数ではない

  for (i = 1; i < sq_num; i++) {
    if (prime[i] == 1)          // prime[i]が素数なら
      for (j = (i+1); (i+1) * j <= NUM; j++)
        prime[(i+1) * j - 1] = 0;  // 素数の倍数は素数ではない
  }

  for (i = 0; i < NUM; i++)
    if (prime[i] == 1)  // 1が立っているものが素数
      printf("%3d", i + 1);
  putchar('¥n');

  return (0);
}

訂正 2010/10/14 ： for (i = 1; i+1 < sq_num; i++)⇒for (i = 1; i < sq_num; i++)

訂正 2011/05/17 ： for (j = 2; (i+1) * j <= NUM; j++)⇒for (j = (i+1); (i+1) * j <= NUM; j++)

エラトステネスの篩のアルゴリズム解説

エラトステネスの篩はある整数Ｎ以下の素数を見つけるアルゴリズムです。

まずはじめに２からＮまでの整数全てにフラグ１を立てる。１が立っている整数が素数である。次に１は素数ではないので０を立てる。これらが前準備。

次に２から順番に√Ｎまでの整数を見ていく。（上の証明よりNは√N以下の約数を持つ）

整数２には１が立っているので２は素数である。素数の倍数は合成数（素数でない整数）であるので、２（素数）の倍数全てに０を立てる。
２の次の整数３には１が立っているので３は素数である。同様にして、３（素数）の倍数は合成数（素数でない整数）であるので３の倍数全てに０を立てる。３の次の整数４には０が立っているので合成数（素数でない整数）０が立っている数は何もせずに次に進む。
４の次の整数５には１が立っているので５は素数である。同様にして、５（素数）の倍数は合成数であるので５の倍数全てに０を立てる。

このようにして２から√Ｎまでの整数を順番に見ていき整数に１が立っている時だけその整数の倍数に０を立てる作業をする。整数に０が立っているときは何もせずに次の整数に着目する。これを繰り返していき最後に１が立っている整数が素数。

ソースコード解説

上のコードでは配列の添え字と調べる整数値が一致していません。prime[i]はi+1という整数と一致します。もし配列の要素と調べる整数値との値を一致させると、コードが読みやすくなる代わりに配列の要素を一つ無駄にしてしまいます。prime[0]を使わないことになりますから。

もう少し詳しく説明します。例えば、i = 2の時を考えてみましょう。ちなみにi = 2のときに着目する整数は3です。prime[0] : prime[1] : prime[2] = 1 : 2 : 3

for (i = 1; i < sq_num; i++) {
  if (prime[i] == 1)          // prime[i]が素数なら
    for (j = (i+1); (i+1) * j <= NUM; j++)
      prime[(i+1) * j - 1] = 0;  // 素数の倍数は素数ではない
}

i = 2の時prime[2] には1が立っているので（初めに全ての要素を１で初期化しています）if (prime[i] == 1)は真になり中にあるfor文が実行されます。j = (i+1)からなのは3の倍数は素数でないということから3の倍数に対応する要素全てに０を立てるのですが、その際3の０倍と３の1倍は見る必要がないからです。また、(i+1)-1倍も見る必要がありません。この場合ですと、３×２は調べないということです。なぜなら、i = 2の時に２×３を消してあるからです。

(i+1) * j <= NUM;これは今 i は2ですが調べたいのは3の倍数なので i+1 としています。NUM以下の３の倍数に0を立てています。

prime[(i+1) * j - 1] = 0; これは例えば i = 2, j = 3のとき(i+1) * jつまり(2+1)×3 = 9に0を立てたいのですが整数9に対応する添え字はprime[8]なので -1 としています。

JavaScript - エラトステネスの篩シミュレータ : Please Comment on My Codeこちらもどうぞ。

追記 2010/10/14 ： ソースコード解説

追記 2011/05/17 ： 下線部分を加筆。