多项式入门

前言

这个傻狗啥都不会还敢来学多项式。

虽然这好像连入门都不算（（（

背板！！！背板！！！

前置的一些定义

点值表示法：将 $n+1$ 个不同的值带入到 $n$ 次多项式中形成 $n+1$ 个点的坐标以代表唯一多项式，支持 $O(n)$ 计算卷积。
单位根：复数意义下 $x^{n}=1$ 即为 $n$ 次单位根。(有 $n$ 个)
本原单位根：$e^{i\frac{2\pi}{n}}$ 称为 $n$ 次单位根的本原单位根，记为 $ω_n$ 。
原根：若 $1\le g\le n-1$ 且 $g$ 模 $n$ 的阶为 $\varphi(n)$ ，则称 $g$ 为 $n$ 的原根。
卷积的一般形式：

c_i=\sum\limits_{i=j\oplus k}a_jb_k

DFT

具体实现

暴力计算点值

b_k=\sum\limits_{i=0}^{n-1} a_{i}k

FFT

实现原理

在通过点值表示法来求卷积的过程中，利用转换和单位根的性质可以将求解的复杂度降为 $O(n\log n)$ 。

具体实现

懒得码那么多又臭又长的证明了，直接给结论：(以下设 $m=\frac{n}{2}$)

令

A(x)=\sum\limits_{i=0}^{n-1} a_ix^i

则

A(\omega_n^{k})=A_0(\omega_m^{k})+\omega_n^{k}A_1(\omega_m^{k})\ (k\le m)

A(\omega_n^{k+m})=A_0(\omega_m^{k})-\omega_n^{k}A_1(\omega_m^{k})\ (k\le m)

其中

A_0=\sum\limits_{i=0}^{m-1} a_{2i}x^i

A_1=\sum\limits_{i=0}^{m-1} a_{2i+1}x^i

就可以递归实现了。

优化(迭代实现)

建出递归树后我们可以发现，最后一层的树是最初的数的二进制下反转后排序的值。

看着很拗口，上例子：

0 1 2 3 4 5 6 7
0 2 4 6, 1 3 5 7
0 4, 2 6, 1 5, 3 7
0, 4, 2, 6, 1, 5, 3, 7
二进制：  
000, 100, 010, 110, 001, 101, 011, 111
翻转后：
000, 001, 010, 011, 100, 101, 110, 111

求出最后一行就可以向上合并了！

这种变换好像学名叫做 $\text{Rader}$ 排序。

$O(n)$ 预处理出 $\text{Rader}$ 排序的映射关系，然后将数组交换就可以 $O(n)$ 求出底层数了！！！

IDFT

具体实现

设

b_k=\sum\limits_{i=0}^{n-1} a_{i}\omega_n^{ik}

即 $k$ 的点值。

则

a_k=\frac{1}{n}\sum\limits_{i=0}^{n-1} b_{i}\omega_n^{-ik}

反正就是 $\text{FFT}$ 的逆变换。

不会证明，摆烂。

IFFT

快速地实现 $\text{IDFT}$ 。

可以发现 $\text{IDFT}$ 和 $\text{DFT}$ 唯一相差的就是 $\omega_m^{-ik}$ ，因此可以和 $\text{FFT}$ 套用同种方法实现。

最后结果记得除以 $n$ 。

NNT

这个 $\text{FFT}$ 又要码又臭又长的复数，又要卡精度，我们考虑用别的东西代替单位根实现。

原根！好东西！！！

再次不加证明地扔出一个结论：

\omega_n\equiv g^{\frac{p-1}{n}}\pmod{p}

证明，狗都不会！！！

然后代换掉 $\text{FFT}$ 中的复数单位根就好了。

Code

$\text{FFT}$ (递归版)

#include<bits/stdc++.h>
#define N 4000005
using namespace std;

inline int read(){
    int x=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){if(ch=='-')w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return x*w;
}

const double phi=acos(-1);
int n,m,lim;
complex<double> f[N],g[N];

void F(complex<double>*a,int n,int opt){
    if(!n) return ;
    complex<double> a0[n],a1[n];
    for(int i=0;i<n;++i)
        a0[i]=a[i<<1],a1[i]=a[i<<1|1];
    F(a0,n>>1,opt);F(a1,n>>1,opt);
    complex<double> t(cos(phi/n),sin(phi/n)*opt),w(1,0);
    for(int i=0;i<n;++i,w*=t)
        a[i]=a0[i]+w*a1[i],a[i+n]=a0[i]-w*a1[i];
}

signed main(){
    n=read();m=read();
    for(int i=0;i<=n;++i) f[i]=read();
    for(int i=0;i<=m;++i) g[i]=read();
    for(lim=1;lim<=n+m;lim<<=1);
    F(f,lim>>1,1);F(g,lim>>1,1);
    for(int i=0;i<lim;++i) f[i]*=g[i];
    F(f,lim>>1,-1);
    for(int i=0;i<=n+m;++i)printf("%.0f ",fabs(f[i].real())/lim);
    return 0;
}

$\text{FFT}$ (迭代版)

#include<bits/stdc++.h>
#define N 2000005
using namespace std;

inline int read(){
    int x=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){if(ch=='-')w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return x*w;
}

const double phi=acos(-1);
int n,m,lim,l;
int r[N];
complex<double> f[N],g[N];

void F(complex<double>*a,int n,int opt){
    for(int i=0;i<lim;++i)
        if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=1;i<lim;i<<=1){
        complex<double> W(cos(phi/i),sin(phi/i)*opt);
        for(int j=0;j<lim;j+=(i<<1)){
            complex<double> w(1,0);
            for(int k=0;k<i;++k,w*=W){
                complex<double> x=a[j+k],y=w*a[j+i+k];
                a[j+k]=x+y; a[j+k+i]=x-y;
            }
        }
    }
}

signed main(){
    n=read();m=read();
    for(int i=0;i<=n;++i) f[i]=read();
    for(int i=0;i<=m;++i) g[i]=read();
    for(lim=1;lim<=n+m;lim<<=1,++l);
    for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    F(f,lim>>1,1);F(g,lim>>1,1);
    for(int i=0;i<lim;++i) f[i]*=g[i];
    F(f,lim>>1,-1);
    for(int i=0;i<=n+m;++i)printf("%.0f ",fabs(f[i].real())/lim);
    return 0;
}

$\text{NNT}$ (只有迭代版)

#include<bits/stdc++.h>
#define N 4000005
#define int long long
using namespace std;

inline int read(){
    int x=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){if(ch=='-')w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return x*w;
}

const int mod = 998244353;//good moshu !!!
const int g1 = 3,g2= 332748118;
int n,m,lim,l;
int r[N];
int f[N],g[N];

int qpow(int x,int k){
    int res=1;
    while(k){
        if(k&1) res=res*x%mod;
        x=x*x%mod;k>>=1;
    }
    return res;
}

void NTT(int *a,int opt){
    for(int i=0;i<lim;++i)
        if(i>r[i]) swap(a[i],a[r[i]]);
    for(int i=1;i<lim;i<<=1){
        int W=qpow(opt==1?g1:g2,(mod-1)/(i<<1));
        for(int j=0;j<lim;j+=(i<<1)){
            int w=1;
            for(int k=0;k<i;++k,w=w*W%mod){
                int x=a[j+k],y=w*a[j+i+k]%mod;
                a[j+k]=(x+y)%mod; a[j+k+i]=(x-y+mod)%mod;
            }
        }
    }
}

signed main(){
    n=read();m=read();
    for(int i=0;i<=n;++i) f[i]=read();
    for(int i=0;i<=m;++i) g[i]=read();
    for(lim=1;lim<=n+m;lim<<=1,++l);
    for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    NTT(f,1);NTT(g,1);
    for(int i=0;i<lim;++i) f[i]*=g[i];
    NTT(f,-1);
    int inv =qpow(lim,mod-2);
    for(int i=0;i<=n+m;++i)printf("%lld ",f[i]*inv%mod);
    return 0;
}

FWT/IFWT

用于处理某一部分特殊的卷积，形如：

c_i=\sum\limits_{i=j| k}a_jb_k

其中的或运算可以改为与/异或。

OR

与FFT大体类似，我们先构造出一个函数，令

F(a)_i=\sum\limits_{j|i=i}a_j

则不难得到：

F(a) \times F(b) = F(c)

那么这是否拥有FFT中点值表示法那样优秀的性质呢，答案是有的，我们可以用类似的方式在 $O(n\log n)$ 的时间复杂度内得到这个生成函数，因为它满足。

F(a)=\text{merge}(F(a_0),F(a_0)+F(a_1))

其中 $a_0$，$a_1$ 的定义与FFT中的不同，单纯指 $a$ 的前后两端。

那么我们可以分治求解(代码比FFT好写很多)。

同理，我们不难推出这个的逆变换：

$$ a=merge(a_0,a_1-a_0) $$

改变符号即可套用同个模板。

And

与或类似的，不难得到：

F(a)=\text{merge}(F(a_0)+F(a_1),F(a_1))

同样可以分治求解(和或的代码差不多)。

逆变换：

$$ a=merge(a_0-a_1,a_1) $$

Xor

这位更是重量级，不能套用上述的方式，要创造一种新的构造方法。

那么我们令：

F(a)_i=\sum\limits_{j\otimes i=0}a_j-\sum\limits_{j\otimes i=1}a_j

不难证明也满足性质：

F(a) \times F(b) = F(c)

同时我们可以证明：

F(a)=\text{merge}(F(a_0)+F(a_1),F(a_1)-F(a_1))

逆推同理。

于是我们就可以打掉模板题)了。

Code

#include<bits/stdc++.h>
#define rep(i,l,r) for(int i=l;i<=r;++i)
#define per(i,r,l) for(int i=r;i>=l;--i)
#define N 200005
#define int long long
using namespace std;

int read(){
	int x=0,w=1; char ch=getchar();
	while(ch>'9'||ch<'0'){if(ch=='-')w=-1;ch=getchar();}
	while(ch>='0'&&ch<='9') x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
	return x*w;
}

const int mod = 998244353;
const int divd = 499122177;

int n,m;
int A[N],B[N],a[N],b[N];

void Mod(int &x,int y){x=(x+y)>mod?(x+y)-mod:x+y;}
void Or(int *f,int opt){
	for(int l=2,k=1;l<=n;l<<=1,k<<=1)
		for(int i=1;i<=n;i+=l) rep(j,0,k-1)
			Mod(f[i+j+k],(opt==1)?f[i+j]:-f[i+j]+mod);
}
void And(int *f,int opt){
	for(int l=2,k=1;l<=n;l<<=1,k<<=1)
		for(int i=1;i<=n;i+=l) rep(j,0,k-1)
			Mod(f[i+j],(opt==1)?f[i+j+k]:-f[i+j+k]+mod);
}
void Xor(int *f,int opt){
	for(int l=2,k=1;l<=n;l<<=1,k<<=1)
		for(int i=1;i<=n;i+=l) rep(j,0,k-1)
			Mod(f[i+j],f[i+j+k]),f[i+j+k]=(f[i+j]-f[i+j+k]*2%mod+mod)%mod,
			f[i+j]=(f[i+j]*opt%mod+mod)%mod, f[i+j+k]=(f[i+j+k]*opt%mod+mod)%mod;
}

signed main(){
	n=read(); n=(1<<n);
	rep(i,1,n) A[i]=a[i]=read(); rep(i,1,n) B[i]=b[i]=read();

	Or(a,1); Or(b,1); rep(i,1,n) a[i]=a[i]*b[i]%mod; Or(a,-1);
	rep(i,1,n) printf("%lld ",a[i]); puts(""); 

	rep(i,1,n) a[i]=A[i],b[i]=B[i];
	And(a,1); And(b,1); rep(i,1,n) a[i]=a[i]*b[i]%mod; And(a,-1);
	rep(i,1,n) printf("%lld ",a[i]); puts("");

	rep(i,1,n) a[i]=A[i],b[i]=B[i];
	Xor(a,1); Xor(b,1); rep(i,1,n) a[i]=a[i]*b[i]%mod; Xor(a,divd);
	rep(i,1,n) printf("%lld ",a[i]); puts("");
	return 0;
}

ChengJY's blog