並列処理でzeta(2)を計算

今回は次の有名な等式(の右辺)をプログラム(C++)で並列処理をして近似計算してみようというのがテーマです。

単純に足し上げる場合と並列処理をした場合で時間がどれくらい異なるかを実験してみました。
ちなみに正確な値は

となっています。
並列化に用いる手法はいろいろとあるようですが、ここでは std::thread を用いるものを考えます。

単純に足し上げるケース

最初に並列処理せずに第1項から順番に1億項まで足し上げるプログラムについて考えます。

#include <thread>
#include <iostream>
#include <vector>
using namespace std;
double zeta2(double n) {
	double sum = 0.0;
	for (long long i = 1; i < n; i++)
	{
		sum +=1.0 / (i*i);
	}
	return sum;
}
int main()
{
	auto start = clock();
	long long N = 100000000;
	cout.precision(14);
	cout << "zeta(2):" << zeta2(N) << endl;
	auto end = clock();
	cout << "time:" << (end - start) / (double)CLOCKS_PER_SEC << endl;
}

実行結果.

zeta(2):1.6449340578346
time:0.623

計測時間は環境によりますが、実行結果を見ると1億項の計算に0.6秒程度かかるようです。

並列処理 (２スレッド)

上の計算を並列化します。一番最初に考えられるのは、偶数項と奇数項に分ける方法です。つまり

として右辺を二つのスレッドに分けて計算します。

#include <thread>
#include <iostream>
#include <vector>
using namespace std;
double zeta2_odd(int n) {
	double sum = 0.0;
	for (int i = 1; 2 * i - 1 < n; i++)
	{
		double d = 2 * i - 1;
		sum += 1.0 / (d * d);
	}
	return sum;
}
double zeta2_even(int n) {
	double sum = 0.0;
	for (int i = 1; 2 * i < n; i++)
	{
		double d = 2 * i;
		sum += 1.0 / (d * d);
	}
	return sum;
}
int main()
{
	auto start = clock();
	long long N = 100000000;
	double sum1 = 0.0;
	double sum2 = 0.0;
	thread th1([&sum1, N] {sum1 += zeta2_odd(N); });
	thread th2([&sum2, N] {sum2 += zeta2_even(N); });
	th1.join();
	th2.join();
	double sum = sum1 + sum2;
	cout.precision(14);
	cout << "zeta(2):" << sum << endl;
	auto end = clock();
	cout << "time:" << (end - start) / (double)CLOCKS_PER_SEC << endl;
}

実行結果.

zeta(2):1.6449340573291
time:0.374

計測時間は 0.623 から 0.374 に改善されました。
ただし、この結果は8コアのパソコンで計測しており、おそらく1コアだとあまり速度は改善されない？のかもしれません。

複数スレッドで計算

２スレッドで計算すると計算速度が向上したので、少し大胆に10スレッドに増やしてみます。計算の分割方法ですが、分母をスレッド数で割ったあまりで分割して計算します。

#include <thread>
#include <iostream>
#include <vector>
using namespace std;
double zeta2(double n, double m, double k) {
	double sum = 0.0;
	for (long long i = 1; m * i - k < n; i++)
	{
		double d = m * i - k;
		sum = sum + 1.0 / (d * d);
	}
	return sum;
}
int main()
{
	auto start = clock();
	long long N = 100000000;
	const int num_thread = 10;
	vector<double> sumV(num_thread);
	double sum = 0.0;
	thread th[num_thread];
	for (int i = 0; i < num_thread; i++)
	{
		th[i] = thread([N, &sumV, i, num_thread] {sumV[i] += zeta2(N, num_thread, i); });
	}
	for (int i = 0; i < num_thread; i++)
	{
		th[i].join();
	}
	for (int i = 0; i < sumV.size(); i++)
	{
		sum += sumV[i];
	}
	cout.precision(14);
	cout << "zeta(2):" << sum << endl;
	auto end = clock();
	cout << "time:" << (end - start) / (double)CLOCKS_PER_SEC << endl;
}

実行結果.

zeta(2):1.6449340569441
time:0.181

最初の計測時間と比較すると、0.623から0.181となり3倍以上速くなったことがわかります。
10スレッド用意したので10倍とはならなかったですが、3倍以上速くなれば充分嬉しいですね。
今回はCPUで並列処理をしましたが、GPUを用いるとどれくらい速くなるのかについても気になるところです。

追記（OpenMPを用いたケース)

上記の場合だとOpenMPを使ったほうがずっと簡潔に書けます。

#include <iostream>
#include <numeric> 
#include <vector>
#include <omp.h>
#include <time.h>
using namespace std;

double zeta2(double n, double m, double k) {
	double sum = 0.0;
	for (long long i = 1; m * i - k < n; i++)
	{
		double d = m * i - k;
		sum = sum + 1.0 / (d * d);
	}
	return sum;
}
int main()
{
	auto start = omp_get_wtime(); 
	int N = 100000000;
	double sum = 0;
	const int num_threads = 10; // スレッドの数.
	vector<double> v(num_threads);
	omp_set_num_threads(num_threads); 
	#pragma omp parallel for 
	for (int i = 0; i < num_threads; i++)
	{
		v[i] = zeta2(N, num_threads, i);
	}
	sum = accumulate(v.begin(), v.end(), 0.0); // 和の計算.
	cout.precision(15);
	cout << sum << endl;
	auto end = omp_get_wtime();
	cout << "time:" << end - start << endl; 
}

実行結果.

1.6449340569441
time:0.18281970010139