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})!$分块打表吧,或者你想写快速阶乘也没人拦你

            分享图片
              1 #include<cstdio>
              2 #include<algorithm>
              3 #include<cstring>
              4 using namespace std;
              5 typedef long long lint;
              6 const lint mo=999948289,G=13;
              7 const int border=32768,N=40069,maxn=32000;
              8 lint fpow(lint a,lint p)
              9 {
             10     lint ret=1;
             11     while(p)
             12     {
             13         if(p&1ll) (ret*=a)%=mo;
             14         (a*=a)%=mo;
             15         p>>=1;
             16     }
             17     return ret;
             18 }
             19 
             20 lint fac[N],ifac[N];
             21 lint wg[N],iwg[N];
             22 lint ta[N],tb[N],tc[N];
             23 int inv[N];
             24 void ntt(lint *a,int len,int tp)
             25 {
             26     lint ilen=fpow(len,mo-2);
             27     for(int i=0;i<len;i++) if(i<inv[i]) swap(a[i],a[inv[i]]);
             28     for(int i=1;i<len;i<<=1)
             29     {
             30         lint w0=(~tp)?wg[i]:iwg[i];
             31         for(int j=0;j<len;j+=(i<<1))
             32         {
             33             lint w=1;
             34             for(int k=0;k<i;k++,(w*=w0)%=mo)
             35             {
             36                 lint w1=a[j+k],w2=w*a[j+k+i]%mo;
             37                 a[j+k]=(w1+w2)%mo,a[j+k+i]=(w1-w2+mo)%mo;
             38             }
             39         }
             40     }
             41     if(tp==-1) for(int i=0;i<len;i++) (a[i]*=ilen)%=mo;
             42 }
             43 lint dp[N];
             44 lint pa[20][2][N];
             45 void clr(int len)
             46 {
             47     memset(ta,0,len*8);
             48     memset(tb,0,len*8);
             49     memset(tc,0,len*8);
             50 }
             51 void work(int l,int r,int dep=0,int pos=0)
             52 {
             53     if(l==r){pa[dep][pos][1]=1,pa[dep][pos][0]=(lint)l-1;return;}
             54     int mm=l+r>>1;
             55     int n=r-l+2;
             56     for(int i=0;i<n;i++) pa[dep+1][0][i]=pa[dep+1][1][i]=0; 
             57     work(l,mm,dep+1,0),work(mm+1,r,dep+1,1);
             58     int len=1,pl=0;
             59     while(len<n) len<<=1,pl++;
             60     clr(len);
             61     for(int i=1;i<len;i++) inv[i]=(inv[i>>1]>>1)|((i&1)<<(pl-1));
             62     for(int i=0;i<n;i++) ta[i]=pa[dep+1][0][i],tb[i]=pa[dep+1][1][i];
             63     ntt(ta,len,1),ntt(tb,len,1);
             64     for(int i=0;i<len;i++) tc[i]=ta[i]*tb[i]%mo;
             65     ntt(tc,len,-1);
             66     for(int i=0;i<len;i++) pa[dep][pos][i]=tc[i];
             67 }
             68 
             69 
             70 lint calc(lint x);
             71 int n,a[100069];
             72 int main()
             73 {
             74     fac[0]=ifac[0]=1;
             75     for(int i=1;i<=maxn;i++) fac[i]=fac[i-1]*i%mo,ifac[i]=fpow(fac[i],mo-2);
             76     for(int i=1;i<border;i<<=1) wg[i]=fpow(G,(mo-1)/(i<<1)),iwg[i]=fpow(wg[i],mo-2);
             77     while(~scanf("%d",&n))
             78     {
             79         for(int i=1;i<=n;i++) scanf("%d",&a[i]);
             80         if(n>=maxn){puts("0");continue;}
             81         memset(pa[0][0],0,sizeof(pa[0][0]));
             82         work(1,n);
             83         lint ans=1;
             84         ans*=calc(n);
             85         (ans*=fpow(ifac[n],n))%=mo;
             86         for(int i=1;i<=n;i++) (ans*=pa[0][0][a[i]])%=mo;
             87         printf("%lld\n",ans);
             88     }
             89     return 0;
             90 }
             91 
             92 lint bfac[160]={
             93     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};
             94 lint calc(lint x)
             95 {
             96     lint b=x/200;
             97     lint ret=bfac[b];
             98     lint g=b*200+1,np=1ll*b*b*40000;
             99     lint pp=x*x;
            100     while(np<pp){np++;(ret*=np)%=mo;}
            101     return ret;
            102 }
            View Code
            相关文章
            相关标签/搜索
            每日一句
              每一个你不满意的现在,都有一个你没有努力的曾经。
            公众号推荐
               一个历史类的公众号,欢迎关注
            一两拨千金
            中彩堂zzyz,ccus时来棋牌游戏下载-时来棋牌游戏平台app最新版下载