MonteCarlo方法求PI，内置两个样本实验。Cesaro每次选取两个随机数求其最大公约数，若为1则计数器加1，最终计数器的值除以总的尝试次数近似于6/PI*PI。Area则是在以(0,0)为中心点，以1为边长的正方形中随机抽取一个坐标点，若其位于以(0,0)为圆心，1为半径的圆中，则计数器加1，最终计数器的值除以总的尝试次数近似于PI/4。

monte_carlo.c

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "mpi.h"

//------Common base------
typedef int bool;
enum { false, true };

struct MonteCarlo{
    int trials;
    bool (*experiment)(void);
};

double monteCarloIter(struct MonteCarlo monteCarlo, int trialsRemaining, int trialsPassed)
{
    if(trialsRemaining == 0) return (double)trialsPassed/monteCarlo.trials;
    else if((*monteCarlo.experiment)() == true)
        return monteCarloIter(monteCarlo, trialsRemaining-1, trialsPassed+1);
    else
        return monteCarloIter(monteCarlo, trialsRemaining-1, trialsPassed);
}

double monteCarloWith(int trials, bool (*experiment)(void))
{
    struct MonteCarlo monteCarlo;
    monteCarlo.trials = trials;
    monteCarlo.experiment = experiment;
    return monteCarloIter(monteCarlo, trials, 0);
}

//------cesaro------
int gcd(int a, int b)
{
    int temp;
    while (b != 0)
    {
        temp = a % b;

        a = b;
        b = temp;
    }
    return a;
}

bool cesaroTest()
{
    if(gcd(rand(), rand()) == 1) return true;
    else return false;
}

//------area------
bool areaTest()
{
    //1*1的正方形
    double x = (double)rand()/RAND_MAX;
    double y = (double)rand()/RAND_MAX;
    double z = x*x + y*y;
    if(z <= 1) return true;
    else return false;
}

//------interface------
double estimatePi(int method, int trials)
{
    if(method==0) return sqrt(6/monteCarloWith(trials, &cesaroTest));
    else if(method==1) return 4*monteCarloWith(trials, &areaTest);
    else return 0.0;
}

int main(int argc, char *argv[])
{
    srand(time(NULL)); //以当前时间为参数初始化随机数种子
    int rank;
    int size;
    double result;
    double sum;

    MPI_Init( 0, 0 );
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    if(argc!=3) {
      printf("参数不足");
    } else {
      int method = atoi(argv[1]);
      int trials = atoi(argv[2]);
      result = estimatePi(method, trials);
      MPI_Reduce(&result, &sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

      if(rank == 0) {
        printf("执行总次数为 %d, 结果为 %0.4f \n", trials*size, sum/size);
      }
    }


    MPI_Finalize();

    return 0;
}

/*------ Usage ------
Linux:
    gcc -Wall monte_carlo.c -o MonteCarlo -lm
    ./MonteCarlo
*/