モンテカルロ法で円周率を求める

二次元平面上、0 < x <= 1, 0 < y <= 1 の領域に、デタラメに点を打ちます。その点と原点(0,0)との距離が1以下であるなら、点は1/4円の内部にあります。

大量の点　(x[i],y[i]) : i = 0, 1...N-1 をデタラメに打ち、1/4円の中に入った点の数をPとすると、P/N はだいたい π/4 になるハズね。CUDAでやってみましょうね。

// CUDA Runtime
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

// cuRAND
#include <curand.h>

// Thrust
#include <thrust/device_ptr.h>
#include <thrust/reduce.h>

#include <iostream>
using namespace std;

__global__ void kernel_is_in_circle(const float* x, const float* y, 
                                    float* count, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    count[i] = ( x[i]*x[i] + y[i]*y[i] < 1.0f ) ? 1.0f : 0.0f; 
  }
}

int main() {
  const unsigned int N = 1000 * 1000;

  float* x = nullptr;
  float* y = nullptr;
  float* count = nullptr;

  cudaMalloc(&x,     N*sizeof(float));
  cudaMalloc(&y,     N*sizeof(float));
  cudaMalloc(&count, N*sizeof(float));

  curandGenerator_t gen;
  // メルセンヌ・ツイスタ
  curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_MTGP32);
  curandSetPseudoRandomGeneratorSeed(gen, 1234ULL); // seedはテキトーな値

  // 0～1の一様乱数列を生成する
  cout << "--- generating uniform random numbers ---" << endl;
  curandGenerateUniform(gen, x, N);
  curandGenerateUniform(gen, y, N);

  // i番目の点(x[i],y[i])が1/4円の中に入っていれば count[i] = 1
  cout << "--- is the point in the circle? ---" < endl;
  unsigned int block = 256;
  unsigned int grid = (N + block -1) / block;
  kernel_is_in_circle<<<grid,block>>>(x, y, count, N);

  // count[N]の総計をNで割って4倍すれば π になるハズ
  cout << "--- accumulating the result(pi) ---" << endl;
  thrust::device_ptr<float> dev(count);
  float pi = thrust::reduce(dev, dev+N) * 4.0f / N;
  cout << pi << endl;

  curandDestroyGenerator(gen);

  cudaFree(x);
  cudaFree(y);
  cudaFree(count);

}

百万個の点を打って勘定したら 3.14116、まぁいいセンいってんじゃないでしょか。

(original: 2015-04-06)