1 #include <thrust/random.h> 2 #include <thrust/iterator/counting_iterator.h> 3 #include <thrust/functional.h> 4 #include <thrust/transform_reduce.h> 5 6 #include <iostream> 7 #include <iomanip> 8 #include <cmath> 9 10 // we could vary M & N to find the perf sweet spot 11 12 __host__ __device__ hash(unsigned int a)13unsigned int hash(unsigned int a) 14 { 15 a = (a+0x7ed55d16) + (a<<12); 16 a = (a^0xc761c23c) ^ (a>>19); 17 a = (a+0x165667b1) + (a<<5); 18 a = (a+0xd3a2646c) ^ (a<<9); 19 a = (a+0xfd7046c5) + (a<<3); 20 a = (a^0xb55a4f09) ^ (a>>16); 21 return a; 22 } 23 24 struct estimate_pi : public thrust::unary_function<unsigned int,float> 25 { 26 __host__ __device__ operator ()estimate_pi27 float operator()(unsigned int thread_id) 28 { 29 float sum = 0; 30 unsigned int N = 10000; // samples per thread 31 32 unsigned int seed = hash(thread_id); 33 34 // seed a random number generator 35 thrust::default_random_engine rng(seed); 36 37 // create a mapping from random numbers to [0,1) 38 thrust::uniform_real_distribution<float> u01(0,1); 39 40 // take N samples in a quarter circle 41 for(unsigned int i = 0; i < N; ++i) 42 { 43 // draw a sample from the unit square 44 float x = u01(rng); 45 float y = u01(rng); 46 47 // measure distance from the origin 48 float dist = sqrtf(x*x + y*y); 49 50 // add 1.0f if (u0,u1) is inside the quarter circle 51 if(dist <= 1.0f) 52 sum += 1.0f; 53 } 54 55 // multiply by 4 to get the area of the whole circle 56 sum *= 4.0f; 57 58 // divide by N 59 return sum / N; 60 } 61 }; 62 main(void)63int main(void) 64 { 65 // use 30K independent seeds 66 int M = 30000; 67 68 float estimate = thrust::transform_reduce(thrust::counting_iterator<int>(0), 69 thrust::counting_iterator<int>(M), 70 estimate_pi(), 71 0.0f, 72 thrust::plus<float>()); 73 estimate /= M; 74 75 std::cout << std::setprecision(3); 76 std::cout << "pi is approximately " << estimate << std::endl; 77 78 return 0; 79 } 80 81