// kadai1-2.cpp
// ex) plot "???.txt" w boxes
#define _USE_MATH_DEFINES
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
 
//ランダムウォークのステップ数
const int STEPS = 1000;
 
const int RANDOMS = 1000; //生成する乱数の個数
const int LOWER_LIMIT = 0; //ヒストグラムの下限と
const int UPPER_LIMIT = 99; //上限
const int DIVISION = 1; //分割数（1の幅をDIVISION等分）
 
//0以上1未満の一様乱数を生成する
double frand()
{
	return rand() / (RAND_MAX+1.0);
}
 
//Box-Muller法で標準正規乱数を生成する
double box_muller()
{
	//Box-Muller法は2個の一様乱数から2個の正規乱数を生成する
	//そのうち一方はすぐに返し、もう一方は次に呼び出されたときに返す
	static bool haveRandom = false;
	static double random;
 
	if(haveRandom){
		//正規乱数はすでにできている
		haveRandom = false;
		return random;
	} else {
		//正規乱数は手元にないので新しく2個の正規乱数を作る
		double u1 = sqrt(-2*log(1-frand()));
		double u2 = 2*M_PI*frand();
 
		//2個のうち一方は次の呼び出しのためにとっておく
		haveRandom = true;
		random = u1 * sin(u2);
 
		//もう一方を返す
		return u1 * cos(u2);
	}
}
 
int main()
{
	double	x, y, d;
	int	i, step, h;
 
	//乱数を初期化する
	srand((unsigned int)time(NULL));
 
	//頻度を記録する配列、0で初期化
	int hist[(UPPER_LIMIT - LOWER_LIMIT) * DIVISION + 1] = {0};
 
	for (i = 0; i < RANDOMS; i++) {
		//STEPSステップのランダムウォークを実施する
		x = y = 0.0;
		for (step = 0; step < STEPS; step++) {
			x += box_muller();
			y += box_muller();
		}
		d = sqrt(x * x + y * y);
		h = (int)d;
		if (h <= UPPER_LIMIT) {
			hist[h]++;
		}
	}
	for (h = LOWER_LIMIT; h <= UPPER_LIMIT; h++) {
		printf("%3d %d\n", h, hist[h]);
	}
 
	return 0;
}
 

Ly8ga2FkYWkxLTIuY3BwCi8vIGV4KSBwbG90ICI/Pz8udHh0IiB3IGJveGVzCiNkZWZpbmUgX1VTRV9NQVRIX0RFRklORVMKI2luY2x1ZGUgPG1hdGguaD4KI2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPHRpbWUuaD4KCi8v44Op44Oz44OA44Og44Km44Kp44O844Kv44Gu44K544OG44OD44OX5pWwCmNvbnN0IGludCBTVEVQUyA9IDEwMDA7Cgpjb25zdCBpbnQgUkFORE9NUyA9IDEwMDA7IC8v55Sf5oiQ44GZ44KL5Lmx5pWw44Gu5YCL5pWwCmNvbnN0IGludCBMT1dFUl9MSU1JVCA9IDA7IC8v44OS44K544OI44Kw44Op44Og44Gu5LiL6ZmQ44GoCmNvbnN0IGludCBVUFBFUl9MSU1JVCA9IDk5OyAvL+S4iumZkApjb25zdCBpbnQgRElWSVNJT04gPSAxOyAvL+WIhuWJsuaVsO+8iDHjga7luYXjgpJESVZJU0lPTuetieWIhu+8iQoKLy8w5Lul5LiKMeacqua6gOOBruS4gOanmOS5seaVsOOCkueUn+aIkOOBmeOCiwpkb3VibGUgZnJhbmQoKQp7CglyZXR1cm4gcmFuZCgpIC8gKFJBTkRfTUFYKzEuMCk7Cn0KCi8vQm94LU11bGxlcuazleOBp+aomea6luato+imj+S5seaVsOOCkueUn+aIkOOBmeOCiwpkb3VibGUgYm94X211bGxlcigpCnsKCS8vQm94LU11bGxlcuazleOBrzLlgIvjga7kuIDmp5jkubHmlbDjgYvjgoky5YCL44Gu5q2j6KaP5Lmx5pWw44KS55Sf5oiQ44GZ44KLCgkvL+OBneOBruOBhuOBoeS4gOaWueOBr+OBmeOBkOOBq+i/lOOBl+OAgeOCguOBhuS4gOaWueOBr+asoeOBq+WRvOOBs+WHuuOBleOCjOOBn+OBqOOBjeOBq+i/lOOBmQoJc3RhdGljIGJvb2wgaGF2ZVJhbmRvbSA9IGZhbHNlOwoJc3RhdGljIGRvdWJsZSByYW5kb207CgoJaWYoaGF2ZVJhbmRvbSl7CgkJLy/mraPopo/kubHmlbDjga/jgZnjgafjgavjgafjgY3jgabjgYTjgosKCQloYXZlUmFuZG9tID0gZmFsc2U7CgkJcmV0dXJuIHJhbmRvbTsKCX0gZWxzZSB7CgkJLy/mraPopo/kubHmlbDjga/miYvlhYPjgavjgarjgYTjga7jgafmlrDjgZfjgY8y5YCL44Gu5q2j6KaP5Lmx5pWw44KS5L2c44KLCgkJZG91YmxlIHUxID0gc3FydCgtMipsb2coMS1mcmFuZCgpKSk7CgkJZG91YmxlIHUyID0gMipNX1BJKmZyYW5kKCk7CgoJCS8vMuWAi+OBruOBhuOBoeS4gOaWueOBr+asoeOBruWRvOOBs+WHuuOBl+OBruOBn+OCgeOBq+OBqOOBo+OBpuOBiuOBjwoJCWhhdmVSYW5kb20gPSB0cnVlOwoJCXJhbmRvbSA9IHUxICogc2luKHUyKTsKCgkJLy/jgoLjgYbkuIDmlrnjgpLov5TjgZkKCQlyZXR1cm4gdTEgKiBjb3ModTIpOwoJfQp9CgppbnQgbWFpbigpCnsKCWRvdWJsZQl4LCB5LCBkOwoJaW50CWksIHN0ZXAsIGg7CgoJLy/kubHmlbDjgpLliJ3mnJ/ljJbjgZnjgosKCXNyYW5kKCh1bnNpZ25lZCBpbnQpdGltZShOVUxMKSk7CgoJLy/poLvluqbjgpLoqJjpjLLjgZnjgovphY3liJfjgIEw44Gn5Yid5pyf5YyWCglpbnQgaGlzdFsoVVBQRVJfTElNSVQgLSBMT1dFUl9MSU1JVCkgKiBESVZJU0lPTiArIDFdID0gezB9OwoKCWZvciAoaSA9IDA7IGkgPCBSQU5ET01TOyBpKyspIHsKCQkvL1NURVBT44K544OG44OD44OX44Gu44Op44Oz44OA44Og44Km44Kp44O844Kv44KS5a6f5pa944GZ44KLCgkJeCA9IHkgPSAwLjA7CgkJZm9yIChzdGVwID0gMDsgc3RlcCA8IFNURVBTOyBzdGVwKyspIHsKCQkJeCArPSBib3hfbXVsbGVyKCk7CgkJCXkgKz0gYm94X211bGxlcigpOwoJCX0KCQlkID0gc3FydCh4ICogeCArIHkgKiB5KTsKCQloID0gKGludClkOwoJCWlmIChoIDw9IFVQUEVSX0xJTUlUKSB7CgkJCWhpc3RbaF0rKzsKCQl9Cgl9Cglmb3IgKGggPSBMT1dFUl9MSU1JVDsgaCA8PSBVUFBFUl9MSU1JVDsgaCsrKSB7CgkJcHJpbnRmKCIlM2QgJWRcbiIsIGgsIGhpc3RbaF0pOwoJfQoKCXJldHVybiAwOwp9Cg==