FFT小結

時間 2019-12-10

標籤 fft 小結简体版

原文原文鏈接

FFT小結！

零，說在前面

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

……使人毒瘤的FFThtml

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

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

那麼話很少說咱們開始FFT的總結.數組

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

你們請善用搜索引擎啦......ide

一，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 }