Description
Solution
- 我们可以模拟这个过程。
- 一种简单的方法是对于奇数长度段直接分，对于偶数长度的段有两个位置可以选，我们可以钦定它选前一个位置，然后再继续往下分，最后再以这个位置在这个段内将所有的在这之后的人的概率对称过去。
- 具体来说，我们可以发现对于所有的当前的段，最大且相同的一些段要一起做，考虑一共有个，那么接下来的个人都填这些段。
- 但是我们注意到，由于偶数位置有两个，而且这两个在所有的位置中都要占概率，所以要用DP来计算，钦定一个位置（奇/偶）不选，剩下的位置选了i个奇，j个偶的概率是，然后就可以计算出这个中的第个在某个奇或某个偶的位置的贡献了。
- 最后对称过去就好了。
- 由于这个分治是类似线段树的，所以一共有层，每一层最坏是的（实际上远远到不了），最后的对称也是的，因为只有在某个操作之后的人才要被这个操作对称。
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#define ll long long
#define maxn 1005
using namespace std;
int n,mo,c[maxn],a[maxn][maxn][2],i,j,k,l,now,p[maxn],q[maxn],nowc,nc[maxn];
ll inv[maxn],ans[maxn][maxn],f[maxn][maxn],g[maxn][maxn];
ll ksm(ll x,ll y){
ll s=1;
for(;y;y/=2,x=x*x%mo) if (y&1)
s=s*x%mo;
return s;
}
int main(){
scanf("%d%d",&n,&mo);
for(i=1;i<=n;i++) inv[i]=ksm(i,mo-2);
c[0]=1,a[0][1][0]=1,a[0][1][1]=n,now=nowc=0;
while (1){
int mx=0;
for(i=1;i<=c[now];i++) mx=max(mx,(a[now][i][1]-a[now][i][0])/2);
if (mx==0){
int cnt=0;
for(i=1;i<=c[now];i++)
cnt+=a[now][i][1]-a[now][i][0]+1;
for(i=1;i<=c[now];i++)
for(j=a[now][i][0];j<=a[now][i][1];j++)
for(k=1;k<=cnt;k++)
ans[nowc+k][j]+=inv[cnt];
break;
}
p[0]=q[0]=0;
for(i=1;i<=c[now];i++) if ((a[now][i][1]-a[now][i][0])/2==mx){
k=(a[now][i][1]-a[now][i][0])/2;
if ((a[now][i][1]-a[now][i][0]+1)&1)
p[++p[0]]=a[now][i][0]+k;
else q[++q[0]]=a[now][i][0]+k;
c[now+1]++;
a[now+1][c[now+1]][0]=a[now][i][0];
a[now+1][c[now+1]][1]=a[now][i][0]+k-1;
c[now+1]++;
a[now+1][c[now+1]][0]=a[now][i][0]+k+1;
a[now+1][c[now+1]][1]=a[now][i][1];
} else {
c[now+1]++;
a[now+1][c[now+1]][0]=a[now][i][0];
a[now+1][c[now+1]][1]=a[now][i][1];
}
for(i=0;i<=p[0];i++) for(j=0;j<=q[0];j++) f[i][j]=g[i][j]=0;
f[0][0]=g[0][0]=1;
for(k=1;k<=p[0]+q[0];k++){
ll sum0=0,sum1=0;
for(i=max(k-1-q[0],0);i<=min(p[0],k-1);i++){ j=k-1-i;
int res=p[0]-i+2*(q[0]-j);
if (i<p[0]) sum0+=f[i][j]*inv[res]%mo;
if (j<q[0]) sum1+=g[i][j]*2*inv[res]%mo;
}
sum0%=mo,sum1%=mo;
for(i=1;i<=p[0];i++) ans[nowc+k][p[i]]+=sum0;
for(i=1;i<=q[0];i++) ans[nowc+k][q[i]]+=sum1;
if (k==p[0]+q[0]) break;
for(i=max(k-1-q[0],0);i<=min(p[0],k-1);i++){j=k-1-i;
int res=p[0]-i+2*(q[0]-j);
if (i<p[0]) {
f[i+1][j]+=f[i][j]*(p[0]-i-1)%mo*inv[res]%mo;
g[i+1][j]+=g[i][j]*(p[0]-i)%mo*inv[res]%mo;
}
if (j<q[0]) {
f[i][j+1]+=f[i][j]*2*(q[0]-j)%mo*inv[res]%mo;
g[i][j+1]+=g[i][j]*2*(q[0]-j-1)%mo*inv[res]%mo;
}
}
}
nowc+=p[0]+q[0],nc[++now]=nowc;
}
for(now--;now>=0;now--){
for(i=1;i<=c[now];i++) if (!((a[now][i][1]-a[now][i][0]+1)&1)){
int mid=(a[now][i][0]+a[now][i][1])/2;
for(j=a[now][i][0];j<=mid;j++){
l=a[now][i][1]-(j-a[now][i][0]);
for(k=nc[now]+1;k<=n;k++){
ll s=(ans[k][j]+ans[k][l])*inv[2]%mo;
ans[k][j]=ans[k][l]=s;
}
}
}
}
for(i=1;i<=n;i++,printf("
"))
for(j=1;j<=n;j++)
printf("%lld ",ans[i][j]);
}