FFT小結！

零，說在前面

（轉載請注明原文地址：http://www.cnblogs.com/LadyLex/p/8335420.html?）

……令人毒瘤的FFT

其實FFT是卷積這個龐大學術課題的一個小應用啦.....

WC的時候聽卷積定理真是上天的體驗

那么話不多說我們開始FFT的總結.

板子我們就不放了!這個是很基礎的,網上也有講解的很詳細的學習資料

大家請善用搜索引擎啦......

一，FFT與NTT的簡單例題以及基礎思想&技巧

首先要提到的就是反轉思想了。

可能有的題目形式并不是卷積，但是我們通過改變枚舉順序，或者反轉數組就可以得到卷積的形式

這樣的一個典型題目就是bzoj2194了，非常基礎的題目……

這里就不推式子了

?

然后是關于字符串匹配的問題……說一種最簡單的做法吧，我們把一個串反過來，然后同種字符設為1，其余的為0，然后做卷積

這樣就可以拿到哪些位置匹配了，然后把每種字符都做一遍就可以得出2個字符串有多長是一樣的……

的確有的題問這個的……

?

以及FFT還可以和斯特林數結合起來……斯特林數大概需要掌握的是容斥的求法以及用fft優化的求法,以及一個重要的轉化式子

$$n^{k}= \sum _{j=0} ^{k}? S(k,j) * j!* C(n,j) $$

除了數本身,它的實際意義也對于推導有幫助.

學習的時候參考了學長Aglove的筆記........https://www.cnblogs.com/joyouth/p/5600541.html

?

此外,集合計數類問題常用的2個思路是"枚舉最后一個集合的大小"或者"枚舉集合的個數"

在很多題目中都有應用!

?

還有一個技巧,是在出現組合數的時候兩邊同除階乘,可以轉化為多項式求逆的做法

很多需要用CDQ+FFT的題目都可以轉化成多項式求逆

比如bzoj3456,bzoj4555這些都可以使用

二，FFT的更多應用

其實我和Cooook寫了一個課件的包括了底下大部分東西(其實只有1~8的內容)……啊我好懶啊我不想再打一遍了

如果有想要課件的同學跟我說下啊……我把密碼放在評論區……

鏈接

一，多項式求逆

二，多項式開根

三，多項式除法&取模

四，多項式取$ln$

五，多項式取$exp$

六，多項式快速插值

七，多項式多點求值

八，指數型母函數

九，任意模數FFT

這個我還沒有看……現在正在學ing，請等我補坑……

十，與卷積定理的結合

這個在WC2018上聽LCA爺爺講了講……聽傻了……

現在我還是不會的…估計一段時間內我也不太可能會…留坑待填

十一，特征多項式與矩陣

我們把一個矩陣A的特征值定義為滿足$Ax=\lambda x$的$\lambda?$

然后特征向量就是上面那個x

特征多項式$g(x)$定義為$|A-xI|$,其中I為單位矩陣

那么這個行列式里面應該是帶著一堆x的,這也是它是多項式的原因.....

這個多項式的根就是這個矩陣的特征值

那么這玩意有什么用?

一個很重要的定理是叫做凱萊哈密頓定理,內容是$g(A)=0$,0代表0矩陣

至于矩陣怎么帶到多項式里面....請把矩陣看成變量

這樣的話我們就可以快速進行矩陣快速冪了.下面說下個人的理解

我們可以通過構造把任意矩陣冪次表示為$A^{n}=k(A)*f(A)+l(A)$的形式,其中k為另外某個多項式,$f(x)$為某多項式,$l$為某矩陣

那么當$f(x)$取$g(x)$的時候,$k*f(x)=k*g(A)=0$

這樣的話,就有$A^{n}=l$

那么我們先求$x^{n} \mod g(x)$的值,設其為$h(x)$,那么$h(A)=l$,我們就得到了$A^{n}$

然后這東西一個很大的用途就是處理常系數線性遞推的問題

設有遞推式:

$$f(n)=\sum _{i=1} ^{k} f_{n-i} a_{i}$$

我們先設系數$a$個數為k,需要求$f$的第n項,

這本來是矩陣乘法的活,但是如果矩陣太大的話時間復雜度就爆炸了......

這時候我們先求出上面提到的$h(x)=x^{n-k} \mod g(x)$,然后我們考慮接下來要干什么

假設我們已經知道了要求的遞推式的前$2k$項,

那么最后第n項應該是等于一開始的一橫條矩陣(設為T)里面的每一項(即$f(k)$到$f(1)$),然后乘$A^{n-k}$矩陣的對應位置

因為我們把矩陣帶入了$h(x)$,即$h(A)=A^{n-k}$

那么我們可以分別對$h(x)$的每一項計算它對答案的貢獻

那么其實假如我們在統計第j項的時候,枚舉到了$f(i)$,那么它如果乘了$A^{j}$對應位置的系數,那么就等價于變成了$f(i+j)$項

然后我們再乘對應的$h(x)$中系數即可更新$f(n-k+i)$的值

暴力進行這個過程,求單個$f(n)$是$O(k)$的,如果求多個(比如最后k個)$f(n)$,我們還可以使用卷積進行優化

最后一個問題,$g(x)$怎么求?我們對于常系數線性遞推有一個結論是

$$g(x)=(-1)^{k} * ( x^{k}-\sum _{i=1} ^{k} a_{i} x^{k-i} )$$

其中$a_{i}$代表對應的第i個系數,這個結論的具體證明和$A$矩陣的特殊形態有關,可是我還不會證明.....

先當作結論記下來記下來

我們的紅太陽Troywar給出了證明，我丟個鏈接就跑

三，FFT/NTT好題選講

后面的題我在不斷做不斷更新中……

一，UOJ#50 鏈式反應

這個題我們發現其實就是一個類似于二叉樹計數的題目。

首先我們發現那個標號的限制完全沒有用……

不難想到一個$O(n^{3})$的dp方程，即枚舉兩個中子打到的子樹大小，然后再看剩下的大小是否符合破壞死光的要求，再乘組合數。這是一個40pts做法。

然后，這個dp出來之后，我們發現兒子的貢獻就是一個卷積的形式啊！

然后把合并寫成指數型母函數的卷積，這樣的復雜度是$O(n^{2}log_{n})$,可以拿到60pts

然后……你會很快的發現，如果我們把死光也看成多項式，那么就是3個卷起來，那么復雜度就成了$O(nlog_{n}^{2})$，用CDQ+FFT可以得到90pts～100pts

那么我個人就是打的這種方法！4s打2e5的CDQ+FFT是沒有問題的！

最有意思的就是那個多項式的維護了。由于我們的式子是這樣的：

$$F'(x)=C(x)*F^{2}(x)+1$$

所以我們要維護一個平方的多項式，代表表示1～區間中點mid平方的結果

每次我們用平方數組更新f，然后用新計算了mid+1~r，我們就把它更新到原來那個平方里面去

怎么更新？提示：$(x+a)^{2}=x^{2}+2ax+a^{2}$

當然，我們還有一種解微分方程的做法，但是鑒于我的數學水平我并不會23333。

有興趣的同學可以看Vfk的題解。

給個代碼：

1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 #include <cstdlib> 5 using namespace std; 6 #define N 400010 7 #define L 1100000 8 #define mod 998244353 9 #define phi 998244352 10 #define g 3 11 #define LL long long 12 inline int max(int a,int b){return a>b?a:b;} 13 inline int min(int a,int b){return a<b?a:b;} 14 char s[N]; 15 int f[N],pf[N],n,logi[N],litlen=1,len=1; 16 int bin[25],poww[L],rev[L],fac[L],inv[L],invf[L]; 17 inline int quick_mod(int di,int mi) 18 { 19 int ret=1; 20 for(;mi;mi>>=1,di=(LL)di*di%mod) 21 if(mi&1)ret=(LL)ret*di%mod; 22 return ret; 23 } 24 inline void upd(int &a,int b){a+=b;if(a>=mod)a-=mod;} 25 inline int sub(int a,int b){a-=b;if(a<0)a+=mod;return a;} 26 inline int dft(int *a,int opt,int le) 27 { 28 register int i,j,d,l,wn,w,tmp; 29 for(i=0;i<le;++i)if(i<rev[i])swap(a[i],a[rev[i]]); 30 for(d=2;d<=le;d<<=1) 31 for(wn=(opt==1)?poww[len/d]:poww[len-len/d],l=d>>1,i=0;i<le;i+=d) 32 for(w=1,j=0;j<l;++j,w=(LL)w*wn%mod) 33 tmp=(LL)a[i+j+l]*w%mod,a[i+j+l]=(a[i+j]-tmp+mod)%mod,a[i+j]=(a[i+j]+tmp)%mod; 34 // tmp=(LL)a[i+j+l]*w%mod,a[i+j+l]=sub(a[i+j],tmp),upd(a[i+j],tmp); 35 if(opt==-1) 36 for(i=0;i<le;++i)a[i]=(LL)a[i]*inv[le]%mod; 37 } 38 int ret[L],C[L],A[L],B[L]; 39 inline void CDQ(int l,int r) 40 { 41 if(l==r) 42 if(l<n){printf("%d\n",(LL)f[l]*fac[l-1]%mod);return;} 43 else {printf("%d\n",(LL)f[l]*fac[l-1]%mod);exit (0);} 44 register int le=r-l+1,i,mi=(l+r)>>1; 45 CDQ(l,mi); 46 for(i=0;i<le;++i) 47 if(i&1)rev[i]=(rev[i>>1]>>1)|(le>>1); 48 else rev[i]=(rev[i>>1])>>1; 49 memset(A,0,le<<2),memset(B,0,le<<2); 50 for(i=l;i<=mi;++i)A[i-l]=pf[i]; 51 for(i=0;i<r-l;++i)B[i]=C[i]; 52 dft(A,1,le),dft(B,1,le); 53 for(i=0;i<le;++i)ret[i]=(LL)A[i]*B[i]%mod; 54 dft(ret,-1,le); 55 for(i=mi+1;i<=r;++i)upd(f[i],ret[i-l-1]); 56 memset(A,0,le<<2),memset(B,0,le<<2); 57 for(i=l;i<=mi;++i)A[i-l]=(LL)f[i]*inv[i]%mod; 58 for(i=0;i<=r-l;++i) 59 if(i<l)B[i]=(LL)f[i]*inv[i]*2%mod; 60 else if(i<=mi)B[i]=(LL)f[i]*inv[i]%mod; 61 dft(A,1,le),dft(B,1,le); 62 for(i=0;i<le;++i)ret[i]=(LL)A[i]*B[i]%mod; 63 dft(ret,-1,le); 64 for(i=mi+1;i<=r;++i)upd(pf[i],ret[i-l]); 65 CDQ(mi+1,r); 66 } 67 int main() 68 { 69 register int i,j; 70 for(bin[0]=i=1;i<=20;++i)bin[i]=bin[i-1]<<1; 71 scanf("%d%s",&n,s); 72 while(len<(n<<1))len<<=1; 73 poww[0]=1,poww[1]=quick_mod(g,phi/len); 74 for(i=2;i<=len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod; 75 for(fac[0]=fac[1]=1,i=2;i<=len;++i)fac[i]=(LL)fac[i-1]*i%mod; 76 for(inv[0]=inv[1]=1,i=2;i<=len;++i)inv[i]=(LL)inv[mod%i]*(mod-mod/i)%mod; 77 for(invf[0]=invf[1]=1,i=2;i<=len;++i)invf[i]=(LL)invf[i-1]*inv[i]%mod; 78 for(i=0;i<n;++i) 79 if(s[i]=='1')C[i]=(LL)invf[i]*inv[2]%mod; 80 f[1]=1,CDQ(1,len>>1); 81 } UOJ#50

二，UOJ#182 a^-1+b problem

這題我在課件里面講了！

然后寫的自己還算滿意吧……再簡單說一下

總之你可以發現每個數再若干次操作之后可以寫成$\frac {Ax+C}{Bx+D}$的形式，x代表每個數的初始值。然后……你把它寫丑一點，可以變成形如$e+f*\frac{1}{x+g}$的形式……但是具體內容請自己推導吧……

那么總和就是$n*e+f*\sum\frac{1}{x+g}$，問題轉化為維護$\sum\frac{1}{x_{i}+g}$

如果我們把g看成變量的話

那么這個東西，通分可以得到等價于$\frac{? \sum _{j} \prod _{i!=j}? ?{A_{i}+x}? ?}? {\prod _{j} (A_{j}+x)}$

注意到分母里的式子（設為$f$）可以倍增加FFT求出來

分子的式子（設為$g$），經過我們的推導，是$f$的導數

你問我怎么推導？

第一種方法，你可以打表

第二種方法，觀察特殊點……我當時只看到了$[x^{n-1}]g(x)=n$

第三種方法，我們推他一下……

啊……我推了半小時推不出來

如果有會的同學請在評論區教教我

然后呢？我們可以套用多項式多點求值的模板，求出兩個函數每一個點的點值，然后計算答案就行了

然后這題讓我知道了指針的常數很小23333

1 #include <cstdio> 2 #include <cstring> 3 #include <ctime> 4 #include <algorithm> 5 using namespace std; 6 #define mod 998244353 7 #define LL long long 8 #define RG register 9 #define N 100010 10 #define M 60010 11 char cB[1<<15],*S=cB,*T=cB; 12 #define getc (S==T&&(T=(S=cB)+fread(cB,1,1<<15,stdin),S==T)?0:*S++) 13 inline int read() 14 { 15 int x=0;RG char c=getc; 16 while(c<'0'|c>'9')c=getc; 17 while(c>='0'&c<='9')x=10*x+(c^48),c=getc; 18 return x; 19 } 20 inline int quick_mod(int di,int mi) 21 { 22 int ret=1; 23 for(;mi;mi>>=1,di=(LL)di*di%mod) 24 if(mi&1)ret=(LL)ret*di%mod; 25 return ret; 26 } 27 #define L (1<<18)+10 28 #define LM (1<<16)+10 29 int logi[L],inv[L],bin[25],rev[L],poww[L],len=1,ro=1; 30 inline void dft(int *a,int n,int opt) 31 { 32 RG int i,j,d=logi[len]-logi[n],wn,l,tmp,*x,*y,*w; 33 for(i=0;i<n;++i)if(i<(rev[i]>>d))swap(a[i],a[rev[i]>>d]); 34 for(d=2;d<=n;d<<=1) 35 for(i=0,l=d>>1,wn=(opt==1)?(len/d):(-len/d);i<n;i+=d) 36 for(w=poww+((opt==1)?0:len),j=0,x=a+i+j,y=x+l;j<l;++j,w+=wn,++x,++y) 37 tmp=(LL)(*w)*(*y)%mod,*y=(*x-tmp+mod)%mod,*x=(*x+tmp)%mod; 38 if(opt==-1)for(i=0;i<n;++i)a[i]=(LL)a[i]*inv[n]%mod; 39 } 40 int n,m,cnt,sum,x[N]; 41 int g[L<<1],h[L<<1]; 42 int sta[LM],ans[M]; 43 struct node{int e,f,g,id;}s[M]; 44 inline void ins(int a,int b,int c,int d,int id) 45 { 46 if(b==0) 47 {ans[id]=(((LL)a*sum)+((LL)c*n))%mod*(LL)quick_mod(d,mod-2)%mod;return;} 48 s[++cnt].f=((LL)b*c-(LL)a*d)%mod, 49 b=quick_mod(b,mod-2); 50 s[cnt].e=(LL)a*b%mod,s[cnt].g=(LL)d*b%mod, 51 b=(LL)b*b%mod,s[cnt].f=(LL)s[cnt].f*b%mod; 52 if(s[cnt].f<0)s[cnt].f+=mod; 53 s[cnt].id=id; 54 } 55 inline int solve0(int *a,int l,int r) 56 { 57 RG int ra=r-l+2; 58 if(ra<=128) 59 { 60 memset(a,0,ra<<2),a[0]=1; 61 RG int i,j,v; 62 for(i=l;i<=r;++i) 63 for(v=x[i],j=i-l;~j;--j) 64 a[j+1]=(a[j+1]+a[j])%mod, 65 a[j]=(LL)a[j]*v%mod; 66 return ra; 67 } 68 RG int i,mi=l+r>>1,la=1; 69 while(la<ra)la<<=1; 70 RG int *f1=a,r1=solve0(f1,l,mi); 71 RG int *f2=a+r1,r2=solve0(f2,mi+1,r); 72 static int tmp1[L],tmp2[L]; 73 memcpy(tmp1,f1,r1<<2),memset(tmp1+r1,0,(la-r1)<<2),dft(tmp1,la,1); 74 memcpy(tmp2,f2,r2<<2),memset(tmp2+r2,0,(la-r2)<<2),dft(tmp2,la,1); 75 for(i=0;i<la;++i)a[i]=(LL)tmp1[i]*tmp2[i]%mod; 76 dft(a,la,-1); 77 return ra; 78 } 79 int *p[L]; 80 inline int solve1(int id,int l,int r) 81 { 82 static int mem[LM<<4],*head=mem; 83 RG int ra=r-l+2; 84 if(ra<=128) 85 { 86 RG int i,j,v,*f=p[id]=head;head+=ra; 87 memset(f,0,ra<<2),f[0]=1; 88 for(i=l;i<=r;++i) 89 for(v=(mod-sta[i])%mod,j=i-l;~j;--j) 90 f[j+1]=(f[j+1]+f[j])%mod,f[j]=(LL)f[j]*v%mod; 91 return ra; 92 } 93 RG int mi=l+r>>1,la=1; 94 while(la<ra)la<<=1; 95 RG int r1=solve1(id<<1,l,mi),*f1=p[id<<1]; 96 RG int r2=solve1(id<<1|1,mi+1,r),*f2=p[id<<1|1]; 97 static int tmp1[LM],tmp2[LM]; 98 memcpy(tmp1,f1,r1<<2),memset(tmp1+r1,0,(la-r1)<<2),dft(tmp1,la,1); 99 memcpy(tmp2,f2,r2<<2),memset(tmp2+r2,0,(la-r2)<<2),dft(tmp2,la,1); 100 RG int i,*f=p[id]=head;head+=ra; 101 for(i=0;i<la;++i)f[i]=(LL)tmp1[i]*tmp2[i]%mod; 102 dft(f,la,-1);return ra; 103 } 104 inline int get_inv(int *a,int *ret,int ra) 105 { 106 if(ra==1){ret[0]=quick_mod(a[0],mod-2);return 1;} 107 static int tmp[L]; 108 RG int r1=get_inv(a,ret,(ra+1)>>1),i,la=1; 109 while(la<(ra<<1))la<<=1; 110 memcpy(tmp,a,ra<<2),memset(tmp+ra,0,(la-ra)<<2); 111 memset(ret+r1,0,(la-r1)<<2); 112 dft(tmp,la,1),dft(ret,la,1); 113 for(i=0;i<la;++i) 114 { 115 ret[i]=((LL)ret[i]*(2ll-((LL)ret[i]*tmp[i]%mod)))%mod; 116 if(ret[i]<0)ret[i]+=mod; 117 } 118 dft(ret,la,-1);return ra; 119 } 120 inline void rev_copy(int *st,int *to,int ra) 121 {for(RG int i=0;i<ra;++i)to[i]=st[ra-i-1];} 122 inline void reverse(int *st,int ra) 123 {for(RG int t,i=0,j=ra-1;i<j;++i,--j)t=st[i],st[i]=st[j],st[j]=t;} 124 inline int get_mod(int *a,int ra,int *p,int rp,int *ret) 125 { 126 while(ra&&!a[ra-1])--ra; 127 while(rp&&!p[rp-1])--rp; 128 if(ra<rp){memcpy(ret,a,ra<<2),memset(ret+ra,0,(rp-ra)<<2);return rp;} 129 static int tmp1[L],tmp2[L]; 130 RG int i,j,re=ra-rp+1,la=1; 131 while(la<(re<<1))la<<=1; 132 memset(tmp1,0,la<<2),rev_copy(p,tmp1,rp),get_inv(tmp1,tmp2,re); 133 memset(tmp2+re,0,(la-re)<<2),dft(tmp2,la,1); 134 rev_copy(a,tmp1,ra),memset(tmp1+re,0,(la-re)<<2),dft(tmp1,la,1); 135 for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod; 136 dft(tmp1,la,-1); 137 la=1;while(la<ra)la<<=1; 138 reverse(tmp1,re),memset(tmp1+re,0,(la-re)<<2),dft(tmp1,la,1); 139 memcpy(tmp2,p,rp<<2),memset(tmp2+rp,0,(la-rp)<<2),dft(tmp2,la,1); 140 for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod; 141 dft(tmp1,la,-1); 142 for(i=0;i<rp;++i)ret[i]=(a[i]-tmp1[i]+mod)%mod; 143 memset(ret+rp,0,(la-rp)<<2); 144 while(rp&&!ret[rp-1])--rp; 145 return rp; 146 } 147 int val0[LM],val[LM]; 148 inline void solve2(int id,int *a,int *b,int l,int r,int deg) 149 { 150 RG int ra=r-l+2; 151 if(deg<=128) 152 { 153 RG int i,j,v,u,F,G; 154 for(i=l;i<=r;val0[i]=F,val[i]=G,++i) 155 for(u=sta[i],v=1,F=G=0,j=0;j<=deg;++j,v=(LL)v*u%mod) 156 F=(F+(LL)v*a[j])%mod,G=(G+(LL)v*b[j])%mod; 157 return; 158 } 159 RG int i,mi=l+r>>1; 160 int r1=get_mod(a,deg,p[id],ra,a+deg);a+=deg; 161 int r2=get_mod(b,deg,p[id],ra,b+deg);b+=deg; 162 ra=min(r1,r2); 163 solve2(id<<1,a,b,l,mi,ra),solve2(id<<1|1,a,b,mi+1,r,ra); 164 } 165 int main() 166 { 167 RG int i,j,opt,v,A=1,B=0,C=0,D=1; 168 n=read(),m=read(); 169 while(len<(n<<1))len<<=1; 170 for(bin[0]=i=1;i<=20;++i)bin[i]=bin[i-1]<<1; 171 for(i=0;i<=20;++i) 172 inv[bin[i]]=quick_mod(bin[i],mod-2),logi[bin[i]]=i; 173 for(int i=0;i<len;++i) 174 if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1); 175 else rev[i]=(rev[i>>1]>>1); 176 poww[0]=1,poww[1]=quick_mod(3,(mod-1)/len); 177 for(i=2;i<=len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod; 178 for(sum=0,i=1;i<=n;++i)x[i]=read(),sum=(sum+x[i])%mod; 179 for(i=1;i<=m;++i) 180 if(read()==1)v=read(),A=(A+(LL)v*B)%mod,C=(C+(LL)v*D)%mod,ins(A,B,C,D,i); 181 else swap(A,B),swap(C,D),ins(A,B,C,D,i); 182 if(cnt) 183 { 184 for(i=1;i<=cnt;++i)sta[i]=s[i].g; 185 sort(sta+1,sta+cnt+1),opt=unique(sta+1,sta+cnt+1)-sta-1; 186 for(i=1;i<=cnt;++i) 187 s[i].g=lower_bound(sta+1,sta+opt+1,s[i].g)-sta; 188 solve0(g,1,n); 189 for(h[n]=0,i=1;i<=n;++i)h[i-1]=(LL)g[i]*i%mod; 190 solve1(1,1,opt),solve2(1,g,h,1,opt,n+1); 191 for(i=1;i<=cnt;++i) 192 ans[s[i].id]=( (LL)s[i].e*n+(LL)s[i].f*val[s[i].g]%mod*quick_mod(val0[s[i].g],mod-2) )%mod; 193 } 194 for(i=1;i<=m;++i)printf("%d\n",ans[i]); 195 } UOJ#182

三，UOJ#86

myy的題真是太毒辣！

首先我們要知道c++11有一個__int128，可以當一個假高精度用

觀察暴力……我們當然可以寫一個盧卡斯暴力判斷

但是我們看，把盧卡斯的暴力寫出來，應該是

$$C(x,n)=C(x\%p,n\%p)*C(x/p,n/p)$$

也就是說可以拆成一堆C的$\prod$

這樣它和p進制數位DP很像……那么我們把數轉成p進制

定義數組$f[i][j]$是前i位，那一堆C乘起來是j的方案數，然后枚舉這一位填什么，計算C乘起來轉移

然后我們看，它類似一個假的卷積

我們通過原根使得它變成加法，再用FFT優化，就可以做到$O(plog_{p}^{n}log_{2}^{p})$的復雜度了

其實我覺得最難想的是變成數位DP這個部分……如果想到數位DP后面就沒那么難了……

丟代碼：

1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 using namespace std; 5 #define LLL __int128 6 #define LL long long 7 #define RG register 8 #define PMAX 30010 9 #define mod 998244353 10 inline int quick_mod(int di,int mi) 11 { 12 RG int ret=1; 13 for(;mi;mi>>=1,di=(LL)di*di%mod) 14 if(mi&1)ret=(LL)ret*di%mod; 15 return ret; 16 } 17 #define L (1<<16)+10 18 int len=1,rev[L],poww[L]; 19 inline void dft(int *a,int opt) 20 { 21 RG int i,j,d,l,*w,wn,tmp,*x,*y; 22 for(i=0;i<len;++i)if(i<rev[i])swap(a[i],a[rev[i]]); 23 for(d=2;d<=len;d<<=1) 24 for(i=0,l=(d>>1),wn=(opt==1?(len/d):(-len/d));i<len;i+=d) 25 for(j=0,w=(poww+(opt==1?0:len)),x=a+i,y=x+l;j<l;++j,++x,++y,w+=wn) 26 tmp=(LL)(*w)*(*y)%mod,*y=(*x+mod-tmp)%mod,*x=(*x+tmp)%mod; 27 if(opt==-1) 28 for(tmp=quick_mod(len,mod-2),i=0;i<len;++i)a[i]=(LL)a[i]*tmp%mod; 29 } 30 int p,phi,root,logv[PMAX],expg[PMAX]; 31 int fac[PMAX],invf[PMAX]; 32 inline int C(int a,int b) 33 {return (a<b)?0:(fac[a]*invf[b]%p*invf[a-b]%p);} 34 inline int quick_p(int di,int mi) 35 { 36 RG int ret=1; 37 for(;mi;mi>>=1,di=di*di%p) 38 if(mi&1)ret=ret*di%p; 39 return ret; 40 } 41 inline void getroot() 42 { 43 RG int i,j,v; 44 for(i=1;i<p;++i) 45 { 46 memset(logv,-1,sizeof(logv)); 47 logv[1]=0,expg[0]=1; 48 for(v=i,j=1;j<phi;v=v*i%p,++j) 49 {if(logv[v]>=0)break;expg[j]=v,logv[v]=j;} 50 if(j==phi){root=i;return;} 51 } 52 } 53 char B[1<<15],*S=B,*T=B; 54 #define getc (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++) 55 template <typename _T> 56 inline _T read() 57 { 58 RG _T x=0;RG char c=getc; 59 while(c<'0'|c>'9')c=getc; 60 while(c>='0'&&c<='9')x=10*x+(c^48),c=getc; 61 return x; 62 } 63 int lenn,bn[300]; 64 int f[300][L],g[300][L]; 65 int ans[PMAX],tmp[L],tmp1[L],tmp2[L]; 66 inline void calc(LLL n,int opt) 67 { 68 if(n<=0)return; 69 RG int i,j,v,k,bin[300],to,bit=0; 70 memset(bin,0,sizeof(bin)); 71 while(n)bin[++bit]=n%p,n/=p; 72 RG int lim=max(bit,lenn); 73 memset(f,0,sizeof(f)),f[0][0]=1; 74 if(opt==1) 75 { 76 memset(g,0,sizeof(g)),g[0][0]=1; 77 for(i=1;i<=lim;++i) 78 { 79 memset(tmp1,0,sizeof(tmp1)); 80 for(j=bn[i];j<p;++j) 81 {v=logv[C(j,bn[i])];if(v!=-1)++tmp1[v];} 82 memcpy(tmp2,g[i-1],phi<<2); 83 memset(tmp2+phi,0,(len-phi)<<2); 84 dft(tmp1,1),dft(tmp2,1); 85 for(j=0;j<len;++j)tmp[j]=(LL)tmp1[j]*tmp2[j]%mod; 86 dft(tmp,-1); 87 for(j=0;j<phi;++j) 88 g[i][j]=(tmp[j]+tmp[j+phi])%mod; 89 } 90 } 91 for(i=1;i<=lim;++i) 92 { 93 memset(tmp1,0,sizeof(tmp1)); 94 for(j=bn[i];j<bin[i];++j) 95 { 96 v=logv[C(j,bn[i])]; 97 if(v!=-1)++tmp1[v]; 98 } 99 memcpy(tmp2,g[i-1],phi<<2); 100 memset(tmp2+phi,0,(len-phi)<<2); 101 dft(tmp1,1);dft(tmp2,1); 102 for(j=0;j<len;++j)tmp[j]=(LL)tmp1[j]*tmp2[j]%mod; 103 dft(tmp,-1); 104 for(j=0;j<phi;++j)f[i][j]=(tmp[j]+tmp[j+phi])%mod; 105 v=logv[C(bin[i],bn[i])]; 106 if(v!=-1)for(k=0;k<phi;++k) 107 to=(v+k)%phi,f[i][to]=(f[i][to]+f[i-1][k])%mod; 108 } 109 if(opt==1)for(i=0;i<phi;++i)ans[expg[i]]=(ans[expg[i]]+f[lim][i])%mod; 110 else for(i=0;i<phi;++i)ans[expg[i]]=(ans[expg[i]]+mod-f[lim][i])%mod; 111 } 112 int main() 113 { 114 RG int i,j,sum; 115 p=read<int>(),phi=p-1,getroot(); 116 while(len<=(p<<1))len<<=1; 117 for(i=0;i<len;++i) 118 if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1); 119 else rev[i]=rev[i>>1]>>1; 120 poww[0]=poww[len]=1,poww[1]=quick_mod(3,(mod-1)/len); 121 for(i=1;i<len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod; 122 for(fac[0]=fac[1]=1,i=2;i<p;++i)fac[i]=fac[i-1]*i%p; 123 invf[p-1]=quick_p(fac[p-1],p-2); 124 for(i=p-1;i;--i)invf[i-1]=invf[i]*i%p; 125 LLL n,l,r; 126 n=read<LLL>(),l=read<LLL>(),r=read<LLL>(); 127 while(n)bn[++lenn]=n%p,n/=p; 128 calc(r,1),calc(l-1,-1); 129 for(sum=(r-l+1)%mod,i=1;i<p;++i)sum=(sum-ans[i]+mod)%mod; 130 for(printf("%d\n",sum),i=1;i<p;++i)printf("%d\n",ans[i]); 131 } UOJ#86

四，UOJ#23

?留坑待填

五，UOJ#316 NOI2017泳池

去年NOI的試題.......和Cooook研究了一晚上+一上午才研究明白

多項式取模的板子還是打不熟......

首先我們要有一個DP的式子......

這個東西由于你不知道那個矩形在哪飄著所以你很難想

但是一個模糊的思路肯定是有的,就是按列轉移

我一開始看這題發現隨著寬度變大高度會縮小....然后最后證明這其實是預處理時間復雜度的保障

我們先用類似容斥的想法,設$s_{k}$為最大子矩形小于等于k的概率,那么答案等于$s_{k}-s_{k-1}$

然后我們考慮剛才的按列轉移,如果我能算出來對于一個$j*1001$的矩形滿足合法要求的概率$legal(j)$,我就枚舉這個$j$,然后在它旁邊填上一個x就行了

大概畫圖長這樣

那么你可以定義$f(i)$為前i列合法,并且最后一列底下是x的概率

然后統計答案的時候我就求出最后$k+1$個f值,然后在他們后面補充一個第一行沒x的合法區間(因為合法解不一定最后一列第一行是x)

那么我們可以寫出式子

$$f(i)= \sum _{j=1} ^{ min(i,k) } f(i-j-1)*(1-p)*legal(j) $$

$$ans=\sum _{i=1} ^{k}? f(n-i) legal(i)$$

當然你要先判掉$n<=k$

這樣我們問題就變成了處理長度為j的區間合法的概率

這東西怎么求呢

我們利用和f的求法類似的求法來預處理

設$g(i,j)$為j列的矩形,前i行全部安全,第$i$行以上未知,并且最后合法的概率

再設$h(i,j)$為j列的矩形,前i行全部安全,第$i$行以上未知,第$j$列的第$i+1$行確定為$x$,并且最后合法的概率

然后我們就用$g$和$h$來回轉移

具體的做法,是從高處往低處轉移

一開始我們可以直接算出$g(k,1)$,然后我們由第$i+1$行更新第$i$行,枚舉第$i+1$行哪個位置是$x$使得我們無法確保第$i+1$行

大概是這樣的:

$$g(i,j)=\sum _{k=1} ^ {K/i}? ?h(i,k)*g(i+1,j-k)$$

然后h的轉移和它類似

$$h(i,j)=\sum _{k=1} ^ {K/i}? ?h(i,k)*g(i+1,j-k-1)*(1-p)*p^{i}$$

這樣,最后求出的$g(1,j)$就是我們剛才定義的$legal(j)$

最后再用上面介紹的矩陣特征多項式優化一下,就可以做到$O(k^{2}+klog_{k}log_{n})$了

一開始我們要先求f的前k項,本題的k很小,所以我選擇了$O(k^2)$暴力,其實我們也可以用CDQ分治或者多項式求逆的解法,復雜度更加優秀

同理,求出最后$k+1$項也可以用$FFT$卷積優化,我也用了暴力

其實這題都不用$FFT$和多項式取模,暴力卷積和取模即可

其實我一開始看暴力取模看不懂,后來發現....其實是類似豎式除法的過程

這樣這道666的題目就被我們解決了

后記:這題讓Cooook有一種加強數據做毒瘤題的沖動.....我及時制止了他23333

那么我們扔代碼

1 #include <cstdio> 2 #include <algorithm> 3 #include <cstring> 4 using namespace std; 5 #define mod 998244353 6 #define K 1010 7 #define L 4096 8 #define RG register 9 #define LL long long 10 inline int quick_mod(int di,int mi) 11 { 12 int ret=1; 13 for(;mi;mi>>=1,di=(LL)di*di%mod) 14 if(mi&1)ret=(LL)ret*di%mod; 15 return ret; 16 } 17 int n,k,p,q; 18 int g[K][K],h[K][K],pp[K],f[K<<1],ff[K<<1]; 19 int tmp[L],tmp1[L],tmp2[L],tmp3[L],c[L],d[L],e[L]; 20 inline int min(int a,int b){return a<b?a:b;} 21 int poww[L],logi[L],rev[L],len=1; 22 inline void dft(int *a,int ra,int opt) 23 { 24 register int i,j,l,d=logi[len]-logi[ra],wn,tmp,*w,*x,*y; 25 for(i=0;i<ra;++i)if(i<(rev[i]>>d))swap(a[i],a[rev[i]>>d]); 26 for(d=2;d<=ra;d<<=1) 27 for(wn=((opt==1)?(len/d):(-len/d)),i=0,l=(d>>1);i<ra;i+=d) 28 for(w=poww+(opt==1?0:len),j=0,x=a+i,y=x+l;j<l;++j,++x,++y,w+=wn) 29 tmp=(LL)(*w)*(*y)%mod,*y=(*x-tmp+mod)%mod,*x=(*x+tmp)%mod; 30 if(opt==-1) 31 for(tmp=quick_mod(ra,mod-2),i=0;i<ra;++i)a[i]=(LL)a[i]*tmp%mod; 32 } 33 inline int get_inv(int *a,int *ret,int ra) 34 { 35 if(ra==1){ret[0]=quick_mod(a[0],mod-2);return 1;} 36 RG int r1=get_inv(a,ret,(ra+1)>>1),i,la=1; 37 while(la<(ra<<1))la<<=1; 38 memcpy(tmp3,a,ra<<2),memset(tmp3+ra,0,(la-ra)<<2); 39 memset(ret+r1,0,(la-r1)<<2);dft(tmp3,la,1),dft(ret,la,1); 40 for(i=0;i<la;++i) 41 ret[i]=((LL)ret[i]*(2ll+mod-((LL)ret[i]*tmp3[i]%mod)))%mod; 42 dft(ret,la,-1);return ra; 43 } 44 inline void rev_copy(int *st,int *to,int ra) 45 {for(RG int i=0;i<ra;++i)to[i]=st[ra-i-1];} 46 inline void reverse(int *st,int ra) 47 {for(RG int t,i=0,j=ra-1;i<j;++i,--j)t=st[i],st[i]=st[j],st[j]=t;} 48 inline int get_mod(int *a,int ra,int *p,int rp,int *ret) 49 { 50 while(ra&&!a[ra-1])--ra; 51 while(rp&&!p[rp-1])--rp; 52 if(ra<rp){memcpy(ret,a,ra<<2),memset(ret+ra,0,(rp-ra)<<2);return rp;} 53 static int tmp1[L],tmp2[L]; 54 RG int i,j,re=ra-rp+1,la=1; 55 while(la<(re<<1))la<<=1; 56 rev_copy(p,tmp1,rp);memset(tmp1+re,0,(la-re)<<2); 57 get_inv(tmp1,tmp2,re);memset(tmp2+re,0,(la-re)<<2); 58 rev_copy(a,tmp1,ra),memset(tmp1+re,0,(la-re)<<2); 59 dft(tmp1,la,1);dft(tmp2,la,1); 60 for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod; 61 dft(tmp1,la,-1); 62 la=1;while(la<ra)la<<=1; 63 reverse(tmp1,re);memset(tmp1+re,0,(la-re)<<2),dft(tmp1,la,1); 64 memcpy(tmp2,p,rp<<2),memset(tmp2+rp,0,(la-rp)<<2),dft(tmp2,la,1); 65 for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod; 66 dft(tmp1,la,-1); 67 for(i=0;i<rp;++i)ret[i]=(a[i]-tmp1[i]+mod)%mod; 68 memset(ret+rp,0,(la-rp)<<2); 69 while(rp&&!ret[rp-1])--rp;return rp; 70 } 71 inline void mul(int *a,int *b,int *p,int ra,int la) 72 { 73 RG int i,j,v,r1=ra+1; 74 memcpy(tmp,a,r1<<2);memset(tmp+r1,0,(la-r1)<<2); 75 memcpy(tmp2,b,r1<<2);memset(tmp2+r1,0,(la-r1)<<2); 76 dft(tmp,la,1),dft(tmp2,la,1); 77 for(i=0;i<la;++i)tmp[i]=(LL)tmp[i]*tmp2[i]%mod; 78 dft(tmp,la,-1); 79 get_mod(tmp,ra<<1+1,p,ra+1,a); 80 } 81 inline int calc(int k) 82 { 83 if(!k)return quick_mod(q,n); 84 RG int i,j,u,lim,la=1,ra,r1=k+2,ret=0; 85 memset(g,0,sizeof(g)),memset(h,0,sizeof(h)); 86 g[k][1]=(LL)pp[k]*q%mod; 87 g[k][0]=h[k][0]=1; 88 for(i=k-1;i>0;--i) 89 { 90 ra=k/i,g[i][0]=h[i][0]=1; 91 for(j=1;j<=ra;++j) 92 for(u=0;u<j;++u) 93 h[i][j]=(h[i][j]+(LL)h[i][u]*g[i+1][j-u-1]%mod*pp[i]%mod*q%mod)%mod; 94 for(j=1;j<=ra;++j) 95 for(u=0;u<=j;++u) 96 g[i][j]=(g[i][j]+(LL)h[i][u]*g[i+1][j-u]%mod)%mod; 97 } 98 memset(f,0,sizeof(f)),ra=k+1,f[0]=1; 99 for(i=1;i<=(ra<<1);++i) 100 for(j=1,lim=min(k+1,i);j<=lim;++j) 101 f[i]=(f[i]+(LL)f[i-j]*g[1][j-1]%mod*q%mod)%mod; 102 if(n<=(ra<<1)) 103 { 104 for(i=max(0,n-k);i<=n;++i) 105 ret=(ret+(LL)f[i]*g[1][n-i]%mod)%mod; 106 return ret; 107 } 108 memset(c,0,sizeof(c)),c[1]=1; 109 memset(d,0,sizeof(d)),d[ra]=1; 110 for(i=0;i<ra;++i)d[i]=(mod-(LL)q*g[1][k-i]%mod)%mod; 111 while(la<(r1<<1))la<<=1; 112 memset(e,0,sizeof(e)),e[0]=1; 113 for(lim=n-ra;lim;mul(c,c,d,ra,la),lim>>=1) 114 if(lim&1)mul(e,c,d,ra,la); 115 memset(ff,0,sizeof(ff)); 116 for(i=1;i<=ra;++i) 117 for(j=0;j<=k;++j) 118 ff[i]=(ff[i]+(LL)e[j]*f[i+j]%mod)%mod; 119 for(i=1;i<=ra;++i) 120 ret=(ret+(LL)ff[i]*g[1][ra-i])%mod; 121 return ret; 122 } 123 int main() 124 { 125 RG int i,x,y; 126 scanf("%d%d%d%d",&n,&k,&x,&y); 127 logi[1]=0; 128 while(len<=(k+1<<1))len<<=1,logi[len]=logi[len>>1]+1; 129 poww[0]=poww[len]=1,poww[1]=quick_mod(3,(mod-1)/len); 130 for(i=2;i<len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod; 131 for(i=0;i<len;++i) 132 if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1); 133 else rev[i]=(rev[i>>1]>>1); 134 135 p=(LL)x*quick_mod(y,mod-2)%mod, 136 q=(LL)(y-x)*quick_mod(y,mod-2)%mod; 137 for(pp[0]=1,pp[1]=p,i=2;i<=k;++i)pp[i]=(LL)pp[i-1]*p%mod; 138 printf("%d\n",(calc(k)-calc(k-1)+mod)%mod); 139 } UOJ#316

六，codeforces553E

?留坑待填

七，bzoj5093

斯特林數真的是太神了!

我們考慮對于每一個點計算它在每種情況下的貢獻

容易發現每個點都是一樣的.....那么我們乘個n就行了

然后考慮 每個點的貢獻應該可以寫成這樣

$$2^{?C(n-1,2)?? ?}* \sum_{i=1} ^{n-1} i^{k}C(n-1,i)???$$

所以現在我們就要求后面那個$ \sum_{i=1} ^{n-1} i^{k}C(n-1,i)?? $

這時候最神的東西就出現了

$$n^{k}= \sum _{j=0} ^{k}? S(k,j) * j! *C(n,j) $$

上午我和Cooook說要自己推,結果推不出來.....

不知道這個式子的實際意義是什么.....

但是似乎 斯特林數的展開之后也是一堆形如$i^{k}$的形式

也許和這有關?

但是不管了,我們繼續推式子,

$$?\sum_{i=1} ^{n-1} i^{k}?C(n-1,i)??? ? ? ?$$

$$=\sum_{i=1} ^{n-1}?C(n-1,i)?? ?\sum _{j=0} ^{k}? S(k,j) * j! *C(n,j)?$$

$$=\sum _{j=0} ^{n-1}? S(k,j) * j!? ?\sum_{i=j} ^{n-1} C(n-1,i)? *C(n,j)?$$

那么后面那個$C*C$的是什么玩意?

你想一想實際意義,相當于n個選i個,再i個選j個

相當于對于每j個元素,都會有$2^{n-j}$中不同的選取方案包含j

所以后面那東西就等于$C(n,j)*2^{n-j}$

所以我們的式子就變成了

$$=\sum _{j=0} ^{n-1}? S(k,j) * j!? C(n,j)*2^{n-j} $$
$$=\sum _{j=0} ^{n-1}? S(k,j) * A(n,j)*2^{n-j}?$$

這就很有意思了.....那個$A(n,j)$和$2^{n-j}$都可以$O(1)$維護

所以我們用fft優化求第k行斯特林數再按照式子打就行了

丟代碼:

1 #include <cstring> 2 #include <cstdio> 3 #include <algorithm> 4 using namespace std; 5 #define K 200010 6 #define LL long long 7 #define RG register 8 #define mod 998244353 9 #define L (1<<19)+10 10 int n,k,fac[K],invf[K],f[L],g[L]; 11 inline int min(int a,int b){return a<b?a:b;} 12 inline int quick_mod(int di,int mi) 13 { 14 int ret=1; 15 for(;mi;mi>>=1,di=(LL)di*di%mod) 16 if(mi&1)ret=(LL)ret*di%mod; 17 return ret; 18 } 19 int len=1,rev[L],poww[L]; 20 inline void dft(int *a,int opt) 21 { 22 RG int i,j,l,d,tmp,*w,wn,*x,*y; 23 for(i=0;i<len;++i)if(i<rev[i])swap(a[i],a[rev[i]]); 24 for(d=2;d<=len;d<<=1) 25 for(i=0,l=(d>>1),wn=((opt==1)?(len/d):(-len/d));i<len;i+=d) 26 for(w=(poww+((opt==1)?0:len)),j=0,x=a+i,y=x+l;j<l;++x,++y,++j,w+=wn) 27 tmp=(LL)(*w)*(*y)%mod,*y=(*x+mod-tmp)%mod,*x=(*x+tmp)%mod; 28 if(opt==-1) 29 for(tmp=quick_mod(len,mod-2),i=0;i<len;++i)a[i]=(LL)a[i]*tmp%mod; 30 } 31 int main() 32 { 33 scanf("%d%d",&n,&k); 34 RG int i,j,v,u,ans=(LL)n*quick_mod( 2,( ((LL)n-1)*(n-2)/2 )%(mod-1) )%mod; 35 while(len<(k<<1))len<<=1; 36 for(i=0;i<len;++i) 37 if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1); 38 else rev[i]=(rev[i>>1]>>1); 39 poww[0]=poww[len]=1,poww[1]=quick_mod(3,(mod-1)/len); 40 for(i=2;i<len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod; 41 for(fac[0]=fac[1]=1,i=2;i<=k;++i)fac[i]=(LL)fac[i-1]*i%mod; 42 for(invf[k]=quick_mod(fac[k],mod-2),i=k;i;--i)invf[i-1]=(LL)invf[i]*i%mod; 43 for(f[0]=0,i=1;i<=k;++i)f[i]=(LL)quick_mod(i,k)*invf[i]%mod; 44 for(i=0,v=1;i<k;++i,v=mod-v)g[i]=(LL)v*invf[i]%mod; 45 dft(f,1),dft(g,1); 46 for(i=0;i<len;++i)f[i]=(LL)f[i]*g[i]%mod; 47 dft(f,-1),--n; 48 RG int tn=n,lim=min(n,k),tmp=0,inv2=499122177; 49 for(v=1,u=quick_mod(2,n),i=0;i<=lim;++i) 50 tmp=(tmp+(LL)f[i]*v%mod*u)%mod,u=(LL)u*inv2%mod,v=(LL)v*tn%mod,--tn; 51 printf("%d\n",(LL)ans*tmp%mod ); 52 } bzoj5093

四，寫在最后

FFT的東西很多很雜也很難……有不少思想和題目都非常優秀的題目！

我現在學到的東西也只是滄海一粟而已，如果有不足歡迎指出！

轉載于:https://www.cnblogs.com/LadyLex/p/8335420.html

總結

以上是生活随笔為你收集整理的FFT小结的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。