[loj6051]PATH

（不妨将下标改为从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)$，可以通过

#include<bits/stdc++.h>
using namespace std;
#define N (1<<21)
#define L 21
#define mod 1004535809
#define PI acos(-1.0)
#define ll long long
#define cd complex<double>
cd a[N],b[N],S[2][L][N];
int n,C,x,ans,fac[N],inv[N],rev[N];
int qpow(int n,ll m){
    int s=n,ans=1;
    while (m){
        if (m&1)ans=(ll)ans*s%mod;
        s=(ll)s*s%mod;
        m>>=1;
    }
    return ans;
}
void init(){
    fac[0]=inv[0]=inv[1]=1;
    for(int i=1;i<N;i++)fac[i]=(ll)fac[i-1]*i%mod;
    for(int i=2;i<N;i++)inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod;
    for(int i=1;i<N;i++)inv[i]=(ll)inv[i-1]*inv[i]%mod;
    for(int i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)+((i&1)<<L-1);
    for(int i=0;i<L;i++){
        cd s=cd(cos(PI/(1<<i)),sin(PI/(1<<i))),ss=conj(s);
        S[0][i][0]=S[1][i][0]=cd(1,0);
        for(int j=1;j<(1<<i);j++){
            S[0][i][j]=S[0][i][j-1]*s;
            S[1][i][j]=S[1][i][j-1]*ss;
        }
    }
}
void fft(cd *a,int p){
    for(int i=0;i<N;i++)
        if (i<rev[i])swap(a[i],a[rev[i]]);
    for(int i=2,ii=0;i<=N;i<<=1,ii++)
        for(int j=0;j<N;j+=i)
            for(int k=0;k<(i>>1);k++){
                cd x=a[j+k],y=a[j+k+(i>>1)]*S[p][ii][k];
                a[j+k]=x+y,a[j+k+(i>>1)]=x-y;
            }
    if (p){
        for(int i=0;i<N;i++)a[i]/=N;
    }
}
int main(){
    init();
    scanf("%d",&n);
    C=1e6,ans=1;
    for(int i=1;i<=n;i++){
        scanf("%d",&x);
        ans=(ll)ans*fac[x]%mod;
        x+=n-i;
        ans=(ll)ans*inv[x]%mod;
        a[x].real(a[x].real()+1);
        b[C-x].real(b[C-x].real()+1);
    }
    fft(a,0);
    fft(b,0);
    for(int i=0;i<N;i++)a[i]*=b[i];
    fft(a,1);
    for(int i=C+1;i<N;i++)ans=(ll)ans*qpow(i-C,floor(a[i].real()+0.5))%mod;
    printf("%d\n",ans);
    return 0;
}