Skip to content

Instantly share code, notes, and snippets.

@LuckyKoala LuckyKoala/monte_carlo.c
Last active Jun 18, 2018

Embed
What would you like to do?
MonteCarlo方法求PI,内置两个样本实验。Cesaro每次选取两个随机数求其最大公约数,若为1则计数器加1,最终计数器的值除以总的尝试次数近似于6/PI*PI。Area则是在以(0,0)为中心点,以1为边长的正方形中随机抽取一个坐标点,若其位于以(0,0)为圆心,1为半径的圆中,则计数器加1,最终计数器的值除以总的尝试次数近似于PI/4。
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.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(void)
{
srand(time(NULL)); //以当前时间为参数初始化随机数种子
for(;;)
{
int method, trials;
printf("------ MonteCarlo PI ------ \n");
printf("样本实验代号(-1为退出, 0为Cesaro,1为Area)\n");
printf("参数示例: 0 10000 即是选择Cesaro样本实验并尝试10000次 \n");
printf("请输入参数: \n");
scanf("%d %d", &method, &trials);
if(method == -1) break;
printf("本轮计算结果为 %f\n", estimatePi(method, trials));
}
return 0;
}
//------ Usage ------
//Linux:
// gcc -Wall monte_carlo.c -o MonteCarlo -lm
// ./MonteCarlo
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.