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的题解。

给个代码：

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cstdlib>
using namespace std;
#define N 400010
#define L 1100000
#define mod 998244353
#define phi 998244352
#define g 3
#define LL long long
inline int max(int a,int b){return a>b?a:b;}
inline int min(int a,int b){return a<b?a:b;}
char s[N];
int f[N],pf[N],n,logi[N],litlen=1,len=1;
int bin[25],poww[L],rev[L],fac[L],inv[L],invf[L];
inline int quick_mod(int di,int mi)
{
    int ret=1;
    for(;mi;mi>>=1,di=(LL)di*di%mod)
        if(mi&1)ret=(LL)ret*di%mod;
    return ret;
}
inline void upd(int &a,int b){a+=b;if(a>=mod)a-=mod;}
inline int sub(int a,int b){a-=b;if(a<0)a+=mod;return a;}
inline int dft(int *a,int opt,int le)
{
    register int i,j,d,l,wn,w,tmp;
    for(i=0;i<le;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(d=2;d<=le;d<<=1)
        for(wn=(opt==1)?poww[len/d]:poww[len-len/d],l=d>>1,i=0;i<le;i+=d)
            for(w=1,j=0;j<l;++j,w=(LL)w*wn%mod)
                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;
                // tmp=(LL)a[i+j+l]*w%mod,a[i+j+l]=sub(a[i+j],tmp),upd(a[i+j],tmp);
    if(opt==-1)
        for(i=0;i<le;++i)a[i]=(LL)a[i]*inv[le]%mod;
}
int ret[L],C[L],A[L],B[L];
inline void CDQ(int l,int r)
{
    if(l==r)
        if(l<n){printf("%d
",(LL)f[l]*fac[l-1]%mod);return;}
        else {printf("%d
",(LL)f[l]*fac[l-1]%mod);exit (0);}
    register int le=r-l+1,i,mi=(l+r)>>1;
    CDQ(l,mi);
    for(i=0;i<le;++i)
        if(i&1)rev[i]=(rev[i>>1]>>1)|(le>>1);
        else rev[i]=(rev[i>>1])>>1;
    memset(A,0,le<<2),memset(B,0,le<<2);
    for(i=l;i<=mi;++i)A[i-l]=pf[i];
    for(i=0;i<r-l;++i)B[i]=C[i];
    dft(A,1,le),dft(B,1,le);
    for(i=0;i<le;++i)ret[i]=(LL)A[i]*B[i]%mod;
    dft(ret,-1,le);
    for(i=mi+1;i<=r;++i)upd(f[i],ret[i-l-1]);
    memset(A,0,le<<2),memset(B,0,le<<2);
    for(i=l;i<=mi;++i)A[i-l]=(LL)f[i]*inv[i]%mod;
    for(i=0;i<=r-l;++i)
        if(i<l)B[i]=(LL)f[i]*inv[i]*2%mod;
        else if(i<=mi)B[i]=(LL)f[i]*inv[i]%mod;
    dft(A,1,le),dft(B,1,le);
    for(i=0;i<le;++i)ret[i]=(LL)A[i]*B[i]%mod;
    dft(ret,-1,le);
    for(i=mi+1;i<=r;++i)upd(pf[i],ret[i-l]);
    CDQ(mi+1,r);
}
int main()
{
    register int i,j;
    for(bin[0]=i=1;i<=20;++i)bin[i]=bin[i-1]<<1;
    scanf("%d%s",&n,s);
    while(len<(n<<1))len<<=1;
    poww[0]=1,poww[1]=quick_mod(g,phi/len);
    for(i=2;i<=len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod;
    for(fac[0]=fac[1]=1,i=2;i<=len;++i)fac[i]=(LL)fac[i-1]*i%mod;
    for(inv[0]=inv[1]=1,i=2;i<=len;++i)inv[i]=(LL)inv[mod%i]*(mod-mod/i)%mod;
    for(invf[0]=invf[1]=1,i=2;i<=len;++i)invf[i]=(LL)invf[i-1]*inv[i]%mod;
    for(i=0;i<n;++i)
        if(s[i]=='1')C[i]=(LL)invf[i]*inv[2]%mod;
    f[1]=1,CDQ(1,len>>1);
}

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

这题我在课件里面讲了！

然后写的自己还算满意吧……再简单说一下

总之你可以发现每个数再若干次操作之后可以写成$frac {Ax+C}{Bx+D}$的形式，x代表每个数的初始值。然后……你把它写丑一点，可以变成形如$e+f*frac{1}{x+g}$的形式……但是具体内容请自己推导吧……

那么总和就是$n*e+f*sumfrac{1}{x+g}$，问题转化为维护$sumfrac{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

#include <cstdio>
#include <cstring>
#include <ctime>
#include <algorithm>
using namespace std;
#define mod 998244353
#define LL long long
#define RG register
#define N 100010
#define M 60010
char cB[1<<15],*S=cB,*T=cB;
#define getc (S==T&&(T=(S=cB)+fread(cB,1,1<<15,stdin),S==T)?0:*S++)
inline int read()
{
    int x=0;RG char c=getc;
    while(c<'0'|c>'9')c=getc;
    while(c>='0'&c<='9')x=10*x+(c^48),c=getc;
    return x;
}
inline int quick_mod(int di,int mi)
{
    int ret=1;
    for(;mi;mi>>=1,di=(LL)di*di%mod)
        if(mi&1)ret=(LL)ret*di%mod;
    return ret;
}
#define L (1<<18)+10
#define LM (1<<16)+10
int logi[L],inv[L],bin[25],rev[L],poww[L],len=1,ro=1;
inline void dft(int *a,int n,int opt)
{
    RG int i,j,d=logi[len]-logi[n],wn,l,tmp,*x,*y,*w;
    for(i=0;i<n;++i)if(i<(rev[i]>>d))swap(a[i],a[rev[i]>>d]);
    for(d=2;d<=n;d<<=1)
        for(i=0,l=d>>1,wn=(opt==1)?(len/d):(-len/d);i<n;i+=d)
            for(w=poww+((opt==1)?0:len),j=0,x=a+i+j,y=x+l;j<l;++j,w+=wn,++x,++y)
                tmp=(LL)(*w)*(*y)%mod,*y=(*x-tmp+mod)%mod,*x=(*x+tmp)%mod;
    if(opt==-1)for(i=0;i<n;++i)a[i]=(LL)a[i]*inv[n]%mod;
}
int n,m,cnt,sum,x[N];
int g[L<<1],h[L<<1];
int sta[LM],ans[M];
struct node{int e,f,g,id;}s[M];
inline void ins(int a,int b,int c,int d,int id)
{
    if(b==0)
        {ans[id]=(((LL)a*sum)+((LL)c*n))%mod*(LL)quick_mod(d,mod-2)%mod;return;}
    s[++cnt].f=((LL)b*c-(LL)a*d)%mod,
    b=quick_mod(b,mod-2);
    s[cnt].e=(LL)a*b%mod,s[cnt].g=(LL)d*b%mod,
    b=(LL)b*b%mod,s[cnt].f=(LL)s[cnt].f*b%mod;
    if(s[cnt].f<0)s[cnt].f+=mod;
    s[cnt].id=id;
}
inline int solve0(int *a,int l,int r)
{
    RG int ra=r-l+2;
    if(ra<=128)
    {
        memset(a,0,ra<<2),a[0]=1;
        RG int i,j,v;
        for(i=l;i<=r;++i)
            for(v=x[i],j=i-l;~j;--j)
                a[j+1]=(a[j+1]+a[j])%mod,
                a[j]=(LL)a[j]*v%mod;
        return ra;
    }
    RG int i,mi=l+r>>1,la=1;
    while(la<ra)la<<=1;
    RG int *f1=a,r1=solve0(f1,l,mi);
    RG int *f2=a+r1,r2=solve0(f2,mi+1,r);
    static int tmp1[L],tmp2[L];
    memcpy(tmp1,f1,r1<<2),memset(tmp1+r1,0,(la-r1)<<2),dft(tmp1,la,1);
    memcpy(tmp2,f2,r2<<2),memset(tmp2+r2,0,(la-r2)<<2),dft(tmp2,la,1);
    for(i=0;i<la;++i)a[i]=(LL)tmp1[i]*tmp2[i]%mod;
    dft(a,la,-1);
    return ra;
}
int *p[L];
inline int solve1(int id,int l,int r)
{
    static int mem[LM<<4],*head=mem;
    RG int ra=r-l+2;
    if(ra<=128)
    {
        RG int i,j,v,*f=p[id]=head;head+=ra;
        memset(f,0,ra<<2),f[0]=1;
        for(i=l;i<=r;++i)
            for(v=(mod-sta[i])%mod,j=i-l;~j;--j)
                f[j+1]=(f[j+1]+f[j])%mod,f[j]=(LL)f[j]*v%mod;
        return ra;
    }
    RG int mi=l+r>>1,la=1;
    while(la<ra)la<<=1;
    RG int r1=solve1(id<<1,l,mi),*f1=p[id<<1];
    RG int r2=solve1(id<<1|1,mi+1,r),*f2=p[id<<1|1];
    static int tmp1[LM],tmp2[LM];
    memcpy(tmp1,f1,r1<<2),memset(tmp1+r1,0,(la-r1)<<2),dft(tmp1,la,1);
    memcpy(tmp2,f2,r2<<2),memset(tmp2+r2,0,(la-r2)<<2),dft(tmp2,la,1);
    RG int i,*f=p[id]=head;head+=ra;
    for(i=0;i<la;++i)f[i]=(LL)tmp1[i]*tmp2[i]%mod;
    dft(f,la,-1);return ra;
}
inline int get_inv(int *a,int *ret,int ra)
{
    if(ra==1){ret[0]=quick_mod(a[0],mod-2);return 1;}
    static int tmp[L];
    RG int r1=get_inv(a,ret,(ra+1)>>1),i,la=1;
    while(la<(ra<<1))la<<=1;
    memcpy(tmp,a,ra<<2),memset(tmp+ra,0,(la-ra)<<2);
    memset(ret+r1,0,(la-r1)<<2);
    dft(tmp,la,1),dft(ret,la,1);
    for(i=0;i<la;++i)
    {
        ret[i]=((LL)ret[i]*(2ll-((LL)ret[i]*tmp[i]%mod)))%mod;
        if(ret[i]<0)ret[i]+=mod;
    }
    dft(ret,la,-1);return ra;
}
inline void rev_copy(int *st,int *to,int ra)
    {for(RG int i=0;i<ra;++i)to[i]=st[ra-i-1];}
inline void reverse(int *st,int ra)
    {for(RG int t,i=0,j=ra-1;i<j;++i,--j)t=st[i],st[i]=st[j],st[j]=t;}
inline int get_mod(int *a,int ra,int *p,int rp,int *ret)
{
    while(ra&&!a[ra-1])--ra;
    while(rp&&!p[rp-1])--rp;
    if(ra<rp){memcpy(ret,a,ra<<2),memset(ret+ra,0,(rp-ra)<<2);return rp;}
    static int tmp1[L],tmp2[L];
    RG int i,j,re=ra-rp+1,la=1;
    while(la<(re<<1))la<<=1;
    memset(tmp1,0,la<<2),rev_copy(p,tmp1,rp),get_inv(tmp1,tmp2,re);
    memset(tmp2+re,0,(la-re)<<2),dft(tmp2,la,1);
    rev_copy(a,tmp1,ra),memset(tmp1+re,0,(la-re)<<2),dft(tmp1,la,1);
    for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod;
    dft(tmp1,la,-1);
    la=1;while(la<ra)la<<=1;
    reverse(tmp1,re),memset(tmp1+re,0,(la-re)<<2),dft(tmp1,la,1);
    memcpy(tmp2,p,rp<<2),memset(tmp2+rp,0,(la-rp)<<2),dft(tmp2,la,1);
    for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod;
    dft(tmp1,la,-1);
    for(i=0;i<rp;++i)ret[i]=(a[i]-tmp1[i]+mod)%mod;
    memset(ret+rp,0,(la-rp)<<2);
    while(rp&&!ret[rp-1])--rp;
    return rp;
}
int val0[LM],val[LM];
inline void solve2(int id,int *a,int *b,int l,int r,int deg)
{
    RG int ra=r-l+2;
    if(deg<=128)
    {
        RG int i,j,v,u,F,G;
        for(i=l;i<=r;val0[i]=F,val[i]=G,++i)
            for(u=sta[i],v=1,F=G=0,j=0;j<=deg;++j,v=(LL)v*u%mod)
                F=(F+(LL)v*a[j])%mod,G=(G+(LL)v*b[j])%mod;
        return;
    }
    RG int i,mi=l+r>>1;
    int r1=get_mod(a,deg,p[id],ra,a+deg);a+=deg;
    int r2=get_mod(b,deg,p[id],ra,b+deg);b+=deg;
    ra=min(r1,r2);
    solve2(id<<1,a,b,l,mi,ra),solve2(id<<1|1,a,b,mi+1,r,ra);
}
int main()
{
    RG int i,j,opt,v,A=1,B=0,C=0,D=1;
    n=read(),m=read();
    while(len<(n<<1))len<<=1;
    for(bin[0]=i=1;i<=20;++i)bin[i]=bin[i-1]<<1;
    for(i=0;i<=20;++i)
        inv[bin[i]]=quick_mod(bin[i],mod-2),logi[bin[i]]=i;
    for(int i=0;i<len;++i)
        if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1);
        else rev[i]=(rev[i>>1]>>1);
    poww[0]=1,poww[1]=quick_mod(3,(mod-1)/len);
    for(i=2;i<=len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod;
    for(sum=0,i=1;i<=n;++i)x[i]=read(),sum=(sum+x[i])%mod;
    for(i=1;i<=m;++i)
        if(read()==1)v=read(),A=(A+(LL)v*B)%mod,C=(C+(LL)v*D)%mod,ins(A,B,C,D,i);
        else swap(A,B),swap(C,D),ins(A,B,C,D,i);
    if(cnt)
    {
        for(i=1;i<=cnt;++i)sta[i]=s[i].g;
        sort(sta+1,sta+cnt+1),opt=unique(sta+1,sta+cnt+1)-sta-1;
        for(i=1;i<=cnt;++i)
            s[i].g=lower_bound(sta+1,sta+opt+1,s[i].g)-sta;
        solve0(g,1,n);
        for(h[n]=0,i=1;i<=n;++i)h[i-1]=(LL)g[i]*i%mod;
        solve1(1,1,opt),solve2(1,g,h,1,opt,n+1);
        for(i=1;i<=cnt;++i)
            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;
    }
    for(i=1;i<=m;++i)printf("%d
",ans[i]);
}

三，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后面就没那么难了……

丢代码：

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define LLL __int128
#define LL long long
#define RG register
#define PMAX 30010
#define mod 998244353
inline int quick_mod(int di,int mi)
{
    RG int ret=1;
    for(;mi;mi>>=1,di=(LL)di*di%mod)
        if(mi&1)ret=(LL)ret*di%mod;
    return ret;
}
#define L (1<<16)+10
int len=1,rev[L],poww[L];
inline void dft(int *a,int opt)
{
    RG int i,j,d,l,*w,wn,tmp,*x,*y;
    for(i=0;i<len;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(d=2;d<=len;d<<=1)
        for(i=0,l=(d>>1),wn=(opt==1?(len/d):(-len/d));i<len;i+=d)
            for(j=0,w=(poww+(opt==1?0:len)),x=a+i,y=x+l;j<l;++j,++x,++y,w+=wn)
                tmp=(LL)(*w)*(*y)%mod,*y=(*x+mod-tmp)%mod,*x=(*x+tmp)%mod;
    if(opt==-1)
        for(tmp=quick_mod(len,mod-2),i=0;i<len;++i)a[i]=(LL)a[i]*tmp%mod;
}
int p,phi,root,logv[PMAX],expg[PMAX];
int fac[PMAX],invf[PMAX];
inline int C(int a,int b)
    {return (a<b)?0:(fac[a]*invf[b]%p*invf[a-b]%p);}
inline int quick_p(int di,int mi)
{
    RG int ret=1;
    for(;mi;mi>>=1,di=di*di%p)
        if(mi&1)ret=ret*di%p;
    return ret;
}
inline void getroot()
{
    RG int i,j,v;
    for(i=1;i<p;++i)
    {
        memset(logv,-1,sizeof(logv));
        logv[1]=0,expg[0]=1;
        for(v=i,j=1;j<phi;v=v*i%p,++j)
            {if(logv[v]>=0)break;expg[j]=v,logv[v]=j;}
        if(j==phi){root=i;return;}
    }
}
char B[1<<15],*S=B,*T=B;
#define getc (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++)
template <typename _T>
inline _T read()
{
    RG _T x=0;RG char c=getc;
    while(c<'0'|c>'9')c=getc;
    while(c>='0'&&c<='9')x=10*x+(c^48),c=getc;
    return x;
}
int lenn,bn[300];
int f[300][L],g[300][L];
int ans[PMAX],tmp[L],tmp1[L],tmp2[L];
inline void calc(LLL n,int opt)
{
    if(n<=0)return;
    RG int i,j,v,k,bin[300],to,bit=0;
    memset(bin,0,sizeof(bin));
    while(n)bin[++bit]=n%p,n/=p;
    RG int lim=max(bit,lenn);
    memset(f,0,sizeof(f)),f[0][0]=1;
    if(opt==1)
    {
        memset(g,0,sizeof(g)),g[0][0]=1;
        for(i=1;i<=lim;++i)
        {
            memset(tmp1,0,sizeof(tmp1));
            for(j=bn[i];j<p;++j)
                {v=logv[C(j,bn[i])];if(v!=-1)++tmp1[v];}
            memcpy(tmp2,g[i-1],phi<<2);
            memset(tmp2+phi,0,(len-phi)<<2);
            dft(tmp1,1),dft(tmp2,1);
            for(j=0;j<len;++j)tmp[j]=(LL)tmp1[j]*tmp2[j]%mod;
            dft(tmp,-1);
            for(j=0;j<phi;++j)
                g[i][j]=(tmp[j]+tmp[j+phi])%mod;
        }
    }
    for(i=1;i<=lim;++i)
    {
        memset(tmp1,0,sizeof(tmp1));
        for(j=bn[i];j<bin[i];++j)
        {
            v=logv[C(j,bn[i])];
            if(v!=-1)++tmp1[v];
        }
        memcpy(tmp2,g[i-1],phi<<2);
        memset(tmp2+phi,0,(len-phi)<<2);
        dft(tmp1,1);dft(tmp2,1);
        for(j=0;j<len;++j)tmp[j]=(LL)tmp1[j]*tmp2[j]%mod;
        dft(tmp,-1);    
        for(j=0;j<phi;++j)f[i][j]=(tmp[j]+tmp[j+phi])%mod;
        v=logv[C(bin[i],bn[i])];
        if(v!=-1)for(k=0;k<phi;++k)
            to=(v+k)%phi,f[i][to]=(f[i][to]+f[i-1][k])%mod;
    }
    if(opt==1)for(i=0;i<phi;++i)ans[expg[i]]=(ans[expg[i]]+f[lim][i])%mod;
    else for(i=0;i<phi;++i)ans[expg[i]]=(ans[expg[i]]+mod-f[lim][i])%mod;
}
int main()
{
    RG int i,j,sum;
    p=read<int>(),phi=p-1,getroot();
    while(len<=(p<<1))len<<=1;
    for(i=0;i<len;++i)
        if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1);
        else rev[i]=rev[i>>1]>>1;
    poww[0]=poww[len]=1,poww[1]=quick_mod(3,(mod-1)/len);
    for(i=1;i<len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod;
    for(fac[0]=fac[1]=1,i=2;i<p;++i)fac[i]=fac[i-1]*i%p;
    invf[p-1]=quick_p(fac[p-1],p-2);
    for(i=p-1;i;--i)invf[i-1]=invf[i]*i%p;
    LLL n,l,r;
    n=read<LLL>(),l=read<LLL>(),r=read<LLL>();
    while(n)bn[++lenn]=n%p,n/=p;
    calc(r,1),calc(l-1,-1);
    for(sum=(r-l+1)%mod,i=1;i<p;++i)sum=(sum-ans[i]+mod)%mod;
    for(printf("%d
",sum),i=1;i<p;++i)printf("%d
",ans[i]);
}

四，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

那么我们扔代码

#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;
#define mod 998244353
#define K 1010
#define L 4096
#define RG register
#define LL long long
inline int quick_mod(int di,int mi)
{
    int ret=1;
    for(;mi;mi>>=1,di=(LL)di*di%mod)
        if(mi&1)ret=(LL)ret*di%mod;
    return ret;
}
int n,k,p,q;
int g[K][K],h[K][K],pp[K],f[K<<1],ff[K<<1];
int tmp[L],tmp1[L],tmp2[L],tmp3[L],c[L],d[L],e[L];
inline int min(int a,int b){return a<b?a:b;}
int poww[L],logi[L],rev[L],len=1;
inline void dft(int *a,int ra,int opt)
{
    register int i,j,l,d=logi[len]-logi[ra],wn,tmp,*w,*x,*y;
    for(i=0;i<ra;++i)if(i<(rev[i]>>d))swap(a[i],a[rev[i]>>d]);
    for(d=2;d<=ra;d<<=1)
        for(wn=((opt==1)?(len/d):(-len/d)),i=0,l=(d>>1);i<ra;i+=d)
            for(w=poww+(opt==1?0:len),j=0,x=a+i,y=x+l;j<l;++j,++x,++y,w+=wn)
                tmp=(LL)(*w)*(*y)%mod,*y=(*x-tmp+mod)%mod,*x=(*x+tmp)%mod;
    if(opt==-1)
        for(tmp=quick_mod(ra,mod-2),i=0;i<ra;++i)a[i]=(LL)a[i]*tmp%mod;
}
inline int get_inv(int *a,int *ret,int ra)
{
    if(ra==1){ret[0]=quick_mod(a[0],mod-2);return 1;}
    RG int r1=get_inv(a,ret,(ra+1)>>1),i,la=1;
    while(la<(ra<<1))la<<=1;
    memcpy(tmp3,a,ra<<2),memset(tmp3+ra,0,(la-ra)<<2);
    memset(ret+r1,0,(la-r1)<<2);dft(tmp3,la,1),dft(ret,la,1);
    for(i=0;i<la;++i)
        ret[i]=((LL)ret[i]*(2ll+mod-((LL)ret[i]*tmp3[i]%mod)))%mod;
    dft(ret,la,-1);return ra;
}
inline void rev_copy(int *st,int *to,int ra)
    {for(RG int i=0;i<ra;++i)to[i]=st[ra-i-1];}
inline void reverse(int *st,int ra)
    {for(RG int t,i=0,j=ra-1;i<j;++i,--j)t=st[i],st[i]=st[j],st[j]=t;}
inline int get_mod(int *a,int ra,int *p,int rp,int *ret)
{
    while(ra&&!a[ra-1])--ra;
    while(rp&&!p[rp-1])--rp;
    if(ra<rp){memcpy(ret,a,ra<<2),memset(ret+ra,0,(rp-ra)<<2);return rp;}
    static int tmp1[L],tmp2[L];
    RG int i,j,re=ra-rp+1,la=1;
    while(la<(re<<1))la<<=1;
    rev_copy(p,tmp1,rp);memset(tmp1+re,0,(la-re)<<2);
    get_inv(tmp1,tmp2,re);memset(tmp2+re,0,(la-re)<<2);
    rev_copy(a,tmp1,ra),memset(tmp1+re,0,(la-re)<<2);
    dft(tmp1,la,1);dft(tmp2,la,1);
    for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod;
    dft(tmp1,la,-1);
    la=1;while(la<ra)la<<=1;
    reverse(tmp1,re);memset(tmp1+re,0,(la-re)<<2),dft(tmp1,la,1);
    memcpy(tmp2,p,rp<<2),memset(tmp2+rp,0,(la-rp)<<2),dft(tmp2,la,1);
    for(i=0;i<la;++i)tmp1[i]=(LL)tmp1[i]*tmp2[i]%mod;
    dft(tmp1,la,-1);
    for(i=0;i<rp;++i)ret[i]=(a[i]-tmp1[i]+mod)%mod;
    memset(ret+rp,0,(la-rp)<<2);
    while(rp&&!ret[rp-1])--rp;return rp;
}
inline void mul(int *a,int *b,int *p,int ra,int la)
{
    RG int i,j,v,r1=ra+1;
    memcpy(tmp,a,r1<<2);memset(tmp+r1,0,(la-r1)<<2);
    memcpy(tmp2,b,r1<<2);memset(tmp2+r1,0,(la-r1)<<2);
    dft(tmp,la,1),dft(tmp2,la,1);
    for(i=0;i<la;++i)tmp[i]=(LL)tmp[i]*tmp2[i]%mod;
    dft(tmp,la,-1);
    get_mod(tmp,ra<<1+1,p,ra+1,a);
}
inline int calc(int k)
{
    if(!k)return quick_mod(q,n);
    RG int i,j,u,lim,la=1,ra,r1=k+2,ret=0;
    memset(g,0,sizeof(g)),memset(h,0,sizeof(h));
    g[k][1]=(LL)pp[k]*q%mod;
    g[k][0]=h[k][0]=1;
    for(i=k-1;i>0;--i)
    {
        ra=k/i,g[i][0]=h[i][0]=1;
        for(j=1;j<=ra;++j)
            for(u=0;u<j;++u)
                h[i][j]=(h[i][j]+(LL)h[i][u]*g[i+1][j-u-1]%mod*pp[i]%mod*q%mod)%mod;
        for(j=1;j<=ra;++j)
            for(u=0;u<=j;++u)
                g[i][j]=(g[i][j]+(LL)h[i][u]*g[i+1][j-u]%mod)%mod;
    }
    memset(f,0,sizeof(f)),ra=k+1,f[0]=1;
    for(i=1;i<=(ra<<1);++i)
        for(j=1,lim=min(k+1,i);j<=lim;++j)
            f[i]=(f[i]+(LL)f[i-j]*g[1][j-1]%mod*q%mod)%mod;
    if(n<=(ra<<1))
    {
        for(i=max(0,n-k);i<=n;++i)
            ret=(ret+(LL)f[i]*g[1][n-i]%mod)%mod;
        return ret;
    }
    memset(c,0,sizeof(c)),c[1]=1;
    memset(d,0,sizeof(d)),d[ra]=1;
    for(i=0;i<ra;++i)d[i]=(mod-(LL)q*g[1][k-i]%mod)%mod;
    while(la<(r1<<1))la<<=1;
    memset(e,0,sizeof(e)),e[0]=1;
    for(lim=n-ra;lim;mul(c,c,d,ra,la),lim>>=1)
        if(lim&1)mul(e,c,d,ra,la);
    memset(ff,0,sizeof(ff));
    for(i=1;i<=ra;++i)
        for(j=0;j<=k;++j)
            ff[i]=(ff[i]+(LL)e[j]*f[i+j]%mod)%mod;
    for(i=1;i<=ra;++i)
        ret=(ret+(LL)ff[i]*g[1][ra-i])%mod;
    return ret;
}
int main()
{
    RG int i,x,y;
    scanf("%d%d%d%d",&n,&k,&x,&y);
    logi[1]=0;
    while(len<=(k+1<<1))len<<=1,logi[len]=logi[len>>1]+1;
    poww[0]=poww[len]=1,poww[1]=quick_mod(3,(mod-1)/len);
    for(i=2;i<len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod;
    for(i=0;i<len;++i)
        if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1);
        else rev[i]=(rev[i>>1]>>1);

    p=(LL)x*quick_mod(y,mod-2)%mod,
    q=(LL)(y-x)*quick_mod(y,mod-2)%mod;
    for(pp[0]=1,pp[1]=p,i=2;i<=k;++i)pp[i]=(LL)pp[i-1]*p%mod;
    printf("%d
",(calc(k)-calc(k-1)+mod)%mod);
}

六，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行斯特林数再按照式子打就行了

丢代码:

#include <cstring>
#include <cstdio>
#include <algorithm>
using namespace std;
#define K 200010
#define LL long long
#define RG register
#define mod 998244353
#define L (1<<19)+10
int n,k,fac[K],invf[K],f[L],g[L];
inline int min(int a,int b){return a<b?a:b;}
inline int quick_mod(int di,int mi)
{
    int ret=1;
    for(;mi;mi>>=1,di=(LL)di*di%mod)
        if(mi&1)ret=(LL)ret*di%mod;
    return ret;
}
int len=1,rev[L],poww[L];
inline void dft(int *a,int opt)
{
    RG int i,j,l,d,tmp,*w,wn,*x,*y;
    for(i=0;i<len;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(d=2;d<=len;d<<=1)
        for(i=0,l=(d>>1),wn=((opt==1)?(len/d):(-len/d));i<len;i+=d)
            for(w=(poww+((opt==1)?0:len)),j=0,x=a+i,y=x+l;j<l;++x,++y,++j,w+=wn)
                tmp=(LL)(*w)*(*y)%mod,*y=(*x+mod-tmp)%mod,*x=(*x+tmp)%mod;
    if(opt==-1)
        for(tmp=quick_mod(len,mod-2),i=0;i<len;++i)a[i]=(LL)a[i]*tmp%mod;
}
int main()
{
    scanf("%d%d",&n,&k);
    RG int i,j,v,u,ans=(LL)n*quick_mod(  2,(  ((LL)n-1)*(n-2)/2  )%(mod-1)  )%mod;
    while(len<(k<<1))len<<=1;
    for(i=0;i<len;++i)
        if(i&1)rev[i]=(rev[i>>1]>>1)|(len>>1);
        else rev[i]=(rev[i>>1]>>1);
    poww[0]=poww[len]=1,poww[1]=quick_mod(3,(mod-1)/len);
    for(i=2;i<len;++i)poww[i]=(LL)poww[i-1]*poww[1]%mod;
    for(fac[0]=fac[1]=1,i=2;i<=k;++i)fac[i]=(LL)fac[i-1]*i%mod;
    for(invf[k]=quick_mod(fac[k],mod-2),i=k;i;--i)invf[i-1]=(LL)invf[i]*i%mod;
    for(f[0]=0,i=1;i<=k;++i)f[i]=(LL)quick_mod(i,k)*invf[i]%mod;
    for(i=0,v=1;i<k;++i,v=mod-v)g[i]=(LL)v*invf[i]%mod;
    dft(f,1),dft(g,1);
    for(i=0;i<len;++i)f[i]=(LL)f[i]*g[i]%mod;
    dft(f,-1),--n;
    RG int tn=n,lim=min(n,k),tmp=0,inv2=499122177;
    for(v=1,u=quick_mod(2,n),i=0;i<=lim;++i)
        tmp=(tmp+(LL)f[i]*v%mod*u)%mod,u=(LL)u*inv2%mod,v=(LL)v*tn%mod,--tn;
    printf("%d
",(LL)ans*tmp%mod );
}

四，写在最后

FFT的东西很多很杂也很难……有不少思想和题目都非常优秀的题目！

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