多项式半家桶模板
多项式半家桶
多项式乘法
考虑两个多项式 $ F(x)=\sum_{i=0}^{i<n} A_i x_i \ \ , G(x)=\sum_{i=0}^{i<m} B_i x_i $ ,求 $ F(x) \times G(x) $。
暴力做,复杂度为 $ O(NM) $ ,不够优美。我们有 $ O((N+M)\log(N+M)) $ 的做法。
一个 \(n\) 次多项式可以由 \(n+1\) 个互不相同的点值表示,即点值表示法,取 \(N+M\) 个点的点值分别得到两个多项式的点值表示法,那么两个多项式的乘积的点值表示法就是两者直接按位相乘,这一部分复杂度为 $ O(N+M) $ 。
那么这里的瓶颈在于快速求 $ F(x),G(x) $ 的点值和将点值转化为多项式。
快速傅里叶变换(FFT)
不妨认为 \(n\) 为2的整数次幂,不符的话就向高位补零。
前置知识:单位根
设 $ x^n=1 $ ,则在复数域下 \(x\) 有 \(n\) 个解,第 \(k\) 个解记为第 \(k\) 个 \(n\) 次单位根 $ w_n^k $ ,其值为 $ e^{i\frac{2k\pi}{n}} = \cos(\frac{2k\pi}{n}) + i \sin(\frac{2k\pi}{n}) $ 。可以在单位圆上直观的表示。
有如下性质:
- $ w^{k+\frac{n}{2}}_n = -w_n^k $
- $ ( w_{2n}^k )^2 = w_{2n}^{2k} = w_n^k $
- $ w_{\frac{n}{2}}^k = w_n^{2k} $
对于多项式 $ F(x)=\sum_{i=0}^{i<n} A_i x_i $ 按奇偶性分开,变成 $ F_0(x) = \sum_{i=0}^{i<\frac{n}{2}} A_{2i} x^{2i} , F_1(x) = x \times \sum_{i=0}^{i<\frac{n}{2}} A_{2i+1} x^{2i} $ ,进行分治求解,最后再合并。
假设我们已经得到 $ F_0,F_1 $ 的点值表示法。
对于 $ k < \frac{n}{2} $ ,有
这样就可以类似于线段树一样分治了。最后 $ F(x) $ 的第 \(i\) 项即为 $ x = w_n^i $ 对应的点值。
快速傅里叶逆变换
接下来考虑把点值转化成多项式的一般形式。
设 $ y_i = F(w_n^i) $ , 定义 $ c_k = \sum_{i=0}^{ i<n } y_i (w_n^{ -k })^i $
容易发现,当 $ j \neq k $ 时,有 $ \sum_{i=0}^{i<n} w_n^{i(j-k)} = 0 $ ,因此
由上述定义可知,$ c_k $ 就是多项式 $ H(x) = \sum_{i=0}^{i<n} y_i x^i $ 在 $ w_n^{-k} $ 时的点值。
所以只要把得到的点值再做一遍FFT即可。
递归版实现
CODE
struct comp{
double x,y;
comp operator + (comp k) {return {x+k.x,y+k.y};}
comp operator - (comp k) {return {x-k.x,y-k.y};}
comp operator * (comp k) {return {x*k.x-y*k.y,x*k.y+y*k.x};}
}f[N],g[N],h[N];
void FFT(comp *x,int len,int tp){
if(len==1)return ;
for(int i=0;i<(len>>1);i++)h[i]=x[i<<1],h[i+(len>>1)]=x[i<<1|1];
for(int i=0;i<len;i++)x[i]=h[i];
FFT(x,len>>1,tp);
FFT(x+(len>>1),len>>1,tp);
comp wn={cos(pi*2/len),sin(pi*2/len)*tp},w={1,0};
for(int i=0;i<len;i++)h[i]=x[i];
for(int i=0;i<(len>>1);i++,w=w*wn){
comp t=w*h[i+(len>>1)];
x[i]=h[i]+t;
x[i+(len>>1)]=h[i]-t;
}
}
while(len<=n+m)len<<=1;
FFT(f,len,1);
FFT(g,len,1);
for(int i=0;i<len;i++)f[i]=f[i]*g[i];
FFT(f,len,-1);
for(int i=0;i<=n+m;i++)printf("%d ",int(f[i].x/len+0.5));
注意数组大小和精度误差。千万不要忘记除多项式长度
迭代版实现
从底层往上层迭代,复杂度不变,但是可以有效减小常数。
观察一下可以发现,第 \(i\) 项多项式系数在底层的位置是它的二进制翻转。(卡常以后再说)
CODE
void init(){
int p=0,d=(len>>1);
tax[p++]=0; tax[p++]=d;
for(int i=2;i<len;i<<=1){
d>>=1;
for(int j=0;j<i;j++)tax[p++]=(tax[j]|d);
}
}
void FFT(comp *x,int tp){
for(int i=0;i<len;i++){
if(tax[i]>i)swap(x[i],x[tax[i]]);
}
for(int i=2;i<=len;i<<=1){
for(int j=i;j<=len;j+=i){
comp wn={cos(pi*2/i),sin(pi*2/i)*tp},w={1,0};
for(int k=0;k<(i>>1);k++,w=w*wn){
comp t=w*x[k+(i>>1)+j-i];
h[k+j-i]=x[k+j-i]+t;
h[k+(i>>1)+j-i]=x[k+j-i]-t;
}
for(int k=j-i;k<j;k++)x[k]=h[k];
}
}
}
快速数论变换(NTT)
这个似乎更常用一点。
前置知识:原根
满足 $ a^{n} \equiv 1 \ \ (\mod m) $ 的最小正整数 \(n\) ,称为 \(a\) 模 \(m\) 的阶,记作 $ \delta_m(a) $ 。
已知欧拉定理:$ gcd(a,m) = 1 , a^{\varphi(m)} \equiv 1 \ \ (\mod m) $ 。
正整数 \(a\) 是 \(m\) 的原根,当且仅当 $ \delta_m(a) = \varphi(m) $
主要用到这一性质:$ a^1 , a^2 , \dots , a^{ \delta_m(a) } $ 在模 \(m\) 的意义下互不相同。(具体证明可以看 oi-wiki )
建议掌握一下求原根的方法。
如果你已经学会FFT,那NTT就很简单,直接拿原根替换单位根即可。避免了精度误差而且常数更小。
模数一般取 \(998244353\) ,最小的原根是 \(3\) 。(稍后会说明为什么常用)
CODE
namespace Poly{
int len,tax[N],L;
ll h[N],dw[2][N];
void init(){
for(int i=0;i<len;i++)
tax[i]=((tax[i>>1]>>1)|((i&1)<<(L-1)));
dw[1][1]=dw[0][1]=1;
for(int i=2;i<len;i<<=1){
ll w1=qpow(G,(mod-1)/i/2),w0=qpow(Gi,(mod-1)/i/2);//?
dw[1][i]=1;dw[0][i]=1;
for(int j=1;j<i;j++){
dw[1][i+j]=dw[1][i+j-1]*w1%mod;
dw[0][i+j]=dw[0][i+j-1]*w0%mod;
}
}
}
void NTT(ll *x,int tp){
static ull pa[N];
for(int i=0;i<len;i++)pa[i]=x[tax[i]];
for(int i=2;i<=len;i<<=1){
if((i>>18)&1){
for(int j=0;j<len;j++)pa[j]%=mod;
}
for(int j=i,bg=0;j<=len;j+=i,bg+=i){
for(int k=bg;k<bg+(i>>1);k++){
ull ch=pa[k+(i>>1)]*dw[tp][(i>>1)+k-bg]%mod;
pa[k+(i>>1)]=pa[k]+mod-ch;
pa[k]+=ch;
}
}
}
for(int i=0;i<len;i++)x[i]=pa[i]%mod;
}
void Mul(ll *aas,ll *fir,int lx,ll *sec,int ly){
static ll pa[N],pb[N];
memcpy(pa,fir,sizeof(ll)*lx);
memcpy(pb,sec,sizeof(ll)*ly);
len=1;L=0;
while(len<lx+ly-1)len<<=1,L++;
memset(pa+lx,0,sizeof(ll)*(len-lx));
memset(pb+ly,0,sizeof(ll)*(len-ly));
init();
NTT(pa,1);NTT(pb,1);
for(int i=0;i<len;i++)aas[i]=pa[i]*pb[i]%mod;
NTT(aas,0);
ll k=qpow(len,mod-2);
for(int i=0;i<len;i++)aas[i]=aas[i]*k%mod;
}
}
注意到代码中打问号的一行,不知道大家有没有疑惑,如果不能整除怎么办。然而事实是,如果真不能整除,那就得换模数。但一般不会这样,因为我们惊奇地发现:$ 998244352 = 2^{23} \times 7 \times 17 $ ,大多数数据都不会卡满所有2因子。
多项式乘法逆
给定 $ G(x) $ ,求 \(F(x)\) 使得 $ F(x) G(x) \equiv 1 \ \ (\mod x^n) $ 。保证 $ [x^0] G(x) \neq 0 $
倍增求解,不妨认为我们已经求出模 $ x^{\lceil \frac{n}{2} \rceil} $ 意义下的乘法逆。
设 $ H(x) G(x) \equiv 1 \ \ (\mod x^{\lceil \frac{n}{2} \rceil}) $
可以用牛顿迭代推
多项式开根
给定 $ G(x) $ ,求 $ F(x) $ 使得 $ F^2(x) \equiv G(x) \mod x^n $
设 $ H^2(x) \equiv G(x) \mod x^{\lceil \frac{n}{2} \rceil} $
需要用到多项式求逆。
可以用牛顿迭代推
多项式对数函数(多项式ln)
求 $ F(x) $ 使得 $ F(x) \equiv \ln G(x) \mod x^n $
前置知识:基础微积分
只能讲一些基本的概念和用法,建议找几篇博客学习,效果可能更好。
微分,在我的理解中,就是对可导函数求导,表示的是函数的瞬时变化率。(如果你会导数就可以跳过了)
举个例子,设 $ f(x) = x^2 $ , \(dx\) 为一极小量, $ f(x) $ 的导数 $ f'(x) = \frac{ ( x+dx )^2 - x^2 }{dx} = \frac{ 2 x \ dx + dx^2 }{dx} = \lim_{ dx \to 0 } ( 2 x + dx ) = 2x $ 。同理,若 $ f(x) = x^3 $ ,则有 $ f'(x) = 3 x^2 $
给出一些常见函数的导数:
$ ( k x )' = k $
$ ( x^a )' = a \times x^{ a-1 } $
$ ( \ln x )' = \frac1{x} $
$ ( a^x )' = a^x \times \ln a $
导数符合以下运算法则:
- $ ( f(x) + g(x) )' = f'(x) + g'(x) $
- $ ( f(x) \ g(x) )' = f'(x) g(x) + f(x) g'(x) $
- $ f( g(x) )' = g'(x) \ f'( g(x) ) $
至于积分,就是微分的逆运算,在函数图像上理解为面积。
对于一个多项式 $ F(x) = \sum A_i x_i $ 来说, $ F'(x) = \sum i A_i x_{i-1} $ 。对其逆向运算即可得到原多项式。
已知 $ F(x) \equiv \ln G(x) \mod x^n $ ,求导得到
最后再对 $ F'(x) $ 积分一下就做完了。
可以用牛顿迭代推
多项式指数函数(多项式exp)
前置知识:多项式牛顿迭代
首先介绍泰勒展开:
其中 $ f^{(i)} $ 为f的 \(i\) 阶导数,上式称为 $ f(x) $ 在 $ x_0 $ 处的泰勒展开。
现在考虑如下问题:求一个多项式 $ f(x) $ ,使得函数 $ g(f(x)) \equiv 0 \mod x^n $
我们倍增地去解决这个问题,设已经得到模 $ x^{ \frac{n}{2} } $ 的解 $ f_0 $ ,对 $ g(f(x)) $ 在 $ f_0 $ 处泰勒展开,得到 $$ g(f(x)) = \sum_{i=0} \frac{ g^{(i)}(f_0) }{i!} ( f(x)-f_0 )^i \mod x^n $$
假设存在一个 $ h(x) $ 满足 $ f(x) = f_0 + x^{\frac{n}{2}} \cdot h(x) $ 。
代入得到:
注意到,当 $ i \ge 2 $ 时 ,$ ( x^{\frac{n}{2}} \cdot h(x) )^i \equiv 0 \mod x^n $ ,化简可知
代入 $ f(x) = f_0 + x^{\frac{n}{2}} \cdot h(x) $
至于 $ h(x) $ 为什么总是存在,先感受一下。
求 $ F(x) \equiv e^{ G(x) } \mod x^n $
代入多项式牛顿迭代
好玩吧
模板
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N=8e5+100;
const ll mod=998244353;
const ll G=3,Gi=332748118;
int m,n;
inline int read(){
int k=0;
char ch=getchar();
while(ch<'0'||ch>'9')ch=getchar();
while(ch>='0'&&ch<='9'){
k=k*10+ch-'0';
ch=getchar();
}
return k;
}
inline void write(int x){
if(!x)return putchar('0'),void();
int tmp[20],len=0;
while(x)tmp[++len]=x%10,x/=10;
for(int i=len;i>=1;i--)putchar(char(tmp[i]+'0'));
}
inline ll qpow(ll x,ll y){
ll res=1;
while(y){
if(y&1)res=res*x%mod;
x=x*x%mod;
y>>=1;
}
return res;
}
ll f[N],t[N],ny[N];
namespace Poly{
int len,tax[N],L;
ll h[N],dw[2][N];
void init(){
for(int i=0;i<len;i++)
tax[i]=((tax[i>>1]>>1)|((i&1)<<(L-1)));
dw[1][1]=dw[0][1]=1;
for(int i=2;i<len;i<<=1){
ll w1=qpow(G,(mod-1)/i/2),w0=qpow(Gi,(mod-1)/i/2);
dw[1][i]=1;dw[0][i]=1;
for(int j=1;j<i;j++){
dw[1][i+j]=dw[1][i+j-1]*w1%mod;
dw[0][i+j]=dw[0][i+j-1]*w0%mod;
}
}
}
void NTT(ll *x,int tp){
static ull pa[N];
for(int i=0;i<len;i++)pa[i]=x[tax[i]];
for(int i=2;i<=len;i<<=1){
if((i>>18)&1){
for(int j=0;j<len;j++)pa[j]%=mod;
}
for(int j=i,bg=0;j<=len;j+=i,bg+=i){
for(int k=bg;k<bg+(i>>1);k++){
ull ch=pa[k+(i>>1)]*dw[tp][(i>>1)+k-bg]%mod;
pa[k+(i>>1)]=pa[k]+mod-ch;
pa[k]+=ch;
}
}
}
for(int i=0;i<len;i++)x[i]=pa[i]%mod;
}
void Mul(ll *aas,ll *fir,int lx,ll *sec,int ly){
static ll pa[N],pb[N];
memcpy(pa,fir,sizeof(ll)*lx);
memcpy(pb,sec,sizeof(ll)*ly);
len=1;L=0;
while(len<lx+ly-1)len<<=1,L++;
memset(pa+lx,0,sizeof(ll)*(len-lx));
memset(pb+ly,0,sizeof(ll)*(len-ly));
init();
NTT(pa,1);NTT(pb,1);
for(int i=0;i<len;i++)aas[i]=pa[i]*pb[i]%mod;
NTT(aas,0);
ll k=qpow(len,mod-2);
for(int i=0;i<len;i++)aas[i]=aas[i]*k%mod;
}
inline void Lim(ll *x,int lx,int ln){
if(lx>ln)memset(x+ln,0,sizeof(ll)*(lx-ln));
}
inline void Add(ll *x,ll *y,int ln){
for(int i=0;i<ln;i++)x[i]=(x[i]+y[i])%mod;
}
inline void Red(ll *x,ll *y,int ln){
for(int i=0;i<ln;i++)x[i]=(x[i]+mod-y[i])%mod;
}
inline void mul(ll *x,ll k,int ln){
for(int i=0;i<ln;i++)x[i]=x[i]*k%mod;
}
inline void Der(ll *x,ll *y,int ln){
for(int i=1;i<ln;i++)x[i-1]=y[i]*i%mod;
}
inline void Int(ll *x,ll *y,int ln){
for(int i=ln;i>0;i--)x[i]=y[i-1]*ny[i]%mod;
x[0]=0;
}
void Inv(ll *aas,ll *x,int lx){
assert(x[0]);
aas[0]=qpow(x[0],mod-2);
static ll pa[N],pb[N];
memset(pa,0,sizeof(pa));memset(pb,0,sizeof(pb));
for(int i=1;i<lx;i<<=1){
memcpy(pa,x,sizeof(ll)*(i<<1));
memcpy(pb,aas,sizeof(ll)*i);
Mul(pa,pa,i<<1,pb,i<<1);Lim(pa,len,i<<1);
Mul(pa,pa,i<<1,pb,i<<1);Lim(pa,len,i<<1);
mul(pb,2,i);Red(pb,pa,i<<1);
memcpy(aas,pb,sizeof(ll)*(i<<1));
}
}
void Sqrt(ll *aas,ll *x,int lx){
aas[0]=1;
static ll pa[N],pb[N],pc[N];
memset(pa,0,sizeof(pa));
memset(pb,0,sizeof(pb));
memset(pc,0,sizeof(pc));
for(int i=1;i<lx*2;i<<=1){
memcpy(pa,x,sizeof(ll)*(i<<1));
memcpy(pb,aas,sizeof(ll)*i);
Inv(pc,pb,i);
Mul(pc,pa,i<<1,pc,i);Lim(pc,i<<2,i<<1);
Add(pc,pb,i<<1);mul(pc,qpow(2,mod-2),i<<1);
memcpy(aas,pc,sizeof(ll)*(i<<1));
}
}
void Ln(ll *aas,ll *x,int lx){
static ll pa[N],pb[N];
memset(pa,0,sizeof(pa));memset(pb,0,sizeof(pb));
Der(pa,x,lx);Inv(pb,x,lx);
Mul(aas,pa,lx,pb,lx);
Int(aas,aas,len);
}
void Exp(ll *aas,ll *x,int lx){
aas[0]=1;
static ll pa[N],pb[N];
memset(pa,0,sizeof(pa));memset(pb,0,sizeof(pb));
for(int i=1;i<lx*2;i<<=1){
memcpy(pa,aas,sizeof(ll)*i);
memcpy(pb,x,sizeof(ll)*(i<<1));
Ln(pa,aas,i);
Red(pb,pa,i<<1);pb[0]++;
Mul(aas,aas,i,pb,i<<1);Lim(aas,len,i<<1);
}
}
}
关于常数
说真的,最开始的时候被多项式的常数震撼到了。写多项式开根的时候实现不精细直接T飞了,以上模板是卡常后的结果,跑的还是不快。
卡常可以看九日的博客 。实现好的多项式板子效率还是可以接受的。