hdu5197 DZY Loves Orzing(FFT+分治)
题目描述:一个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 ",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 }