素数分布（1）
手始めに、ガウスの故事にならって、素数分布を求めよう。
1000ずつ区切って、そのうちの素数の数を数える。
1001〜1000万を1000ずつ区切ることにする。
const int nWidth = 128;
const int nHeight = 128;
const int nData = nWidth * nHeight * 4;
float *data = new float[nData];
const int nIntervals = 10000;
const int nRep = 1000; // must be even
for(int i = 1; i < nIntervals; i++) {
data[i*4] = (float)(i * nRep + 1);
}
として、4バイトの最初に1001、2001、…と入れておく。
データは2のべきにならなければいけないため余る。
そこには0を入れておき、0なら計算しないようにする。
data[0] = 0;
for(int i = nIntervals; i < nData / 4; i++) {
data[i*4] = 0;
}
シェーダプログラムは、
#extension GL_ARB_texture_rectangle : enable
#extension GL_EXT_gpu_shader4 : enableuniform sampler2DRect texUnit0; // data
uniform int nRep;
bool isPrime(in int n) {
for(int p = 3; p * p <= n; p += 2) {
if(n % p == 0)
return false;
}
return true;
}void main(void) {
int nPrimes = 0;
vec4 color = texRECT(texUnit0, gl_FragCoord.xy);
int iStart = int(color.x);
if(iStart != 0) {
for(int k = iStart; k < iStart + nRep; k += 2)
nPrimes += int(isPrime(k));
}
gl_FragColor = vec4(nPrimes, 0, 0, 0);
}
3以上の奇数で割っていき、1000個の自然数のうちの素数の数を数える。
一方、CPU側は、
bool isPrime(int n) {
for(int p = 3; p * p <= n; p += 2) {
if(n % p == 0)
return false;
}
return true;
}...
for(int i = 1; i < nIntervals; i++) {
int counter = 0;
for(int k = i * nRep + 1; k < i * nRep + nRep; k += 2) {
counter += (int)isPrime(k);
}
printf("%7d %3d\n", i * nRep, counter);
}
結果は、
GPU : 2.2sec
CPU : 12.0sec
と思ったように速くならなかった。
GPUは非常に分岐に弱い。
同時に走るスレッドのうち、ある数のスレッドがまとまって実行される。
そのスレッドは同じ命令を実行することになる。
すなわち、分岐して違う命令を実行することになると、
どちらかのスレッドはそこで待っているのと同じことになる。
だから、全てのスレッドがなるべく同じ動作をするように組まなければならない。
そこで、分岐しないように、
あるところで非素数となっても割り算を続けるようにした。
bool isPrime(in int n) {
bool ret = true;
for(int p = 3; p * p <= n; p += 2) {
bool indivisible = n % p != 0;
ret = ret && indivisible;
}
return ret;
}
あるいは、
bool isPrime(in int n) {
int counter = 0;
for(int p = 3; p * p <= n; p += 2) {
counter += int(n % p == 0);
}
return counter == 0;
}
しかし、かえってやや遅くなってしまった。
割り算が遅いのだろうか。
