FFT小结

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的东西很多很杂也很难……有不少思想和题目都非常优秀的题目!

我现在学到的东西也只是沧海一粟而已,如果有不足欢迎指出!

posted @ 2018-02-27 19:18 LadyLex 阅读(...) 评论(...) 编辑 收藏