hdu5197 DZY Loves Orzing(FFT+分治)

hdu5197 DZY Loves Orzing(FFT+分治)

hdu

题目描述：一个n*n的矩阵里填入1~n^2的数，要求每一排从前往后能看到a[i]个数(类似于身高阻挡视线那种)，求方案数。

思路：

考虑往一排里填入n个数。

经过简单推导发现正好有j个能被看到的方案数答案是$Sigma_{i=1}^{n}(x+(i-1))$的$x^{j}$项系数。

这个用分治FFT搞一搞就会变成$nlog^{2}n$的了。

之后再乘上一个$frac{(n^{2})!}{(n!)^{n}}$

woc这什么大数。。。

然而$n^2$超过模数的时候答案就是0了。

所以答案不为0的n被限制在了不到32000。

于是NTT也跑得快多了。

上面那个$(n^{2})!$分块打表吧，或者你想写快速阶乘也没人拦你

#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long lint;
const lint mo=999948289,G=13;
const int border=32768,N=40069,maxn=32000;
lint fpow(lint a,lint p)
{
    lint ret=1;
    while(p)
    {
        if(p&1ll) (ret*=a)%=mo;
        (a*=a)%=mo;
        p>>=1;
    }
    return ret;
}

lint fac[N],ifac[N];
lint wg[N],iwg[N];
lint ta[N],tb[N],tc[N];
int inv[N];
void ntt(lint *a,int len,int tp)
{
    lint ilen=fpow(len,mo-2);
    for(int i=0;i<len;i++) if(i<inv[i]) swap(a[i],a[inv[i]]);
    for(int i=1;i<len;i<<=1)
    {
        lint w0=(~tp)?wg[i]:iwg[i];
        for(int j=0;j<len;j+=(i<<1))
        {
            lint w=1;
            for(int k=0;k<i;k++,(w*=w0)%=mo)
            {
                lint w1=a[j+k],w2=w*a[j+k+i]%mo;
                a[j+k]=(w1+w2)%mo,a[j+k+i]=(w1-w2+mo)%mo;
            }
        }
    }
    if(tp==-1) for(int i=0;i<len;i++) (a[i]*=ilen)%=mo;
}
lint dp[N];
lint pa[20][2][N];
void clr(int len)
{
    memset(ta,0,len*8);
    memset(tb,0,len*8);
    memset(tc,0,len*8);
}
void work(int l,int r,int dep=0,int pos=0)
{
    if(l==r){pa[dep][pos][1]=1,pa[dep][pos][0]=(lint)l-1;return;}
    int mm=l+r>>1;
    int n=r-l+2;
    for(int i=0;i<n;i++) pa[dep+1][0][i]=pa[dep+1][1][i]=0; 
    work(l,mm,dep+1,0),work(mm+1,r,dep+1,1);
    int len=1,pl=0;
    while(len<n) len<<=1,pl++;
    clr(len);
    for(int i=1;i<len;i++) inv[i]=(inv[i>>1]>>1)|((i&1)<<(pl-1));
    for(int i=0;i<n;i++) ta[i]=pa[dep+1][0][i],tb[i]=pa[dep+1][1][i];
    ntt(ta,len,1),ntt(tb,len,1);
    for(int i=0;i<len;i++) tc[i]=ta[i]*tb[i]%mo;
    ntt(tc,len,-1);
    for(int i=0;i<len;i++) pa[dep][pos][i]=tc[i];
}


lint calc(lint x);
int n,a[100069];
int main()
{
    fac[0]=ifac[0]=1;
    for(int i=1;i<=maxn;i++) fac[i]=fac[i-1]*i%mo,ifac[i]=fpow(fac[i],mo-2);
    for(int i=1;i<border;i<<=1) wg[i]=fpow(G,(mo-1)/(i<<1)),iwg[i]=fpow(wg[i],mo-2);
    while(~scanf("%d",&n))
    {
        for(int i=1;i<=n;i++) scanf("%d",&a[i]);
        if(n>=maxn){puts("0");continue;}
        memset(pa[0][0],0,sizeof(pa[0][0]));
        work(1,n);
        lint ans=1;
        ans*=calc(n);
        (ans*=fpow(ifac[n],n))%=mo;
        for(int i=1;i<=n;i++) (ans*=pa[0][0][a[i]])%=mo;
        printf("%lld
",ans);
    }
    return 0;
}

lint bfac[160]={
    1,902199578,655434774,588857280,495770768,69882273,553982098,334078355,33146971,638472211,758245769,819694289,212913989,674505681,621807178,52420569,535922477,808220737,910946087,665159051,208303589,824486272,851100796,36810321,352031293,146240630,950654769,83962140,688846899,876526361,855642854,941799736,79240392,127370706,182824403,918730448,7023806,763878567,185845423,313214126,285253420,693669080,371386848,478395563,890609360,658191029,506004018,91639581,780064049,914814533,848366675,816348053,589401095,509135319,5446319,31619815,253732202,31857169,381860443,165388954,340902365,960303088,888954496,990221261,563929977,876772121,287079597,609490658,652825564,188993794,549908577,448482523,951503233,573959686,967072042,138776107,103551474,507659875,978747064,744346128,6336599,284817804,56458444,949826314,426241014,131445497,844320792,317915054,729308605,623307135,415416377,694294158,132884549,468178276,709909378,215890898,517268673,756318464,591277985,591506763,315367468,727640178,83321953,520348269,358779796,726705701,333860770,574957987,792831525,568945082,894644404,585412875,460795027,380631820,408052670,276231384,641853065,627050439,424695359,849416232,257485155,78105895,255310610,546026921,968605310,537749882,484855581,889099301,930219111,978770851,22465690,316568861,457634470,176979627,346624108,121259373,302342099,722399773,64910873,521979194,576421397,629834183,904179107,507024536,317861768,202275207,572048958,636515142,125373754,68190273,464124136,875197733,313248262,908784630,959938689,532154720,996633237,161891031,288552800,0};
lint calc(lint x)
{
    lint b=x/200;
    lint ret=bfac[b];
    lint g=b*200+1,np=1ll*b*b*40000;
    lint pp=x*x;
    while(np<pp){np++;(ret*=np)%=mo;}
    return ret;
}