特卡洛法又称随机抽样技术
是一种应用随机数进行仿真试验的方法。
用该方法计算π的基本思路是:
根据圆面积的公式: s=πR2 ,当R=1时,S=π。
由于圆的方程是:x2+y2=1(X2为X的平方的意思),因此1/4圆面积为X轴、y轴和上述方程所包围的部分。
如果在1*1的矩形中均匀地落入随机点,则落入1/4园中的点的概率就是1/4圆的面积。其4倍,就是圆面积。
由于半径为1,该面积的值为π的值。
#include <pthread.h>
#include <sys/time.h>
#include <stdlib.h>
#include <stdio.h>
#define MaxThreadNum 32
#define kSamplePoints 1000
#define kSpace 1
void *compute_pi(void *);
inline double WallTime();
int total_hits, hits[MaxThreadNum][kSpace];
int sample_points_per_thread, num_threads;
int main(void)
{
int i;
double time_start, time_end;
pthread_t p_threads[MaxThreadNum];
pthread_attr_t attr;
pthread_attr_init(&attr);
pthread_attr_setscope(&attr, PTHREAD_SCOPE_SYSTEM);
printf("Enter num_threads\n");
scanf("%d", &num_threads);
time_start = WallTime();
total_hits = 0;
sample_points_per_thread = kSamplePoints / num_threads;
for(i=0; i<num_threads; i++)
{
hits[i][0] = i;
pthread_create(&p_threads[i], &attr, compute_pi, (void *)&hits[i]);
}
for(i=0; i<num_threads; i++)
{
pthread_join(p_threads[i], NULL);
total_hits += hits[i][0];
}
double pi = 4.0 * (double)total_hits / kSamplePoints;
time_end = WallTime();
printf("Elasped time: %lf, Pi: %lf\n", time_end - time_start, pi);
return 0;
}
void *compute_pi(void * s)
{
unsigned int seed;
int i;
int *hit_pointer;
double rand_no_x, rand_no_y;
hit_pointer = (int *)s;
seed = *hit_pointer;
local_hits = 0;
for(i=0; i<sample_points_per_thread; i++)
{
rand_no_x = (double)(rand_r(&seed))/(double)(RAND_MAX);
rand_no_y = (double)(rand_r(&seed))/(double)(RAND_MAX);
if((rand_no_x - 0.5)*(rand_no_x - 0.5) + (rand_no_y - 0.5) * (rand_no_y - 0.5) < 0.25)
{
(*hit_pointer)++;
}
seed *= i;
}
pthread_exit(0);
}
inline double WallTime()
{
struct timeval tv;
struct timezone tz;
gettimeofday(&tv, &tz);
double currTime = (double)tv.tv_sec + (double)tv.tv_usec/1000000.0;
return currTime;
}