Project Euler 問題10（素数の和）を解いてみる

Project Euler で紹介されている問題10を解いてみました。素直なプログラム例と、エラトステネスの篩を使ったプログラムを試して、アルゴリズムにより計算時間の差を紹介いたします。

問題10（素数の和）

解答例

問題5で使った関数 is_prime を再利用して、200万以下の自然数について、素数か判定して素数であれば加えていきます。

求める数は long int の範囲を超えるため、unsigned long long int を使います。以下は、C の実装です。

#include <stdio.h>

typedef unsigned long long int ull;

#define N 2000000

int is_prime(int n)
{
	int answer = 1;
	int i;

	if (n == 1) {
		answer = 0; /* 1 is not a prime number. */
	}
	for (i = 2; i * i <= n; ++i) {
		if (n % i == 0) {
			answer = 0;
			break;
		}
	}

	return answer;
}

int main(void)
{
	ull answer = 0;
	int i;

	for (i = 2; i <= N; ++i)  {
		if (is_prime(i) == 1) {
			answer += (ull)i;
		}
	}

	printf("Problem010: %lld\n", answer);

	return 0;
}

再利用している関数（is_prime）も含めてソースコードを載せているため、長いように見えますが、main 本体の実装は簡易で素朴な実装です。

考察

驚きました。素朴な実装は、数秒程度の時間がかかると思っていました、かなり速くプログラムが終了しました。cygwin コンソールの time コマンドを使った計測で user time 0.328 秒でした。Project Euler が正解者向けに公開している解説では、500万までの素数の和を求める場合で数秒との記述があり、500万の場合も確かめましが、user time 1.218 秒でした。
paiza.IO で試したところ、200万でも500万でも動作しましたが、1000万で Timeout しました（2022年8月7日の結果です）。

2012年から2013年にかけて Project Euler を楽しんでいたのですが、そのときは、数秒かかっていたと記憶しています。10年で個人用に使うPCもだいぶ速くなったようです。

このような状況ですが、この問題は、エラトステネスの篩（ふるい）が使える典型的な問題ですので、エラトステネスの篩を紹介して、その実装を示します。

エラトステネスの篩の考え方

以下の手順で素数表を作ります。

1. 素数表すべてを値Trueで埋める。
2. もしその素数表に0と1があれば、その値をFalseにする。
3. 素数表の先頭から見ていき、Trueがあれば、それは素数になる（p と置く）。
   （補足　最初の素数は2となるはずです。）
4. p の2倍以上の倍数の値をFlaseにする。
5. 操作3と捜査4を p が素数表の大きさの平方根になるまで繰り返す。
6. 素数表の値がTrueになっているところが素数になる。

古代ギリシャのエラトステネスが考案したとされており、この名前がついています。素数表のために確保するメモリが必要になりますが、効率的なアルゴリズムです。

エラトステネスの篩の動作が分かるアニメ

このアルゴリズムを分かりやすく表現したGIFアニメーションがありましたので紹介します。120までの素数を求めています。$\sqrt{120}=10.954…$ であるため、これ以下の素数である 2、3、5、7 の場合の操作を行えば、残りが素数となります。

この画像は、ドイツ語版ウィキペディアで公開されています。クリエィティブ・コモンズ  表示-継承 3.0 非移植 ライセンスのもとに利用を許諾されています。

プログラム実装

素数表は奇数分だけ用意する実装も可能ですが、実装が簡素となるように自然数をそのまま添え字に使える配列にしました（is_prime_table）。関数 generate_is_prime_table でエラトステネスの篩を用いて素数表を生成します。マクロ関数 is_prime （実体は配列の参照）を用意して、プログラムとして同じように使えるようにしました。
そのため、main 関数は元のものと generate_is_prime_table 呼び出ししか差分がありません（背景色を変更しました）。

#include <stdio.h>

typedef unsigned long long int ull;

#define N 2000000
unsigned char is_prime_table[N+1];

#define is_prime(n) (is_prime_table[(n)])

void generate_is_prime_table(int n)
{
	int i, j;

	for (i = 2; i <= n; ++i) {
		is_prime_table[i] = 1;
	}

	for (i = 2; i * i <= n; ++i) {
		if (is_prime_table[i] == 1) {
			for (j = i + i; j <= n; j += i) {
				is_prime_table[j] = 0;
			}
		}
	}

	return;
}

int main(void)
{
	ull answer = 0;
	int i;

	generate_is_prime_table(N);

	for (i = 2; i <= N; ++i)  {
		if (is_prime(i) == 1) {
			answer += (ull)i;
		}
	}

	printf("Problem010: %lld\n", answer);

	return 0;
}

測定結果

以下は、私の個人PC（プロセッサ：Intel Core i7-10700 @2.90GHz、メモリ 32GB）での測定結果です。キャッシュや他のプログラムの影響もあるため、毎回同じとはなりませんが、速度の目安になるかと思います。測定環境は、cygwin コンソールの time コマンドで測定した user time です。

エラトステネスの篩は、素数表を確保するためメモリは必要となりますが、文字通り桁違いの速さであることが分かります（500万の場合が1千万の場合より遅くなっているのは、実行環境のばらつきだと思われます）。

最後に

問題7の最後に述べた、エラトステネスの篩をここで紹介できました。素朴な方法でも解答は求めることはできますが、効率的なアルゴリズムを用いることにより、劇的に速くなる事例の一つでしょうか。

Project Euler の問題を10問紹介しました。ここまでを振り返る記事を書くことにします。