（不妨将下标改为从1开始）
参考loj2265中关于杨表的相关知识
构造一个$n$行且第$i$行有$a_{i}$个格子的杨表，依次记录其每一次增加的时间（范围为$[1,\sum_{i=1}^{n}a_{i}]$）
不难发现，条件即变为要求得到的杨表为标准杨表
另一方面，每一个标准杨表都对应一组方案，因此合法方案数即为$f_{a}$
关于$f_{a}$的计算，根据性质3.2，即有$f_{a}=\frac{(\sum_{i=1}^{n}a_{i})!}{\prod_{1\le i\le n,1\le j\le a_{i}}h_{a}(i,j)}$
假设第$j$列最后一个格子在第$k$行，那么$h_{a}(i,j)$即为$(a_{i}-j)+(k-i)+1$
不妨先枚举$i$，然后从1到$a_{i}$依次增大$j$，再不断调整（减小）$k$使得其合法（$a_{k}\ge j$）
在这一过程中，$(j,k)$这个二元组的变化是连续的，即每一次要么$j$减小1、要么$k$增加1，因此$h_{a}(i,j)$也即是一个连续区间，通过最小值和最大值可得其为$(a_{i}+n-i)!$
但这并不是所要求的，因为每一次其实只需要计算最终的$k$，而这将过程中的$k$也计算了
具体的，即每一次$k$将要减小时，都要去掉这一组的贡献，注意到此时必然有$j=a_{k}+1$且从$m$到$i+1$每一个$k$都恰好发生一次，也即为$\prod_{k=i+1}^{n}(a_{i}-a_{k}-1)+(k-i)+1$
将其整理，即记$b_{i}=a_{i}+n-i$，则有$f_{a}=\frac{\prod_{1\le i<j\le n}(b_{i}-b_{j})}{\prod_{i=1}^{n}b_{i}!}(\sum_{i=1}^{n}a_{i})!$
再考虑总方案数，实际上即忽略标准杨表中每一列单调递增的条件，先对其任意排列再对每一行重新排列，不难得到方案数为$\frac{(\sum_{i=1}^{n}a_{i})!}{\prod_{i=1}^{n}a_{i}!}$
将两者相除即为概率，因此答案为$\prod_{1\le i<j\le n}(b_{i}-b_{j})\prod_{i=1}^{n}\frac{a_{i}!}{b_{i}!}$
对前者使用fft计数即可，注意到$b_{i}$严格单调递减，只需要保留次数$>0$的项即可（当然为了避免负次数会让所有次数都加上$a_{1}+n$，因此也即$>a_{1}+n$），后者显然容易计算
总复杂度为$o(n\log n)$，可以通过
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N (1<<21) 4 #define L 21 5 #define mod 1004535809 6 #define PI acos(-1.0) 7 #define ll long long 8 #define cd complex<double> 9 cd a[N],b[N],S[2][L][N]; 10 int n,C,x,ans,fac[N],inv[N],rev[N]; 11 int qpow(int n,ll m){ 12 int s=n,ans=1; 13 while (m){ 14 if (m&1)ans=(ll)ans*s%mod; 15 s=(ll)s*s%mod; 16 m>>=1; 17 } 18 return ans; 19 } 20 void init(){ 21 fac[0]=inv[0]=inv[1]=1; 22 for(int i=1;i<N;i++)fac[i]=(ll)fac[i-1]*i%mod; 23 for(int i=2;i<N;i++)inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod; 24 for(int i=1;i<N;i++)inv[i]=(ll)inv[i-1]*inv[i]%mod; 25 for(int i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)+((i&1)<<L-1); 26 for(int i=0;i<L;i++){ 27 cd s=cd(cos(PI/(1<<i)),sin(PI/(1<<i))),ss=conj(s); 28 S[0][i][0]=S[1][i][0]=cd(1,0); 29 for(int j=1;j<(1<<i);j++){ 30 S[0][i][j]=S[0][i][j-1]*s; 31 S[1][i][j]=S[1][i][j-1]*ss; 32 } 33 } 34 } 35 void fft(cd *a,int p){ 36 for(int i=0;i<N;i++) 37 if (i<rev[i])swap(a[i],a[rev[i]]); 38 for(int i=2,ii=0;i<=N;i<<=1,ii++) 39 for(int j=0;j<N;j+=i) 40 for(int k=0;k<(i>>1);k++){ 41 cd x=a[j+k],y=a[j+k+(i>>1)]*S[p][ii][k]; 42 a[j+k]=x+y,a[j+k+(i>>1)]=x-y; 43 } 44 if (p){ 45 for(int i=0;i<N;i++)a[i]/=N; 46 } 47 } 48 int main(){ 49 init(); 50 scanf("%d",&n); 51 C=1e6,ans=1; 52 for(int i=1;i<=n;i++){ 53 scanf("%d",&x); 54 ans=(ll)ans*fac[x]%mod; 55 x+=n-i; 56 ans=(ll)ans*inv[x]%mod; 57 a[x].real(a[x].real()+1); 58 b[C-x].real(b[C-x].real()+1); 59 } 60 fft(a,0); 61 fft(b,0); 62 for(int i=0;i<N;i++)a[i]*=b[i]; 63 fft(a,1); 64 for(int i=C+1;i<N;i++)ans=(ll)ans*qpow(i-C,floor(a[i].real()+0.5))%mod; 65 printf("%d\n",ans); 66 return 0; 67 }