生日悖论

今天http://weibo.com/2887339314/BcqXD9OKz 发了个问题：300个人，至少5个人同一天生
日的概率，博主用蒙特卡罗方法算出来了结果，一时兴起，写了个算精确结果的代码
（似乎结果有点问题。另外浮点数可能有精度问题）。


考虑计算问题的反面：至多4个人在同一天生日，难度在于可能性太多。所以考虑状态压缩。
用一个5元组。(A[0], A[1], A[2], A[3], A[4])表示这365天中有A[0]天，没有不论什么人的生日。
有A[1]天。这些天中有一个人的生日，其他相似。5元组被映射到一个浮点数上，表示相应
事件发生的概率。

从N=1起。不断更新5元组。
N = 1的时候，(364, 1, 0, 0, 0) = 1.0
N = 2的时候。(363, 2, 0, 0, 0) = 364./365 (364, 0, 1, 0, 0) = 1./365
N = 3的时候。(362, 3, 0, 0, 0) = 364./365*363/365
             (363, 1, 1, 0, 0) = 1./365**363/365
             (363, 1, 1, 0, 0) = 364./365*2/365
             (364, 0, 0, 1, 0) = 1./365*1/365
             当中同样的tuple出现多次表示每次出现的值应该累加。

以此类推，考虑到5元组中的元素和一定是365，所以能够将5元组化为4元组。

// 递推计算
#include <cstdio>
#include <iostream>
#include <unordered_map>
using namespace std;
typedef long long int64;

const int AT_LEAST = 5;
const int AT_MOST = AT_LEAST - 1;
const int N = 300;
const int DAY_OF_YEAR = 365;
const int64 one = 1;
const int64 two = 1LL << 9;
const int64 three = 1LL << 18;
const int64 four = 1LL << 27;
const int64 five = 1LL << 36;

const int64 FLAG[] = {0, one, two, three, four, five};

int main()
{
	unordered_map<int64, double> mem;
	mem[1] = 1;
	
	for (int i = 2; i <= N; ++i)
	{
		cerr << i << " " << mem.size() << endl;
		unordered_map<int64, double> next;
		for (auto& curr : mem)
		{
			int how[AT_MOST+1];
			int sum = 0;
			for (int64 val = curr.first, i = 1; i <= AT_MOST; ++i)
			how[i] = val & 511, val >>= 9, sum += how[i];
			
			next[curr.first+one] += curr.second * (DAY_OF_YEAR - sum) / DAY_OF_YEAR;
			
			for (int i = 1; i < AT_MOST; ++i) if (how[i] > 0)
			{
				next[curr.first-FLAG[i]+FLAG[i+1]] += curr.second * how[i] / DAY_OF_YEAR;
			}
		}
		next.swap(mem);
	}
	
	double ans = 0;
	for (auto& curr : mem) ans += curr.second;
	printf("%.16f\n", 1-ans);
	
	return 0;
}

// 蒙特卡罗验证
#include <cstdio>
#include <cstdlib>
#include <ctime>
using namespace std;

int main()
{
	int total = 0;
	srand(time(NULL));
	for (int i = 0; i < 1000000; ++i)
	{
		int data[365] = {0};
		for (int j = 0; j < 300; ++j)
		++data[rand()%365];
		int ok = 0;
		for (int j = 0; j < 365; ++j)
		if (data[j] >= 5) {ok = 1; break;}
		total += ok;
	}
	printf("%.16f\n", 1.*total / 1000000);
	return 0;
}