多项式整合包

update:

学习多项式!!

1.快速傅里叶变换(FFT)

当然是从 $F F T$ 开始了:)

1.1 DFT & IDFT

1.1.1 系数与点值

系数表达: $F (x) = \sum_{i = 0}^{n} F [i] x^{i}$

$F [i]$ 叫做多项式 $F (x)$ 的系数。
这样我们求 $C (x) = F (x) * G (x)$ 的系数:

C [k] = \sum_{i = 0}^{k} F [i] G [k - i]

直接算的复杂度显然是 $O (n^{2})$ 的。

然后我们要知道: n+1个点值(有序数对)可以来确定一个n次多项式 $F (x)$ (具体的可以看 拉格朗日插值) 。
所以多项式的 点值表达 就是用 $n + 1$ 的点值来表示多项式。
点值表达： $F (x) = (x_{0}, y_{0}), (x_{1}, y_{1}), . . ., (x_{n}, y_{n})$ 。

我们利用这样的表示法可以快速的求出两个 多项式的乘法。
只需将 $F (x)$ 与 $G (x)$ 的 点值对应值直接乘 即可。

但是 $F (x) * G (x)$ 明明有 $2 n$ 项，怎么办。

我们直接找 $2 n + 1$ 个点就可以了。

这样我们就完成了 $F (x) * G (x)$ 的点值表达: $(x_{0}, F_{x_{0}} * G_{x_{0}}), (x_{1}, F_{x_{1}} * G_{x_{1}}), . . ., (x_{2 n}, F_{x_{2 n}} * G_{x_{2 n}})$ 。
复杂度仅需 $O (2 n)$ 。

1.1.2 FFT的流程

上述我们可以知道在多项式乘法中 点值表达 的复杂度很优，但是我们一般的多项式都是系数表达 ，我们如何实现点值与系数之间的转换呢？: )

“把系数表达转换为点值表达”的算法叫做 DFT

“把点值表达转换为系数表达”的算法叫做 IDFT(DFT的逆运算)

传承老图：

从一个多项式的系数表达确定其点值表达的过程称为求值
而求值运算的逆运算(也就是从一个多项式的点值表达确定其系数表达)被称为插值.

1.2 单位根的性质

这一点讲如何利用一些有神奇性质的 $x$ 。

复数不讲。

n次单位根： 方程 $x^{n} = 1$ 的 $n$ 个解。
几何意义：将单位圆均分为 $n$ 份，我们将幅角为 $\frac{k}{n}$ 长度为 $1$ 的根称为 $ω_{n}^{k}$
复数表示： $ω_{n}^{k} = c o s \frac{2 k π}{n} + s i n \frac{2 k π}{n} i$

单位根的性质：

$ω_{n}^{j} * ω_{n}^{i} = ω_{n}^{i + j}$
$ω_{2 n}^{2 k} = ω_{n}^{k}$ ，根据几何或者直接表示为复数易证。
$ω_{n}^{k + n} = ω_{n}^{k}$ ，几何上相当于绕了一个圆周
$ω_{n}^{k + \frac{n}{2}} = ω_{n}^{k + \frac{n}{2}} = ω_{n}^{k} * ω_{n}^{\frac{n}{2}} = - ω_{n}^{k}$ ，几何上相当于转了半个圆周

这些性质有大用

1.3 单位根的巧用

一个多项式 $F (x) = \sum_{i = 0}^{n} F [i] x^{i}$ 。
我们把它拆成平均的两份：

F^{0} (x) = F [0] + F [2] x + . . . + F [n - 1] x^{\frac{n}{2} - 1}

F^{1} (x) = F [1] + F [3] x + . . . + F [n] x^{\frac{n}{2} - 1}

则可以得到 $F (x) = F^{0} (x^{2}) + x F^{1} (x^{2})$
这样就可以分治了！！
我们把神奇的单位根带入进去：

$1. k > \frac{n}{2}$ 代入 $ω_{n}^{k}$
$F (ω_{n}^{k}) = F^{0} ((ω_{n}^{k})^{2}) + ω_{n}^{k} F^{1} ((ω_{n}^{k})^{2}) = F^{0} (ω_{\frac{n}{2}}^{k}) + ω_{n}^{k} F^{1} (ω_{\frac{n}{2}}^{k})$
$2. k < \frac{n}{2}$ 代入 $ω_{n}^{k + \frac{n}{2}}$
$F (ω_{n}^{k + \frac{n}{2}}) = F^{0} ((ω_{n}^{k + \frac{n}{2}})^{2}) + ω_{n}^{k + \frac{n}{2}} F^{1} ((ω_{n}^{k + \frac{n}{2}})^{2}) = F^{0} (ω_{\frac{n}{2}}^{k}) - ω_{n}^{k} F^{1} (ω_{\frac{n}{2}}^{k})$

两式的差别只是正负号。

有什么用呢？可以快速求值，我们利用 $1$ 式可以用 $F^{0} (x)$ 与 $F^{1} (x) O (n)$ 求出 $ω_{n}^{0}, ω_{n}^{1}, . . ., ω_{n}^{\frac{n}{2} - 1}$ 的点值，类似利用 $2$ 式的可以求出 $ω_{n}^{\frac{n}{2}}, . . ., ω_{n}^{n - 1}$ 的点值。

最终得出来了 $F (x)$ 的点值！这是一个分治过程，即可以 $O (n l o g n)$ 的复杂度完成求值 (即DFT)

至于IDFT。
结论：只需带入 $ω_{n}^{- k}$ 最后除 $n$ 即可，证明先咕。

1.4 FFT 的实现

上面听不懂没关系，背代码就行了:(

我们上面的方法是递归分治的，会造成大量数组拷贝，是不优的。
有一种方法可以避免，就是迭代FFT(即从小合并)。
但这样我们需要找出每个点最后所在的位置，结论就是原数字的二进制翻转(很奇怪)，可以 $O (n)$ 求。

for(int i = 1;i < k;i++)rev[i] = (rev[i>>1]>>1) | ((i&1)<<(b-1));

剩下的直接看代码吧:)

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define db double
const int N = 2e6+1e5;//2^20 = 1048576
const db pi = acos(-1.0);

inline int read(){
    int x = 0,f = 1;char c = getchar();
    for(;c < '0' || c > '9';c = getchar())if(c == '-')f = -1;
    for(;c >= '0' && c <= '9';c = getchar())x = (x<<1) + (x<<3) + c-'0';
    return x * f;
}

struct cp{
    db x,y;
    cp(){x = y = 0;}
    cp(db x,db y):x(x),y(y){}
    inline cp operator + (const cp a)const{return cp(a.x+x,a.y+y);}
    inline cp operator - (const cp a)const{return cp(x-a.x,y-a.y);}
    inline cp operator * (const cp a)const{return cp(x*a.x-y*a.y,x*a.y+y*a.x);}
}F[N],G[N];

int n,m;
int rev[N],b;

void FFT(int n,cp *a,int op){
    for(int i = 0;i < n;i++)if(i < rev[i])swap(a[i],a[rev[i]]);
    for(int len = 1;len <= (n>>1);len <<= 1){
        cp w1(cos(pi/len),sin(pi/len)*op);
        for(int i = 0;i <= n - (len<<1);i += (len<<1)){
            cp w(1,0);
            for(int j = i;j < i+len;j++){
                cp x = a[j],y = w*a[j+len];
                a[j] = x+y,a[j+len] = x-y;
                w = w * w1;
            }
        }
    }
}
int main(){
    n = read(),m = read();
    for(int i = 0;i <= n;i++)F[i].x = read();
    for(int i = 0;i <= m;i++)G[i].x = read();
    int k = 1;
    while(k <= n+m)k <<= 1,b++;
    for(int i = 1;i < k;i++)rev[i] = (rev[i>>1]>>1) | ((i&1)<<(b-1));
    FFT(k,F,1),FFT(k,G,1);
    for(int i = 0;i <= k;i++)F[i] = F[i] * G[i];
    FFT(k,F,-1);
    for(int i = 0;i <= n+m;i++)printf("%d ",(int)(F[i].x/k+0.5));

	return 0;

}

全文背默！！！

1.5 其他优化

咕咕咕

1.6 快速数论变换(NTT)

学原根回来哩:)
因为FFT引用了复数单位根，所以精度要求很高，跑的可能会慢。
所以我们需要NTT。

1.6.1 阶与原根

详见数论学习笔记1.3
我就不搬了。

1.6.2 原根类似的性质

(以下 $g$ 代表模 $p$ 的最小原根)
由原根性质2可得 $g^{1}, g^{2}, . . ., g^{φ (p)}$ 模 $p$ 两两不同余。

我们令 $ω_{n} = g^{\frac{p - 1}{n}}$ ，你问不整除怎么办，我们就找能整除的 $p$ 。

所以我们有以下性质：

$ω_{n}^{n} = (g^{\frac{p - 1}{n}})^{n} = g^{p - 1} \equiv 1 (\mod p)$
$ω_{2 n}^{2 k} = (g^{\frac{p - 1}{2 n}})^{2 k} = (g^{\frac{p - 1}{n}})^{k} = ω_{n}^{k}$
$ω_{n}^{k + n} = (g^{\frac{p - 1}{n}})^{k + n} \equiv (g^{\frac{p - 1}{n}})^{k} = ω_{n}^{k} (\mod p)$
$ω_{n}^{k + \frac{n}{2}} = (g^{\frac{p - 1}{n}})^{k + \frac{n}{2}} \equiv - (g^{\frac{p - 1}{n}})^{k} = - ω_{n}^{k} (\mod p)$

第四条证明一下：

由费马小定理知 $g^{p - 1} - 1 \equiv 0 (\mod p)$ ，即 $(g^{\frac{p - 1}{2}} - 1) (g^{\frac{p - 1}{2}} + 1) \equiv 0 (\mod p)$ ，即 $g^{\frac{p - 1}{2}} \equiv \pm 1$ ，又因为 $g^{i}$ 互不相同且 $g^{φ (p)} \equiv 1 (\mod p)$ ，得出 $g^{\frac{p - 1}{2}} \equiv - 1 (\mod p)$ ，即得证。

这样原根与单位根的性质就完全相同了！！: )

我们只需要找到一个模数 $p$ 使得任意 $n$ 都整除 $p - 1$ 即可，而一般 $n$ 都是 $2$ 的次幂。

一般NTT模数：
$998244353 = 119 * 2^{23} + 1$ ， $g$ 为 $3$ 。
$1004535809 = 479 * 2^{21} + 1$ ， $g$ 为 $3$ 。

大佬的NTT模数表

模版

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define F(i,a,b) for(int i = a;i <= b;i++)
#define D(i,a,b) for(int i = a;i >= b;i--)
#define in inline
#define pb(x) push_back(x)
#define db double  
const db pi = acos(-1);
const int N = 2e6+1e5,mod = 998244353,g = 3;

ll read(){
    ll x = 0,f = 1;char c = getchar();
    for(;c < '0' || c > '9';c = getchar())if(c == '-')f = -1;
    for(;c >= '0' && c <= '9';c = getchar())x = (x<<3) + (x<<1) + c-'0';
    return x * f;
}

int ksm(int x,int y = mod-2){
    int ans = 1;
    while(y){
        if(y & 1)ans = 1ll * ans * x % mod;
        x = 1ll * x * x % mod,y >>= 1;
    }return ans;
}


int b,n,m;
int rev[N],invg = ksm(g),invn;
ll F[N],G[N];

void NTT(int n,ll *a,int op){
    F(i,0,n)if(i < rev[i])swap(a[i],a[rev[i]]);
    for(int len = 1;len <= (n>>1);len <<= 1){
        ll tg = ksm(op?g:invg,(mod-1)/(len<<1));
        for(int i = 0;i + (len<<1) <= n;i += (len<<1)){
            ll gg = 1;
            for(int j = i;j < i + len;j++){
                ll x = a[j],y = 1ll * gg * a[j+len] % mod;
                a[j] = (x + y) % mod,a[j+len] = (x - y + mod) % mod;
                gg = 1ll * gg * tg % mod;
            }
        }
    }
}
int main(){
    n = read(),m = read();
    F(i,0,n)F[i] = read();
    F(i,0,m)G[i] = read();
    int t = 1;
    while(t <= n+m)t <<= 1,b++;
    F(i,1,t)rev[i] = ((rev[i>>1]>>1) | (i&1)<<(b-1));
    NTT(t,F,1),NTT(t,G,1);
    F(i,0,t)F[i] = F[i] * G[i] % mod;
    NTT(t,F,0);invn = ksm(t);
    F(i,0,n+m)printf("%d ",int(F[i] * invn % mod));
    printf("\n");
    

    return cerr<<endl<<"Time:"<<clock()<<"ms"<<endl,0;

}