Skip to content

Instantly share code, notes, and snippets.

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

MonteCarlo方法求PI，内置两个样本实验。Cesaro每次选取两个随机数求其最大公约数，若为1则计数器加1，最终计数器的值除以总的尝试次数近似于6/PI*PI。Area则是在以(0,0)为中心点，以1为边长的正方形中随机抽取一个坐标点，若其位于以(0,0)为圆心，1为半径的圆中，则计数器加1，最终计数器的值除以总的尝试次数近似于PI/4。
 #include #include #include #include //------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
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.