並列処理を使うと速くなることがある。
具体的にはSSEを使うと速くなることがある。
http://www.intel.co.jp/jp/download/index.htm
SSEというのは、floatを4つ並べて一つの命令でまとめて計算する技術である。例えば、
x[0] /= y[0];
x[1] /= y[1];
x[2] /= y[2];
x[3] /= y[3];
なんて計算が一発でできるということである。
例として、オイラーによる公式、
を用いてπの計算をしてみよう。
ふつうに書くとこんな感じだろうか。
const int N = 100000000; // 1e8
double sum = 0;
for(int i = N; i > 0; i--)
sum += 1 / ((double)i * i);
fprintf(stdout, "%.8f\n", sqrt(sum * 6));
これで、1.6秒。
どうして2段のループにしなければならないのかよくわからないが、インラインアセンブラでSSEを使うと次のようになる。
float a_sum[4] = { 0, 0, 0, 0 };
float x[4];
float y[4] = { -4, -4, -4, -4 };
float z[4] = { 1, 1, 1, 1 };
for(int i = 10; i > 0; i--) {
int m = N / 40;
x[0] = m * i * 4; x[1] = x[0] - 1;
x[2] = x[0] - 2; x[3] = x[0] - 3;
__asm {
mov eax, m ; eax = m
movups xmm1, [x] ; xmm1[k] = x[k]
movups xmm2, [y] ; xmm2[k] = y[k]
movups xmm0, [a_sum] ; xmm0[k] = a_sum[k]
loop_for:
movups xmm3, [z] ; xmm3[k] = z[k]
movups xmm4, xmm1 ; xmm4[k] = xmm1[k]
mulps xmm4, xmm4 ; xmm4[i] *= xmm4[k]
divps xmm3, xmm4 ; xmm3[k] /= xmm4[k]
addps xmm0, xmm3 ; xmm0[k] += xmm3[k]
addps xmm1, xmm2 ; xmm1[k] -= 4
movups a_sum, xmm0 ; a_sum[k] = xmm0[k]
dec eax ; eax--
jnz loop_for ; if(eax != 0) goto loop_for
}
}
float sum = 0;
for(int k = 0; k < 4; k++)
sum += a_sum[k];
fprintf(stdout, "%.8f\n", sqrt(sum * 6));
__asmの中がアセンブラ。;のあとはコメントでC++ならこういう風に書くところというところ。kは0〜3をとる。とにかく大きさ4のfloatの配列に値を入れて演算する。
これで0.4秒と目論見どおり1/4になった。
