モンテカルロ法で円周率を求める。

モンテカルロ法で円周率を求めるプログラムを簡単に書いてみた。
円周率の計算にはもっとよい方法があるのでこの方法は使われていない。

1辺の長さが1の正方形とそれに内接する4分の1円(扇形)がある。

正方形の面積 : 扇形の面積 = 1 : π / 4

この正方形のなかに，ランダムにn個の点を落としたときr個の点が円内に落ちたとすれば，

1 : π / 4 = n : r
π = 4 * r / n

となる。
落とすn個の点を多くするほど正確な値に近づく。
円周率の正確さは落とす点の個数よりも乱数が一様であるかに依存するらしい。

入力は落とす点の個数n。

#include<iostream>
#include<vector>
#include<algorithm>
#include<cstdio>
#include<climits>
#include<cmath>

using namespace std;

typedef unsigned int uint;

//乱数 XOR Shift
uint xor128(void) { 
  static uint x = 123456789;
  static uint y = 362436069;
  static uint z = 521288629;
  static uint w = 88675123; 
  uint t;
  
  t=x^(x<<11);
  x=y; y=z; z=w;
  return w=(w^(w>>19))^(t^(t>>8)); 
}

int main(void){
	double n,r=0;
	
	cin >> n;
	
	for(int i=0;i<n;i++){
		double x=(double)xor128()/UINT_MAX;
		double y=(double)xor128()/UINT_MAX;
		
		if(sqrt(x*x+y*y)<1.0)r++;
	}
	
	printf("%.8f\n",4.0*r/n);
	
	return 0;
}