多项式整合
前言
刚开始学，会慢慢补起来的。
多项式乘法
FFT 快速傅里叶变换
简介
FFT 是高效的多项式乘法算法，其时间复杂度为 \(O(n \log n)\)，但其缺点是算法过程中有多数浮点数运算，容易爆精度。
前置知识
虚数
我们定义虚数单位 \(i\)，满足 \(i^2=-1\)。
自然地，我们定义复数 \(c = a+bi \ (a,b\in \mathbb{R} )\)，以及复数的运算。
设 \(c_1 = a+bi,c_2 = c+di\)，则有：
-
\(c_1+c_2=(a+c)+(b+d)i\)
-
\(c_1-c_2=(a-c)+(b-d)i\)
-
\(c_1c_2=(ac-bd)+(ad+bc)i\)
定义由实轴与虚轴组成的平面为复平面，自然地，我们可以将任意复数 \(c=a+bi\) 看作向量 \((a,b)\)，定义其模长为 \(|c|=\sqrt{a^2+b^2}\)，以正实轴为始边，向量 \((a,b)\) 为终边的角 \(\theta\) 称为复数 \(c\) 的幅角，于是得到复数的三角表示 \(c=|c|(\cos \theta + i\sin \theta )\)。
定义共轭复数 \(\bar c =a -bi\)，一对共轭复数的乘积为实数。
重新观察复数的乘法 \(c_1c_2=(ac-bd)+(ad+bc)i\)，我们有：
-
\(|c_1|=\sqrt{a^2+b^2},|c_2|=\sqrt{c^2+d^2},|c_1c_2|=\sqrt{(ac)^2+(bd)^2+(ad)^2+(bc)^2}=|c_1||c_2|\)，也就是模长相乘。
-
\(\tan \theta_1 = \frac{b}{a},\tan \theta_2 = \frac{d}{c},\tan (\theta_1 + \theta_2)= \frac{bc+ad}{ac-bd}\)，记 \(\theta_0\) 为 \(c_1c_2\) 的幅角，则有 \(\tan \theta_0 = \frac{bd+ad}{ac-bd} = \tan (\theta_1 + \theta_2)\)，也就是 \(\theta_0=\theta_1+\theta_2\)。
也就是说，两复数相乘，模长相乘，幅角相加。
单位根
定义方程 \(x^n=1 \ (n \in \mathbb{N^+})\) 的解为 \(n\) 次单位根，记为 \(\omega_n\)，由代数基本定理可知共有 \(n\) 个不同的 \(n\) 次单位根，考虑其性质：
-
\(\omega_n^{i+n} = \omega_n^i\)。
-
\(w_n^i\) 均为单位根。
所以，所有 \(n\) 次单位根在复平面上是构成正 \(n\) 边形的顶点，其中一个点是 \((1,0)\)。
由此，我们知道：
-
\(\omega_n^i = \omega_{2n}^{2i}\)。
-
若 \(n\) 是偶数，则 \(\omega_n^i = - \omega_n^{i+\frac{n}{2}}\)。
算法
考虑一个多项式的表示方式，在平时使用中，我们大多使用系数表示法，也就是直接存储多项式中每一项的系数，但是我们用此种方法进行多项式乘法时，我们最多做到 \(O(n^2)\) 的时间复杂度，说明，我们需要换一种方法优化。
我们知道任何一个 \((n-1)\) 次多项式，都可以由 \(n\) 个点唯一确定，通过记录 \(n\) 个点及其对应的值来唯一表示一个多项式的方法，我们称之为点值表示法，倘若我们已经知道了这些点值，我们只需要将每个点的取值相乘就行了，这样就会被优化到 \(O(n)\) 的时间复杂度。
现在问题就转化为了取哪些点，如何快速计算点值，以及如何快速将点值表示转化为系数表示。
FFT
以下默认 \(n\) 为 \(2\) 的正整数次幂。
FFT 是快速将系数表示转化为点值表示的算法。
我们将视角移到复数上，我们发现单位根有着非常美妙的性质，于是我们考虑求出所有 \(f(\omega_n^i)\) 的值。
我们不妨设 \(f(x)= a_0+a_1x+a_2x^2+ \cdots +a_{n-1}x^{n-1}\)，我们将每一项按次数奇偶分类，有：
将 \(B(x)\) 中每一项提一个 \(x\) 出来，有：
我们将单位根 \(\omega_n^k\) 带入式子，记 \(m\) 为 \(\frac{n}{2}\)。
- \(0 \le k <m\)，讨论 \(\omega_n^k\)：
- \(0 \le k <m\)，讨论 \(\omega_n^{k+m}\)：
我们发现，如果我们已经预处理出了 \(A(\omega_m^k),B(\omega_m^k)\)，那么我们可以 \(O(1)\) 求出 \(f(\omega_n^k),f(\omega_n^{k+m})\)，我们发现 \(A,B\) 的次数减半了，于是我们可以递归求解，时间复杂度 \(O(n \log n)\)。
IFFT
IFFT 是将点值表达转化为系数表达的算法。
我们现在有 \(f(x)= a_0+a_1x+a_2x^2+ \cdots +a_{n-1}x^{n-1}\) 的不同单位根的函数值 \(b_0,b_1,\cdots,b_{n-1}\)，我们有式子：
那么有结论：
将前一个式子代入即可证明，我们发现实际上就是代入 \(\omega_n^{-k}\) 进 FFT 去计算，本质与 FFT 相同。
代码实现
int n,m;
const double pi=3.141592653589793;
complex<double>f[4000005],g[4000005];
void FFT(complex<double> *a,int n,int opt){
if(n==1)return;
int m=n>>1;
complex<double> al[m+5],ar[m+5];
for(int i=0;i<m;i++)al[i]=a[i<<1],ar[i]=a[i<<1|1];
FFT(al,m,opt),FFT(ar,m,opt);
complex<double> W(cos(pi/m),sin(pi/m)*opt),w(1,0);
for(int i=0;i<m;i++,w*=W){
a[i]=al[i]+w*ar[i],a[i+m]=al[i]-w*ar[i];
}
}
signed main(){
n=read(),m=read();
for(int i=0;i<=n;i++)f[i].real(read());
for(int i=0;i<=m;i++)g[i].real(read());
for(m+=n,n=1;n<=m;n<<=1);
FFT(f,n,1),FFT(g,n,1);
for(int i=0;i<n;i++)f[i]*=g[i];
FFT(f,n,-1);
for(int i=0;i<=m;i++){
printf("%.0lf ",fabs(f[i].real())/n);
}
return 0;
}
NTT 快速数论变换
NTT 主要适用于值域较小或特定模数（如 \(998244353,469762049,1004535809\)）的情况使用，优点是没有精度误差，缺点是只能在特定模数下使用。
前置知识
原根
此处知识我其实还不是很熟。
定义（阶）：如果 \((a,m)=1\)，存在最小的 \(n\) 使得 \(a^n \equiv1 \ (mod \ m)\)，则称 \(n\) 是 \(a\) 模 \(m\) 的阶，记为 \(\delta_m (a)\)。
由于欧拉定理可知 \(a^{\varphi(m)} \equiv 1 \ (mod\ m)\)，所以一定存在一个阶。
定义（原根）：若 \((a,m)=1,\delta_m (a)=\varphi(m)\)，则称 \(a\) 是模 \(m\) 的原根。
（性质）模 \(m\) 意义下，\(a^1,a^2,\cdots,a^{\delta_m(a)}\) 两两不同余。
显然，如果 \(a^i \equiv a^j \ (mod \ m)\)，则 \(a^{|i-j|} \equiv 1 \ (mod\ m)\)，矛盾。
NTT
其实就是 FFT 将虚数单位根换成了模意义下的原根，在 \(m=998244353\) 时，原根 \(g=3\)，那么模 \(m\) 意义下的 \(n\) 次单位根即为 \(g^\frac{m-1}{n}\)。
而为什么需要、像 \(998244353\) 这样的特定模数呢？
因为，\(998244353=2^{23} \times 7 \times 17 +1\)，由于 \(n\) 是 \(2\) 的正整数次幂，所以 \(n \le 2^{23}\) 范围足够大。
而且 \(998244353\) 的原根好找啊，所以为了方便，一般模数都取 \(998244353\)，当然，也可以通过中国剩余定理完成任意模数 NTT。
多项式初等函数
多项式乘法求逆
算法
多项式乘法的逆元：现有 \(n-1\) 次多项式 \(f(x)\)，若函数 \(g(x)\) 满足 \(f(x)g(x) \equiv 1 \pmod {x^n}\)，则 \(g(x)\) 是 \(f(x)\) 的逆元。
这里讨论 \(\bmod \ x^n\) 下的逆元。
求解逆元的算法主要依靠倍增实现，我们现在要求 \(mod \ x^n\) 意义下的逆元，若现在我们已经得出了 \(mod \ x^{\frac{n}{2}}\) 的逆元，考虑如何推导。
设 \(f(x)g(x) \equiv 1 \pmod {x^n}\)，\(f(x)g'(x) \equiv 1 \pmod {x^{\frac{n}{2}}}\)，则：
两边同时平方：
同时乘上 \(f(x)\)，有：
移项即可得到递推式：
利用 NTT 快速多项式乘法即可，主定理分析即可得到时间复杂度 \(O(n\log n)\)。
Code
#include <bits/stdc++.h>
#define ll long long
#define ull unsigned long long
using namespace std;
inline int read(){
int ans=0,w=1;
char ch=getchar();
while(ch<48||ch>57){
if(ch=='-')w=-1;
ch=getchar();
}
while(ch>=48&&ch<=57){
ans=(ans<<1)+(ans<<3)+ch-48;
ch=getchar();
}
return w*ans;
}
const int mod=998244353,gn=3,inv=332748118;
ll Pow(ll x,ll y){
ll ans=1,push=x;
while(y){
if(y&1)ans=ans*push%mod;
push=push*push%mod;
y>>=1;
}
return ans;
}
int n;
int a[400005],p[400005];
int f[400005],g[400005],d1[400005],d2[400005],res[400005];
int rev[400005];
inline void init(int k){
for(int i=0;i<(1<<k);i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
}
inline void NTT(int *a,int *r,int n,int opt){
for(int i=0;i<n;i++)r[i]=a[i];
for(int i=0;i<n;i++){
if(i<rev[i])swap(r[i],r[rev[i]]);
}
for(int l=2;l<=n;l<<=1){
int W=Pow((opt==1)?gn:inv,(mod-1)/l);
for(int i=0;i<n;i+=l){
int w=1;
for(int j=i;j<i+(l>>1);j++,w=1ll*w*W%mod){
int p=r[j],q=r[j+(l>>1)];
r[j]=(p+1ll*w*q)%mod,r[j+(l>>1)]=(p-1ll*w*q%mod+mod)%mod;
}
}
}
if(opt==-1)for(int i=0;i<n;i++)r[i]=1ll*r[i]*Pow(n,mod-2)%mod;
}
inline void mul(int n,int *f,int *g,int *res){
init(__lg(n));
NTT(f,d1,n,1),NTT(g,d2,n,1);
for(int i=0;i<n;i++)d1[i]=1ll*d1[i]*d2[i]%mod;
NTT(d1,res,n,-1);
for(int i=0;i<n;i++)d1[i]=d2[i]=0;
}
int b[400005],c[400005];
void get_inv(int n,int *f,int *res){
memset(f,0,sizeof f);
int len=1,m=4;
f[0]=a[0];
res[0]=Pow(f[0],mod-2);
while(len<n){
memset(b,0,sizeof b);
mul(m,res,res,b);
for(int i=len;i<2*len;i++)f[i]=a[i];
mul(m,b,f,b);
for(int i=len;i<2*len;i++)res[i]=(2ll*res[i]-b[i]+mod)%mod;
len<<=1;m<<=1;
}
}
signed main(){
n=read();
for(int i=0;i<n;i++)a[i]=read();
get_inv(n,p,res);
for(int i=0;i<n;i++)cout<<res[i]<<' ';
return 0;
}
浙公网安备 33010602011771号