分析：
我们发现每个人操作时，不管前面的人怎么挂，当前这个人的\(d\)为定值
那么我们可以尝试以某种特定的方案模拟整个过程，再通过一些变换求出概率
假设目前的最大距离为\(d\)，可以挂的空位置有\(m\)段，说明从这个人开始的后\(m\)个人都是这个\(d\)，他们的概率可以一起算
\(m\)段的长度有奇数和偶数，奇数时可以直接挂在中点，偶数时中间两个点都可以挂
假设p段奇数，q段偶数，我们可以列出DP式\(f_{i,j,k}\)表示第\(i\)个人要挂奇数段时，奇数段剩\(j\)个，偶数段剩\(k\)个的概率
由于\(i+j+k=m\)，\(i\)可以省去
用一个奇数段的概率为\(\frac{j}{j+2k}\)，偶数段为\(\frac{2k}{j+2k}\)
表示第\(i\)个人要挂偶数段时，设一个DP式\(g_{i,j}\)同样处理
复杂度为\(O(n^2)\)
但是一个人选择挂偶数段时，两个点各自会产生新的情况，但是我们发现在这一段当中，挂左边和挂右边产生的情况是一一对称的
我们强行确定挂左边，并在中间记录一个对称轴，全部处理完之后，对称的位置取一次平均数就好了
特别注意，当最后\(d=1\)时，剩下的位置都可以选并且互不影响，不要进行上面的DP，直接暴力处理并结束就好了
复杂度式子冗杂，但整体上是不超过\(O(n^2)\)的
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<iostream>
#include<map>
#include<string>
#define maxn 1005
using namespace std;
inline int getint()
{
int num=0,flag=1;char c;
while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;
while(c>='0'&&c<='9')num=num*10+c-48,c=getchar();
return num*flag;
}
int n,MOD;
int P[maxn],cur;
int vis[maxn];
int T[maxn*maxn],L[maxn*maxn],M[maxn*maxn],R[maxn*maxn],cnt;
int ans[maxn][maxn];
int f[maxn][maxn],g[maxn][maxn];
int inv[maxn];
inline int ksm(int num,int k)
{
int ret=1;
for(;k;k>>=1,num=1ll*num*num%MOD)if(k&1)ret=1ll*ret*num%MOD;
return ret;
}
int main()
{
n=getint(),MOD=getint();inv[0]=1;
for(int i=1;i<=n;i++)inv[i]=ksm(i,MOD-2);
for(int t=1;t<=n;t++)
{
P[cur=0]=0;
for(int i=1;i<=n;i++)if(vis[i])P[++cur]=i;
P[++cur]=n+1;
int mx=0,tim=0,tot[2]={0};
for(int i=cur;i;i--)P[i]=P[i]-P[i-1];
for(int i=1;i<=cur;i++)
if(P[i]/2>mx)mx=P[i]/2,tim=1+(P[i]&1),tot[0]=tot[1]=0,tot[P[i]&1]++;
else if(P[i]/2==mx)tim+=1+(P[i]&1),tot[P[i]&1]++;
int p=ksm(tim,MOD-2);
if(mx==1)
{
for(int i=1;i<=n;i++)if(!vis[i])ans[t][i]=p;
continue;
}
for(int i=0;i<=tot[0];i++)for(int j=0;j<=tot[1];j++)f[i][j]=g[i][j]=0;
f[0][0]=g[0][0]=1;
for(int tt=0;tt<tot[0]+tot[1];tt++)
{
int sumf=0,sumg=0;
for(int i=max(0,tt-tot[1]);i<=min(tt,tot[0]);i++)
{
int j=tt-i,num=tot[0]-i+2*(tot[1]-j);
sumf=(sumf+1ll*f[i][j]*inv[num])%MOD;
sumg=(sumg+1ll*g[i][j]*inv[num])%MOD;
if(i<tot[0])
f[i+1][j]=(f[i+1][j]+1ll*f[i][j]*(tot[0]-i-1)%MOD*inv[num])%MOD,
g[i+1][j]=(g[i+1][j]+1ll*g[i][j]*(tot[0]-i)%MOD*inv[num])%MOD;
if(j<tot[1])
f[i][j+1]=(f[i][j+1]+1ll*f[i][j]*2*(tot[1]-j)%MOD*inv[num])%MOD,
g[i][j+1]=(g[i][j+1]+1ll*g[i][j]*2*(tot[1]-j-1)%MOD*inv[num])%MOD;
}
int now=0;
for(int i=1;i<=cur;now+=P[i++])if(P[i]/2==mx)
{
if(P[i]&1)
{
ans[t+tt][now+mx]=ans[t+tt][now+mx+1]=sumg;
T[++cnt]=t+1,L[cnt]=now+1,M[cnt]=now+mx,R[cnt]=now+P[i]-1;
}
else ans[t+tt][now+mx]=sumf;
vis[now+mx]=1;
}
}
t+=tot[0]+tot[1]-1;
}
for(int i=cnt;i;i--)
{
for(int j=T[i];j<=n;j++)for(int k=L[i];k<=M[i];k++)
ans[j][k]=ans[j][R[i]+L[i]-k]=1ll*(ans[j][k]+ans[j][R[i]+L[i]-k])*inv[2]%MOD;
}
for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)printf("%d%c",ans[i][j],j==n?'\n':' ');
}
