Birthday Problem Revisited

A long long time ago in a galaxy far far away, the Wimwians, inhabitants of planet WimWi, discovered an unmanned drone that had landed on their planet. On examining the drone, they uncovered a device that sought the answer for the so called "Birthday Problem". The description of the problem was as follows:

If people on your planet were to enter a very large room one by one, what will be the expected number of people in the room when you first find 3 people with Birthdays within 1 day from each other.

The description further instructed them to enter the answer into the device and send the drone into space again. Startled by this turn of events, the Wimwians consulted their best mathematicians. Each year on Wimwi has 10 days and the mathematicians assumed equally likely birthdays and ignored leap years (leap years in Wimwi have 11 days), and found 5.78688636 to be the required answer. As such, the Wimwians entered this answer and sent the drone back into space.

After traveling light years away, the drone then landed on planet Joka. The same events ensued except this time, the numbers in the device had changed due to some unknown technical issues. The description read:

If people on your planet were to enter a very large room one by one, what will be the expected number of people in the room when you first find 3 people with Birthdays within 7 days from each other.

With a 100-day year on the planet, the Jokars (inhabitants of Joka) found the answer to be 8.48967364 (rounded to 8 decimal places because the device allowed only 8 places after the decimal point) assuming equally likely birthdays. They too entered the answer into the device and launched the drone into space again.

This time the drone landed on planet Earth. As before the numbers in the problem description had changed. It read:

If people on your planet were to enter a very large room one by one, what will be the expected number of people in the room when you first find 4 people with Birthdays within 7 days from each other.

What would be the answer (rounded to eight places after the decimal point) the people of Earth have to enter into the device for a year with 365 days? Ignore leap years. Also assume that all birthdays are equally likely and independent of each other.

题解

我们把人走进房间的过程看成在跳状态，那么根据期望的线性性，期望次数就是每个不合法的状态被走到的概率之和。

假设某状态第 \(i\) 天过生日的人数为 \(cnt_i\)，那么这个状态被走到的概率为
\[ \frac{(\sum cnt_i)!}{\prod cnt_i!}(\frac{1}{365})^{\sum cnt_i} \]
考虑用DP计算所有不合法的状态的概率总和。分子不好做，那么我们就在状态里记录一下分子。另外为了满足题目中的限制，我们还需要多加一维表示最近7天里的人数情况。

然而事情并没有这么简单，你写个爆搜发现你怎么都凑不出来第一个样例的5.78688636。

然后我去搜题解，不料发现了这个：https://projecteuler.chat/viewtopic.php?t=6496

Q: Could anyone explain what does mean: "3 people with Birthdays within 7 days from each other"? Does it mean e.g. birthdays Jan 1, Jan 8 and Jan 15 or 3 birthdays within a week like Jan 1, Jan 3 and Jan 7?

A: Actually, it means within a week plus one day, since 0 days apart would mean on the same day.

So if, for example, someone would have its Birthday on some Friday, another person on Friday a week later and someone in between these two, the condition would be fulfilled.

Q: Have another question. Let the year length be 5 days. If I want to get birthdays that are 1 day apart, are they just {(1,2), (2, 3), (3, 4), (4,5)} or should we also count (5,1) (day 5 of previous year and day 1 of the current year)?

A: Yes, you should also count the 'wrap-around' possibilities, so 31st December is 1 day away from 1st January. As you may have guessed, the timing of the release of this problem was at least partially influenced by the new year.

原来同一天有两个人叫“within 0 day”，并且一年365天是个环！

据此容易写出能过第一个样例的暴力程序

__int128 fac[30],pw[30];
LD ans;

void dfs(int x,int sum){
    static int cnt[10];
    if(x==10){
        LL res=fac[sum];
        for(int i=0;i<10;++i) res/=fac[cnt[i]];
        ans+=(LD)res/pw[sum];
        return;
    }
    int tot=x>0?cnt[x-1]:0;
    if(x==9) tot=max(tot,cnt[0]);
    for(int i=0;i<=2-tot;++i) cnt[x]=i,dfs(x+1,sum+i);
}
int main(){
    fac[0]=1;
    for(int i=1;i<=20;++i) fac[i]=fac[i-1]*i;
    pw[0]=1;
    for(int i=1;i<=20;++i) pw[i]=pw[i-1]*10;
    dfs(0,0);
    printf("%.8Lf\n",ans);
    return 0;
}

每个人应该在往前考虑7天的时候同时往后考虑7天。之前之所以可以在序列上只考虑前面7天进行DP，是因为在第 \(t\) 天往后考虑7天等价于在第 \(t+7\) 天往前考虑7天。

在一年是个环的情况下，对最后几天考虑的时候便不能采用上述的懒做法。于是我们需要记录初始7天的人数情况。在考虑最后7天的时候，某一天往后7天有的位置还没有填怎么办？直接算成0显然不对，这一天往后7天里没有填的位置显然有一个总和限制，因此我们还需要额外对最后7天记一个表示人数选择限制的状态。

然后我写出了能过第二个样例的程序

int sta[36+1][7],sum[36+1],idx;
int H[16384];

int code(int a[]){
    int ans=0;
    for(int i=0;i<7;++i) ans=ans*3+a[i];
    return ans;
}
void dfs(int x,int tot){
    static int cnt[7];
    if(x==7){
        ++idx,copy(cnt,cnt+7,sta[idx]),sum[idx]=tot;
        H[code(cnt)]=idx;
        return;
    }
    for(int i=0;i<=2-tot;++i) cnt[x]=i,dfs(x+1,tot+i);
}

//bool vis[100+1][100][36+1];
CO int fac[4]={1,1,2,6};
LD dp[100+1][26+1][36+1][36+1][3];
LD pre[26+1];

int main(){
    dfs(0,0);
    printf("%d\n",idx);
//  vis[0][0][1]=1;
//  for(int i=0;i<100;++i)for(int j=0;j<100;++j)
//      for(int s=1;s<=idx;++s)if(vis[i][j][s])
//          for(int k=0;k<=2-sum[s];++k){
//              static int nw[7];
//              copy(sta[s],sta[s]+6,nw+1),nw[0]=k;
//              vis[i+1][j+k][H[code(nw)]]|=vis[i][j][s];
//          }
//  int mj=0;
//  for(int j=0;j<100;++j)for(int s=1;s<=idx;++s)
//      if(vis[100][j][s]) mj=max(mj,j);
//  printf("%d\n",mj);
    for(int s=1;s<=idx;++s){
        dp[7][sum[s]][s][s][2]=1;
        for(int k=0;k<7;++k) dp[7][sum[s]][s][s][2]/=fac[sta[s][k]];
    }
    for(int i=7;i<100;++i)for(int j=0;j<=26;++j)
        for(int s=1;s<=idx;++s)for(int t=1;t<=idx;++t)
            for(int l=0;l<=2;++l)if(dp[i][j][s][t][l]){
                int tot=sum[s];
                if(i>=93){
                    int ano=0;
                    for(int k=0;k<i-92;++k) ano+=sta[t][6-k];
                    tot=max(tot,ano);
                    for(int k=0;k<=min(2-tot,l);++k){
                        static int nw[7];
                        copy(sta[s],sta[s]+6,nw+1),nw[0]=k;
                        dp[i+1][j+k][H[code(nw)]][t][min(2-k-ano,l-k)]+=dp[i][j][s][t][l]/fac[k];
                    }
                }
                else{
                    for(int k=0;k<=2-tot;++k){
                        static int nw[7];
                        copy(sta[s],sta[s]+6,nw+1),nw[0]=k;
                        dp[i+1][j+k][H[code(nw)]][t][l]+=dp[i][j][s][t][l]/fac[k];
                    }
                }
            }
    pre[0]=1;
    for(int i=1;i<=26;++i) pre[i]=pre[i-1]*i/100;
    LD ans=0;
    for(int j=0;j<=26;++j)for(int s=1;s<=idx;++s)
        for(int t=1;t<=idx;++t)for(int l=0;l<=2;++l)
            if(dp[100][j][s][t][l]) ans+=dp[100][j][s][t][l]*pre[j];
    printf("%.8Lf\n",ans);
    return 0;
}

最后写这道题的第三个实际问题的时候，我发现即使我精打细算开空间，这个五维状态的DP数组还是开不下。解决办法是滚动第一维。

int sta[36+1][7],sum[36+1],idx;
int H[16384];

int code(int a[]){
    int ans=0;
    for(int i=0;i<7;++i) ans=ans*3+a[i];
    return ans;
}
void dfs(int x,int tot){
    static int cnt[7];
    if(x==7){
        ++idx,copy(cnt,cnt+7,sta[idx]),sum[idx]=tot;
        H[code(cnt)]=idx;
        return;
    }
    for(int i=0;i<=2-tot;++i) cnt[x]=i,dfs(x+1,tot+i);
}

//bool vis[100+1][100][36+1];
CO int fac[4]={1,1,2,6};
LD dp[100+1][26+1][36+1][36+1][3];
LD pre[26+1];

int main(){
    dfs(0,0);
    printf("%d\n",idx);
//  vis[0][0][1]=1;
//  for(int i=0;i<100;++i)for(int j=0;j<100;++j)
//      for(int s=1;s<=idx;++s)if(vis[i][j][s])
//          for(int k=0;k<=2-sum[s];++k){
//              static int nw[7];
//              copy(sta[s],sta[s]+6,nw+1),nw[0]=k;
//              vis[i+1][j+k][H[code(nw)]]|=vis[i][j][s];
//          }
//  int mj=0;
//  for(int j=0;j<100;++j)for(int s=1;s<=idx;++s)
//      if(vis[100][j][s]) mj=max(mj,j);
//  printf("%d\n",mj);
    for(int s=1;s<=idx;++s){
        dp[7][sum[s]][s][s][2]=1;
        for(int k=0;k<7;++k) dp[7][sum[s]][s][s][2]/=fac[sta[s][k]];
    }
    for(int i=7;i<100;++i)for(int j=0;j<=26;++j)
        for(int s=1;s<=idx;++s)for(int t=1;t<=idx;++t)
            for(int l=0;l<=2;++l)if(dp[i][j][s][t][l]){
                int tot=sum[s];
                if(i>=93){
                    int ano=0;
                    for(int k=0;k<i-92;++k) ano+=sta[t][6-k];
                    tot=max(tot,ano);
                    for(int k=0;k<=min(2-tot,l);++k){
                        static int nw[7];
                        copy(sta[s],sta[s]+6,nw+1),nw[0]=k;
                        dp[i+1][j+k][H[code(nw)]][t][min(2-k-ano,l-k)]+=dp[i][j][s][t][l]/fac[k];
                    }
                }
                else{
                    for(int k=0;k<=2-tot;++k){
                        static int nw[7];
                        copy(sta[s],sta[s]+6,nw+1),nw[0]=k;
                        dp[i+1][j+k][H[code(nw)]][t][l]+=dp[i][j][s][t][l]/fac[k];
                    }
                }
            }
    pre[0]=1;
    for(int i=1;i<=26;++i) pre[i]=pre[i-1]*i/100;
    LD ans=0;
    for(int j=0;j<=26;++j)for(int s=1;s<=idx;++s)
        for(int t=1;t<=idx;++t)for(int l=0;l<=2;++l)
            if(dp[100][j][s][t][l]) ans+=dp[100][j][s][t][l]*pre[j];
    printf("%.8Lf\n",ans);
    return 0;
}