Linux下多線程通過蒙特卡洛法來求取pi值

特卡洛法又稱隨機抽樣技術

是一種應用隨機數進行仿真試驗的方法。

用該方法計算π的基本思路是：

根據圓面積的公式： 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;
}