前にも取り上げたが、コラッツの問題というのは、どの自然数も、奇数なら3倍して1を足す、偶数なら2で割る、という操作を繰り返すといつか1になるか、というものである。例えば、
7→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1
となる。
数値実験では、今のところ1にならない例は見つかっていないらしい。発散しないのは直感的にわかるが、ループをまったくしないというのは驚異的である。なにか未知の構造が整数にはあるのだろうか。いや、ガロア以来の構造がどーとかいう考え方では、この問題には太刀打ちできないのかもしれない。
今回は、GPUで230-1までの自然数についてこれを確認する。まず、32ビットの整数だと簡単にオーバーフローするので、32ビット整数を2つ使って計算する。これは簡単。偶数は直ちにより小さい奇数に辿りつくので、偶数は確認しなくてよい。また、1に辿りつかなくても自分より小さくなればそれ以降は確認しなくてよい。
__device__ bool collatz_core(unsigned int n) {
if(n == 1)
return true;
unsigned int a[2];
a[0] = n;
a[1] = 0;
while(a[0] >= n || a[1] != 0) {
if(a[0] & 1) {
if(a[1] >= 0x7fffffff / 3 * 2 - 3)
return false;
// multiply by 3 / 2
int bit = a[1] & 1;
a[1] += a[1] >> 1;
a[0] += (a[0] >> 1) - 1;
a[0] |= bit << 30;
// add 2
a[0] += 2;
}
else {
// divide by 2
int bit = a[1] & 1;
a[1] >>= 1;
a[0] >>= 1;
a[0] |= bit << 30;
}
int bit = (a[0] >> 31) & 1;
a[0] &= 0x7fffffff;
a[1] += bit;
}
return true;
}__global__ void collatz(int *data, int width, int height, int nRep) {
unsigned int xIndex = blockIdx.x * BLOCK_DIM + threadIdx.x;
unsigned int yIndex = blockIdx.y * BLOCK_DIM + threadIdx.y;
if( (xIndex < width) && (yIndex < height)) {
unsigned int index = yIndex * width + xIndex;
int x = data[index];
for(int i = x; i < x + nRep; i += 2) {
if(!collatz_core(i)) {
data[index] = -i;
break;
}
}
}
}
GPU内で呼ぶ関数は、頭に"__device__"をつける決まりらしい。その関数は、自分より小さくなったらtrueを返し、オーバーフローしそうになったらfalseを返す。falseが返ってきたら、負数にしてCPUに送り返す。
上のコードは2重ループになっているが、前回のように単独のループにしてもほとんど速度に変わりはなかった。
CPU GPU
〜226-1 -> 1304.8ms 58.7ms
〜228-1 -> 5220.5ms 171.0ms
〜230-1 -> 21706.4ms 586.4ms
だいたい40倍近く速くなっている。乗除がないとループの影響なく速くなるのだろうか。
// collatz.cu
#include
#include
#include
#include#define BLOCK_DIM 16
__device__ bool collatz_core(unsigned int n) {
if(n == 1)
return true;
unsigned int a[2];
a[0] = n;
a[1] = 0;
while(a[0] >= n || a[1] != 0) {
if(a[0] & 1) {
if(a[1] >= 0x7fffffff / 3 * 2 - 3)
return false;
// multiply by 3 / 2
int bit = a[1] & 1;
a[1] += a[1] >> 1;
a[0] += (a[0] >> 1) - 1;
a[0] |= bit << 30;
// add 2
a[0] += 2;
}
else {
// divide by 2
int bit = a[1] & 1;
a[1] >>= 1;
a[0] >>= 1;
a[0] |= bit << 30;
}
int bit = (a[0] >> 31) & 1;
a[0] &= 0x7fffffff;
a[1] += bit;
}
return true;
}__global__ void collatz(int *data, int width, int height, int nRep) {
unsigned int xIndex = blockIdx.x * BLOCK_DIM + threadIdx.x;
unsigned int yIndex = blockIdx.y * BLOCK_DIM + threadIdx.y;
if( (xIndex < width) && (yIndex < height)) {
unsigned int index = yIndex * width + xIndex;
int x = data[index];
for(int i = x; i < x + nRep; i += 2) {
if(!collatz_core(i)) {
data[index] = -i;
break;
}
}
}
}void runTest( int argc, char** argv);
int main( int argc, char** argv) {
runTest( argc, argv);
cutilExit(argc, argv);
}void runTest( int argc, char** argv) {
const unsigned int size_x = 256;
const unsigned int size_y = 256;
const unsigned int nData = size_x * size_y;
const unsigned int mem_size = sizeof(int) * nData;
const int nRep = 1024;
unsigned int timer;
cutCreateTimer(&timer);
if( cutCheckCmdLineFlag(argc, (const char**)argv, "device") )
cutilDeviceInit(argc, argv);
else
cudaSetDevice( cutGetMaxGflopsDeviceId() );
// 初期化
int *h_data = (int *)malloc(nData * sizeof(int));
h_data[0] = 0;
for(int i = 0; i < nData; i++) {
h_data[i] = i * nRep + 1;
}
cutStartTimer(timer);
// データをGPUに転送
int *d_data;
cutilSafeCall(cudaMalloc( (void**)&d_data, mem_size));
cutilSafeCall(cudaMemcpy(d_data, h_data, mem_size, cudaMemcpyHostToDevice));
// setup execution parameters
dim3 grid(size_x / BLOCK_DIM, size_y / BLOCK_DIM, 1);
dim3 threads(BLOCK_DIM, BLOCK_DIM, 1);
// GPUで計算実行
collatz<<< grid, threads >>>(d_data, size_x, size_y, nRep);
cudaThreadSynchronize();
cutilCheckMsg("Kernel execution failed");
// データをGPUから転送
cutilSafeCall(cudaMemcpy(h_data, d_data, mem_size, cudaMemcpyDeviceToHost));
float optimizedTime = cutGetTimerValue(timer);
cutStopTimer(timer);
for(int i = 0; i < nData; i++) {
if(h_data[i] < 0) {
printf("%d\n", -h_data[i]);
}
}
printf("time: %0.3f ms\n\n", optimizedTime);
// cleanup memory
free(h_data);
cutilSafeCall(cudaFree(d_data));
cutilCheckError( cutDeleteTimer(timer));
cudaThreadExit();
}
#include
#includebool collatz(unsigned int n) {
if(n == 1)
return true;
unsigned int a[2];
a[0] = n;
a[1] = 0;
while(a[0] >= n || a[1] != 0) {
if(a[0] & 1) {
// if(a[1] >= 0x7fffffff / 3 * 2 - 3)
if(a[1] >= 0x3ffffff / 3 * 2 - 3)
return false;
// multiply by 3 / 2
int bit = a[1] & 1;
a[1] += a[1] >> 1;
a[0] += (a[0] >> 1) - 1;
a[0] |= bit << 30;
// add 2
a[0] += 2;
}
else {
// divide by 2
int bit = a[1] & 1;
a[1] >>= 1;
a[0] >>= 1;
a[0] |= bit << 30;
}
int bit = (a[0] >> 31) & 1;
a[0] &= (1 << 31) - 1;
a[1] += bit;
}
return true;
}int main() {
const unsigned int size_x = 1024;
const unsigned int size_y = 1024;
const unsigned int nData = size_x * size_y;
const unsigned int mem_size = sizeof(int) * nData;
const int nRep = 1024;
int *h_data = (int *)malloc(nData * sizeof(int));
for(int i = 0; i < nData; i++) {
h_data[i] = i * nRep + 1;
}
LARGE_INTEGER nFreq, nStart, nEnd;
QueryPerformanceFrequency(&nFreq);
QueryPerformanceCounter(&nStart);
for(int i = 0; i < nData; i++) {
int x = h_data[i];
for(int n = x; n < x + nRep; n += 2) {
if(!collatz(n))
printf("%d\n", n);
}
}
QueryPerformanceCounter(&nEnd);
printf("time: %0.3f ms\n\n",
(nEnd.QuadPart - nStart.QuadPart) * 1e3 / nFreq.QuadPart);
free(h_data);
}
