（不妨将下标改为从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 } 

View Code