多项式重工业修炼日志

没错,就是这样,用来记录修炼进度。
主要是用来存板子,有些理解了比较长时间的地方会大概讲讲(不过越写越成过程推导了)。
在形式幂级数上定义各种运算,有助于我们更快地解决一些问题,也能更好地为生成函数运算提供快速实现方式。


FFT (Fast Fourier Transform)

一切的开始。
DFT 只需要记住这两个式子就行了:

[F(omega_n^k)=F_l(omega_{n/2}^k)+omega_n^kF_r(omega_{n/2}^k) ]

[F(omega_n^{k+n/2})=F_l(omega_{n/2}^k)-omega_n^kF_r(omega_{n/2}^k) ]

IDFT 的思想即是单位根反演。记 (G={ m DFT}(F)),即

[G[k]=sum_{i=0}^{n-1}F[i](omega_n^k)^iLongleftrightarrow nF[k]=sum_{i=0}^{n-1}G[i](omega_n^{-k})^i ]

//luogu3803
#include <bits/stdc++.h>
using namespace std;

const int N=2e6+3e5;
const double pi=acos(-1);
struct cplx
{
    double x,y;
    cplx(double _x=0,double _y=0) {x=_x,y=_y;}
    cplx operator+(const cplx &C) {return cplx(x+C.x,y+C.y);}
    cplx operator-(const cplx &C) {return cplx(x-C.x,y-C.y);}
    cplx operator*(const cplx &C) {return cplx(x*C.x-y*C.y,x*C.y+y*C.x);}
}f[N],g[N];
int n,m,rv[N];

void FFT(cplx* f,int type)
{
    for(int i=0;i<n;++i) if(i<rv[i]) swap(f[i],f[rv[i]]);
    for(int p=2;p<=n;p<<=1)
    {
        int len=p>>1;
        cplx w1=cplx(cos(pi/len),type*sin(pi/len));
        for(int i=0;i<n;i+=p)
        {
            cplx wk=cplx(1,0);
            for(int k=i;k<i+len;++k)
            {
                cplx tmp=wk*f[k+len];
                f[k+len]=f[k]-tmp,f[k]=f[k]+tmp;
                wk=wk*w1;
            }
        }
    }
    if(type==-1) for(int i=0;i<n;++i) f[i].x/=n;
}

int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;++i) scanf("%lf",&f[i].x);
    for(int i=0;i<=m;++i) scanf("%lf",&g[i].x);
    for(m+=n,n=1;n<=m;n<<=1);
    for(int i=1;i<n;++i)
        rv[i]=(rv[i>>1]>>1)|((i&1)?n>>1:0);
    FFT(f,1); FFT(g,1);
    for(int i=0;i<n;++i) f[i]=f[i]*g[i];
    FFT(f,-1);
    for(int i=0;i<=m;++i) printf("%d ",(int)(f[i].x+0.49));
    return 0;
}

拉格朗日插值(Lagrange Interpolation)

众所周知,(n+1) 个不同点可以唯一确定一个 (n) 次多项式 (F(x))
解多项式有一个比较 naive 的做法是待定系数后高斯消元,但这么做复杂度高达 (O(n^3)) 且还有精度方面的问题。
而拉格朗日插值是一种比较优秀的方法,复杂度只有 (O(n^2))。特别地,如果给出的点横坐标连续,复杂度可以进一步降为 (O(n))
公式很简单:

[F(x)=sum_{i=0}^ny_iprod_{j e i}frac{x-x_j}{x_i-x_j} ]

可以发现对于每一个 ((x_i,y_i)),和式中仅有 (y_i) 那一项后面的乘积为 1,其他项均为零,满足这个多项式过这 (n+1) 个点的条件。于是如果我们要算 (F(k)) 的话直接带入即可。
有时候可能会有加点操作,如果仍然 (O(n^2)) 重构不太明智,可以考虑变个形

[F(x)=prod_{i=1}^n(x-x_i)sum_{i=1}^nfrac{omega_i}{x-x_i} ]

其中

[omega_i=frac{y_i}{prodlimits_{j e i}(x_i-x_j)} ]

这样每次加点时只需维护 (omega_i) 即可。

#include <bits/stdc++.h>
using namespace std;
const int N=2005,P=998244353;
int n,k,x[N],y[N],ans;

int qpow(int bas,int p)
{
    int res=1;
    for(;p;p>>=1,bas=1LL*bas*bas%P)
        if(p&1) res=1LL*res*bas%P;
    return res;
}

int main()
{
    n=IO::read(),k=IO::read();
    for(int i=0;i<n;++i) x[i]=IO::read(),y[i]=IO::read();
    for(int i=0;i<n;++i)
    {
        int up=1,down=1;
        for(int j=0;j<n;++j)
        if(i!=j) up=1LL*up*(k-x[j]+P)%P,down=1LL*down*(x[i]-x[j]+P)%P;
        ans=(ans+1LL*up*qpow(down,P-2)%P*y[i]%P)%P;
    }
    printf("%d",ans);
    return 0;
}

NTT (Number Theory Transform)

FFT 的缺点很明显:精度问题,三角函数太慢。
我们需要找到一种代替单位根的东西,它必须在形式上符合单位根的运算。
考虑原根这个东西,具体的理论不想展开说了(反正也用不上),只需记住对于形如 (p=k2^t+1) 的素数,设其模 (p) 意义下的原根为 (g),则 (g^0,g^1,dots,g^{p-1}) 互不相同。
容易验证,只需考虑 (omega^1_n=g^{frac{p-1}{n}}),则所有对于单位根的运算对原根也成立。
常见的 NTT 模数:(998244353=119 imes 2^{23}+1,2281701377=17 imes 2^{27}+1,1004535809=479 imes 2^{21}+1),它们的原根都是 (3)
代码与 FFT 对比基本没变,只是三角函数换成了原根。
所以我为啥比 FFT 慢啊/kk

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int G=3,P=998244353,N=1350000;

ll qpow(ll b,ll p=P-2)
{
    ll res=1;
    for(;p;p>>=1,b=b*b%P)
        if(p&1) res=res*b%P;
    return res;
}

const ll invG=qpow(G);
int n,m,tr[N*2];
ll f[N*2],g[N*2],invn;

void NTT(ll f[],bool op)
{
    for(int i=0;i<n;++i)
        if(i<tr[i]) swap(f[i],f[tr[i]]);
    for(int p=2;p<=n;p<<=1)
    {
        int len=p>>1,tG=qpow(op?G:invG,(P-1)/p);
        for(int k=0;k<n;k+=p)
        {
            ll bfy=1;
            for(int l=k;l<k+len;++l)
            {
                int tt=bfy*f[len+l]%P;
                f[len+l]=f[l]-tt,f[l]=f[l]+tt;
                if(f[len+l]<0) f[len+l]+=P;
                if(f[l]>P) f[l]-=P;
                bfy=bfy*tG%P;
            }
        }
    }
}

int main()
{
    scanf("%d%d",&n,&m); ++n,++m;
    for(int i=0;i<n;++i) scanf("%d",f+i);
    for(int i=0;i<m;++i) scanf("%d",g+i);
    for(m+=n,n=1;n<m;n<<=1);
    for(int i=0;i<n;++i)
        tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
    NTT(f,1); NTT(g,1);
    for(int i=0;i<n;++i) f[i]=f[i]*g[i]%P;
    NTT(f,0); invn=qpow(n);
    for(int i=0;i<m-1;++i) printf("%lld ",f[i]*invn%P);
    return 0;
}

多项式求逆(乘法逆)

(F(x)equiv G^{-1}(x)pmod {x^n})
我们考虑倍增求逆,记 (F=G^{-1}pmod {x^n},F'=G^{-1}pmod {x^{n/2}}),现在我们已经求出了 (F'),考虑推出 (F)

[egin{aligned}&F-F'equiv0pmod {x^{n/2}} \ iff& (F-F')^2equiv0pmod {x^n} \ iff& F^2-2FF'+F'^2equiv0pmod {x^n} \ iff& F-2F'+GF'^2equiv0pmod {x^n} \ iff& Fequiv 2F'-GF'^2pmod {x^n}end{aligned} ]

由主定理,复杂度为 (T(n)=T(n/2)+O(nlog n)=O(nlog n))
由于我们需要频繁地进行一些操作,适度地封装显得非常有必要。我的多项式全是从 cmd_block 大佬那学来的,然而他的码风实在是……一言难尽。所以以下是根据我的理解封装的,大家可以自行探索合适的封装方式。我抱跳瓜的大腿了,届时会白嫖他的板子,所以代码先咕着。

多项式开方

(F(x)^2equiv G(x)pmod {x^n})
继续考虑倍增,假设已求得 (F'^2equiv Gpmod {x^{n/2}}),于是有

[egin{aligned}&(F-F')(F+F')equiv 0 pmod{x^{n/2}} \ iff& F-F'equiv 0 pmod {x^{n/2}} \ iff& (F-F')^2equiv 0pmod{x^n} \ iff& F^2-2FF'+F'^2equiv 0 pmod{x^n} \ iff& G-2FF'+F'^2equiv 0pmod{x^n} \ iff& Fequiv frac{G+F'^2}{2F'}pmod{x^n}end{aligned} ]

可以看到推法和求逆很像。对于常数项,答案是模意义下的二次剩余,一般来讲做题时常数项都是 (1),如果不是可以骂死出题人

多项式求导(求积分)

这不用说了吧,这都不会做什么多项式啊。

多项式牛顿迭代

给多项式 (G),求 (F) 使得 (G(F(x))equiv 0pmod {x^n})
还是考虑倍增,若求得 (G(F'(x))equiv 0pmod {x^{n/2}}),直接泰勒展开:

[G(F(x))equivsum_{ige 0} frac{(F(x)-F'(x))^iG^{(n)}(F'(x))}{i!}pmod {x^n} ]

由于 (F') 最高次项是 (x^{n/2}) 的,所以对于 (ige 2) 的情况模意义为 (0)
于是有:

[egin{aligned}&G(F(x))equiv G(F'(x))+(F(x)-F'(x))G^{(1)}(F'(x))pmod {x^n} \ iff& F(x)equiv F'(x)-frac{G(F'(x))}{G^{(1)}(F'(x))}pmod{x^n}end{aligned} ]

(然后你会发现这就是普通的牛顿迭代式子)
有了多项式牛迭,我们可以无脑求解很多式子。

多项式 ln

即求 (F(x)=ln G(x))
考虑两边同时求导。但要注意两边的主元是不一样的,右边实际上是一个复合函数求导,应用链式法则,有(注意这里多项式带撇表示求导)

[egin{aligned}& F'(x)= frac{{ m d}ln G(x)}{{ m d}G(x)}G'(x) \ iff& F'(x)= G^{-1}(x)G'(x) \ iff& F(x)= int G^{-1}(x)G'(x)end{aligned} ]

四部曲,求逆,求导,乘起来,求积分。

多项式 exp

即求 (F(x)=exp G(x))
两边先同时取对数,变成 (ln F(x)equiv G(x)pmod{x^n}),移项得 (ln F(x)-G(x)equiv 0pmod {x^n})。如果设 (H(F(x))=ln F(x)-G(x)),我们发现实际要求的其实是这个函数的零点。
所以我们考虑用牛顿迭代来倍增求解。设已求得 (F'(x)equivexp G(x)pmod{x^{n/2}}),则有

[egin{aligned}&F(x)equiv F'(x)-frac{H(F'(x))}{H^{(1)}(F'(x))}pmod {x^n} \ iff& F(x)equiv F'(x)(1-ln F'(x)+G(x))pmod {x^n}end{aligned} ]

注意初值,当 (F) 为常数项时值为 (1)

多项式快速幂

即求 (F^k(x))
考虑到 (F^k(x)=exp( kln F(x))),先 (ln) 回来,再 (exp) 回去,完事。复杂度 (O(nlog n))(但常数巨大)。

多项式带余除法

严格的描述是:给定一个 (n) 次多项式 (F(x)) 和一个 (m) 次多项式 (G(x)) ,请求出多项式 (Q(x),R(x)),满足以下条件:

  • (Q(x)) 次数为 (n−m)(R(x)) 次数小于 (m)
  • (F(x)equiv Q(x)G(x)+R(x)pmod {x^n})

我们记 (F_r(x)=sum_{i=0}^n[x^{n-i}]F(x)x^i),即 (F_r)(F) 系数翻转后的多项式。很容易得到一个等价的式子:

[F_r(x)=x^nF(x^{-1}) ]

开始推式子:

[F(x)=Q(x)G(x)+R(x) \ F(x^{-1})=Q(x^{-1})G(x^{-1})+R(x^{-1}) \ x^nF(x^{-1})=x^{m}Q(x^{-1})x^{n-m}G(x^{-1})+x^nR(x^{-1}) \ F_r(x)=Q_r(x)G_r(x)+x^nR(x^{-1}) ]

式子两边同时 (mod x^{n-m+1}),余数多项式就没了:

[F_r(x)equiv Q_r(x)G_r(x)pmod {x^{n-m+1}} \ Q_r(x)equiv F_r(x)G_r^{-1}(x)pmod {x^{n-m+1}} ]

(Q(x)) 算了出来,(R(x)) 也就可以很容易算出来了。

分治 FFT

给定序列 (g_{1,dots,n-1}),求序列 (f_{0,dots,n-1})。满足 (f_i=sum_{j=1}^if_{i-j}g_j,f_0=1)
发现 (f) 的每一项都依赖于前面的项,考虑用 CDQ 分治解决这个过程。
假设当前已经算出了所有 (f_{lsim mid}) 的值,考虑对 (forall mid+1leq ileq r)(f_i) 所需要加的贡献是 (sum_{j=l}^{mid}f_jg_{i-j}),这刚好是一个卷积的形式。所以我们对于当前区间 ([l,r]),只需将 (f_{lsim mid})(g_{0sim r-l}) 求卷积即可。总复杂度 (O(nlog^2 n))
考虑在模意义下这道题怎么做。对于 ({f_n})({g_n}) 的生成函数 (F(x)=sum_{ige 0}f_ix^i,G(x)=sum_{ige 0}g_ix^i),考虑其卷积

[F(x)G(x)equiv sum_{ige 0}sum_{j=0}^if_{i-j}g_jx^iequiv F(x)-f_0pmod {x^n} \ F(x)equiv frac{1}{1-G(x)}pmod{x^n} ]

直接多项式求逆即可,复杂度 (O(nlog n))竟然比分治快

原文地址:https://www.cnblogs.com/wzzyr24/p/12854613.html