前言
这里的全家桶目前只包括了(ln,exp,sqrt)。还有一些类似于带余数模，快速幂之类用的比较少的有时间再更，(NTT)这种前置知识这里不多说。
还有一些基本的导数和微积分内容要了解，建议不懂的可以先去翻翻高二数学书。
之后多项式算法基本是一环扣一环的，所以前面的看不懂对于后面的理解会造成很大影响。
本博客涉及内容偏浅
Tips
这里是一些我个人的模板书写习惯
- 习惯相关的问题：默认将读入的(n)变为(2)的整数次幂形式，目前为止这样的做法都不会影响正确性
- 正确性相关的问题：模板书写应满足使用的中间变量不重复，如在求(ln)和逆元时不应该使用重复的中间变量（可能会导致信息丢失）
- 正确性相关的问题：在每次进行任意操作前应保证在操作值域内不会有上次的信息参与（要清零）
参考资料
正题
多项式求逆
题目大意
给出一个多项式(F)，求出一个(G)使得
[F(x)*G(x)equiv1(mod x^n)
]
分析
利用经典的倍增思想，假设我们已知多项式(G'(x))满足
[F(x)G'(x)equiv1(mod x^{frac{n}{2}})
]
又有
[F(x)G(x)equiv 1(mod x^{frac{n}{2}})
]
就有了
[F(x)(G(x)-G'(x))equiv 0(mod x^{frac{n}{2}})Rightarrow G(x)-G'(x)equiv 0(mod x^{frac{n}{2}})
]
然后两边同时平方，后面的模数同理也要平方
[G(x)^2-G(x)G'(x)+G'(x)^2equiv 0(mod x^n)
]
再乘上一个(F(x))
[F(x)G(x)^2-F(x)G(x)G'(x)+F(x)G'(x)^2equiv 0(mod x^n)
]
又因为(F(x)G(x)equiv 1(mod x^n))，所以就有
[G(x)-G'(x)+F(x)G'(x)^2equiv0(mod x^n)Rightarrow G(x)equiv G'(x)-F(x)G'(x)^2(mod x^n)
]
然后倍增就好了，时间复杂度是类似于(T(n)=T(frac{n}{2})+nlog n)的形式所以是(O(nlog n))的。
code
比较远古的代码所以码风有点不同
// luogu-judger-enable-o2
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
using namespace std;
const ll N=1e6+100,XJQ=998244353,G=3,Gi=332748118;
const double Pi=acos(-1);
ll n,m,l,a[N<<2],b[N<<2],c[N<<2],r[N<<2];
ll power(ll x,ll b)
{
ll ans=1;
while(b)
{
if(b&1) ans=ans*x%XJQ;
x=x*x%XJQ;b>>=1;
}
return ans;
}
void ntt(ll *f,ll n,ll op)
{
for(ll i=0;i<n;i++)
if(i<r[i])swap(f[i],f[r[i]]);
for(ll p=2;p<=n;p<<=1){
ll len=p>>1,tmp=power(G,(XJQ-1)/p);
for(ll k=0;k<n;k+=p){
ll buf=1;
for(ll i=k;i<k+len;i++){
ll tt=buf*f[len+i]%XJQ;
f[len+i]=(f[i]-tt+XJQ)%XJQ;
f[i]=(f[i]+tt)%XJQ;
buf=buf*tmp%XJQ;
}
}
}
if(op==1) return;
int inv=power(n,XJQ-2);reverse(f+1,f+n);
for(int i=0;i<n;i++) f[i]=f[i]*inv%XJQ;
}
void work(ll *a,ll *b,ll l)
{
if(l==1){b[0]=power(a[0],XJQ-2);return;}
work(a,b,(l+1)>>1);
ll cnt;
for(cnt=1;cnt<(l<<1);cnt<<=1);
for(ll i=1;i<cnt;i++)
r[i]=(r[i>>1]>>1)|((i&1)?cnt>>1:0);
for(ll i=0;i<l;i++) c[i]=a[i];
for(ll i=l;i<cnt;i++) c[i]=0;
ntt(c,cnt,1);ntt(b,cnt,1);
for(ll i=0;i<cnt;i++)
b[i]=(2-b[i]*c[i]%XJQ+XJQ)%XJQ*b[i]%XJQ;
ntt(b,cnt,-1);
for(ll i=l;i<cnt;i++) b[i]=0;
}
int main()
{
scanf("%lld",&n);
for(ll i=0;i<n;i++)
scanf("%lld",&a[i]);
work(a,b,n);
for(ll i=0;i<n;i++)
printf("%lld ",b[i]);
}
多项式导数相关
后面就要开始用到高二的知识了
多项式求导
后面定义(f')表示多项式(f)的求导，定义(a_i=f(x)[i])，那么有
[f'(x)=sum_{i=0}^na_{i+1}(i+1)x^i
]
大体就是把所有位乘上(i)再往前移动，对导数有了解的话应该能理解，这里就不给出推导了
多项式积分
同理定义(a_i=f(x)[i])，那(f)的积分(g)就有
[g(x)=sum_{i=0}^nfrac{a_{i-1}}{i}x^i
]
和上面同理，需要知道积分和求导是逆运算。
多项式复合
定义复合函数
[F(G(x))=sum_{i=0}F[i]G(x)^i
]
看上去没什么用，之后再说
泰勒公式
这里开始就暂时不和多项式有多挂钩了。
[f(x)=lim_{n o infty}sum_{i=0}^{n}frac{f^{(n)}(x_0)}{i!}(x-x_0)^i+R_n(x)
]
其中(f^{(n)})表示(f)的(n)阶求导，(R_n(x))是余项。
看上去没有什么用，是牛迭的基础，可以直接记牛迭的结论。
到时候会有泰勒公式的常用写法
牛顿迭代
这里就不用泰勒公式的推导了，直接感性点理解快速过一下。
我们知道导数
[f'(x)=lim_{Delta x o0}frac{f(x+Delta x)-f(x)}{Delta x}
]
是可以理解为求某个函数在((x,f(x)))处的切点，而牛顿迭代正是利用这个原理求函数的近似零点，可以先看一张生动的图（来自维基百科）。
就是先找一个点((x,f(x)))，然后求它在函数图像上的切线，之后这条切线与(x)有交的位置(x')再带入(x)之后继续这个过程。
这个过程中求得的(x)在不断逼近原点，这样就可以求出一个函数(0)点的近似解。
当然牛顿迭代显然并不是对所有函数都适用的，但是对于我们需要解决的多项式问题来说足够了。
然后要上泰勒公式了
[f(x)=lim_{n o infty}sum_{i=0}^{n}frac{f^{(n)}(x_0)}{i!}(x-x_0)^i+R_n(x)
]
这里只拿(i=1)的情况来展开一下，再定义一个(phi(x)approx f(x))就是
[f(x)approx phi(x)=f'(x_0)(x-x_0)-f(x_0)
]
然后如果求(f(x)=0)就是近似的求(phi(x)=0)也就是
[f'(x_0)(x-x_0)-f(x_0)=0
]
就有
[x=x_0-frac{f(x_0)}{f'(x_0)}Rightarrow x_{n+1}=x_{n}-frac{f(x_n)}{f'(x_n)}
]
这个递推式子。
然后就可以快速近似求解了。
多项式复合零点
那么现在就是牛顿迭代的实战时间了，题目是给出一个多项式(G(x))，要求求出一个(f(x))使得(G(f(x))equiv 0(mod x^n))。
拿多项式来泰勒展开推导或者直接用上面的牛迭结论。设(f_t)满足(G(f_t)equiv 0(mod x^{2^t}))，那么就有结论
[f_tequiv f_{t-1}-frac{G(f_{t-1})}{G'(f_{t-1})}(mod {x^{2^t}})
]
这里还是推导一下吧，先把(f)给泰勒展开了（这里换成了一个比较常用的写法）
[G(f_t)equiv G(f_{t-1})+sum_{i=1}^inftyfrac{G^{(i)}(f_{t-1})}{i!}*(f_t-f_{t-1})^i(mod x^{2^t})
]
又因为(G(f_{t-1})equiv 0(mod x^{2^{t-1}}))和(G(f_t)equiv 0(mod x^{2^t}))。所以(f_t-f_{t-1})的前(2^{t-1})项都是(0)，那么((f_t-f_{t-1})^2)的前(2^t)都是0，也就当(igeq 2)的时候后面的项都被模掉了，所以式子就变得很简单了。
[G(f_t)equiv G(f_{t-1})+G'(f_{t-1})*(f_t-f_{t-1})(mod x^{2^t})
]
把(G(f_t)=0)带入就有
[f_t=f_{t-1}-frac{G(f_{t-1})}{G'(f_{t-1})}(mod x^{2^t})
]
式子到这里就得根据具体情况化简了，然后练练手？
多项式开根
题目大意
给出一个多项式(F)，求一个多项式(G)满足
[G(x)^2=F(x)(mod x^{n})
]
分析
下面的(G)和上面的要求的(G)不同
如果我们求出一个多项式(G(f)=f^2-F(x))。如果(G(f)equiv 0(mod x^n))的解就是(f(x)equiv sqrt{F(x)}(mod x^n))的解了。
之后直接代前面多项式零点求值的东西就有
[f_t=f_{t-1}-frac{G(f_{t-1})}{G'(f_{t-1})}Rightarrow f_t=f_{t-1}-frac{f_{t-1}^2-F(x)}{2f_{t-1}}
]
嗯，那个(2f_{t-1})是对(G)手动求导的结果
时间复杂度也是类(T(n)=T(frac{n}{2})+nlog n)的形式所以还是(O(nlog n))
code
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const ll N=6e5+10,P=998244353,inv2=(P+1)/2;
ll n,a[N],b[N],r[N];
ll t1[N],t2[N],t3[N],t4[N];
ll power(ll x,ll b){
ll ans=1;
while(b){
if(b&1)ans=ans*x%P;
x=x*x%P;b>>=1;
}
return ans;
}
ll GetL(ll len){
ll n=1;
while(n<=len)n<<=1;
for(ll i=0;i<n;i++)
r[i]=(r[i>>1]>>1)|((i&1)?(n>>1):0);
return n;
}
void NTT(ll *f,ll n,ll op){
for(ll i=0;i<n;i++)
if(i<r[i])swap(f[i],f[r[i]]);
for(ll p=2;p<=n;p<<=1){
ll len=p>>1,tmp=power(3,(P-1)/p);
if(op==-1)tmp=power(tmp,P-2);
for(ll k=0;k<n;k+=p){
ll buf=1;
for(ll i=k;i<k+len;i++){
ll tt=f[i+len]*buf%P;
f[i+len]=(f[i]-tt+P)%P;
f[i]=(f[i]+tt)%P;
buf=buf*tmp%P;
}
}
}
if(op==-1){
ll invn=power(n,P-2);
for(ll i=0;i<n;i++)
f[i]=f[i]*invn%P;
}
return;
}
void GetInv(ll *f,ll *g,ll n){
if(n==1){g[0]=power(f[0],P-2);return;}
GetInv(f,g,n>>1);ll l=GetL(n);
for(ll i=0;i<n;i++)t1[i]=f[i],t2[i]=g[i];
for(ll i=n;i<l;i++)t1[i]=t2[i]=0;
NTT(t1,l,1);NTT(t2,l,1);
for(ll i=0;i<l;i++)t1[i]=t1[i]*t2[i]%P*t2[i]%P;
NTT(t1,l,-1);
for(ll i=0;i<n;i++)g[i]=(2*g[i]-t1[i]+P)%P;
return;
}
void Sqrt(ll *f,ll *g,ll n){
if(n==1){g[0]=1;return;}
Sqrt(f,g,n>>1);ll l=GetL(n<<1);
for(ll i=0;i<n;i++)t3[i]=f[i],t4[i]=0;
for(ll i=n;i<l;i++)t3[i]=t4[i]=0;
GetInv(g,t4,n);l=GetL(n<<1);
NTT(t3,l,1);NTT(t4,l,1);
for(ll i=0;i<l;i++)t3[i]=t3[i]*t4[i]%P;
NTT(t3,l,-1);
for(ll i=0;i<n;i++)g[i]=(g[i]+t3[i])*inv2%P;
return;
}
signed main()
{
scanf("%lld",&n);
for(ll i=0;i<n;i++)
scanf("%lld",&a[i]);
ll m=GetL(n);
Sqrt(a,b,m);
for(ll i=0;i<n;i++)
printf("%lld ",b[i]);
return 0;
}
多项式ln
题目大意
给出一个多项式(F)，求一个多项式(G)满足
[G(x)equiv ln(F(x))(mod x^n)
]
分析
这个题不用牛顿迭代，对(G)求个导，因为是复合函数直接展开
[G'(x)equiv ln'(F(x))*F(x)(mod x^n)Rightarrow G'(x)equivfrac{F(x)}{F'(x)}(mod x^n)
]
这个推导要用到的有(ln'(x)=frac{1}{x})。
就是算(frac{F(x)}{F'(x)})再积分就好了
时间复杂度(O(nlog n))
code
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const ll N=4e5+10,P=998244353;
ll n,r[N],f[N],g[N];
ll t1[N],t2[N],t3[N],t4[N];
ll power(ll x,ll b){
ll ans=1;
while(b){
if(b&1)ans=ans*x%P;
x=x*x%P;b>>=1;
}
return ans;
}
ll GetL(ll len){
ll n=1;
while(n<=len)n<<=1;
for(ll i=0;i<n;i++)
r[i]=(r[i>>1]>>1)|((i&1)?(n>>1):0);
return n;
}
void NTT(ll *f,ll n,ll op){
for(ll i=0;i<n;i++)
if(i<r[i])swap(f[i],f[r[i]]);
for(ll p=2;p<=n;p<<=1){
ll len=p>>1,tmp=power(3,(P-1)/p);
if(op==-1)tmp=power(tmp,P-2);
for(ll k=0;k<n;k+=p){
ll buf=1;
for(ll i=k;i<k+len;i++){
ll tt=f[i+len]*buf%P;
f[i+len]=(f[i]-tt+P)%P;
f[i]=(f[i]+tt)%P;
buf=buf*tmp%P;
}
}
}
if(op==-1){
ll invn=power(n,P-2);
for(ll i=0;i<n;i++)
f[i]=f[i]*invn%P;
}
return;
}
void GetInv(ll *f,ll *g,ll n){
if(n==1){g[0]=power(f[0],P-2);return;}
GetInv(f,g,n>>1);ll m=GetL(n);
for(ll i=0;i<n;i++)t1[i]=f[i],t2[i]=g[i];
for(ll i=n;i<m;i++)t1[i]=t2[i]=0;
NTT(t1,m,1);NTT(t2,m,1);
for(ll i=0;i<m;i++)t1[i]=t1[i]*t2[i]%P*t2[i]%P;
NTT(t1,m,-1);
for(ll i=0;i<n;i++)g[i]=(2*g[i]-t1[i]+P)%P;
return;
}
void GetD(ll *f,ll *g,ll n){
for(ll i=0;i<n;i++)
g[i]=f[i+1]*(i+1)%P;
g[n-1]=0;return;
}
void GetJ(ll *f,ll *g,ll n){
for(ll i=1;i<n;i++)
g[i]=f[i-1]*power(i,P-2)%P;
g[0]=0;return;
}
void GetLn(ll *f,ll *g,ll n){
n=GetL(n);
GetD(f,t3,n);GetInv(f,t4,n);
n=GetL(n);NTT(t3,n,1);NTT(t4,n,1);
for(ll i=0;i<n;i++)t3[i]=t3[i]*t4[i]%P;
NTT(t3,n,-1);GetJ(t3,g,n);
return;
}
signed main()
{
scanf("%lld",&n);
for(ll i=0;i<n;i++)
scanf("%lld",&f[i]);
GetLn(f,g,n);
for(ll i=0;i<n;i++)
printf("%lld ",g[i]);
return 0;
}
多项式exp
题目大意
给出多项式(F)，求一个多项式(G)满足
[G(x)equiv e^{F(x)}(mod x^n)
]
解题思路
这个应该是最麻烦的了，和开根一样的思路
定义一个复合函数(G(f)=ln(f)-F(x))那么当(G(f)=0)的解就是答案了。
然后同理直接上倍增加泰勒展开
[f_tequiv f_{t-1}-frac{G(f_{t-1})}{G'(f_{t-1})}(mod x^{2^t})Rightarrow f_tequiv f_{t-1}-(ln(f_{t-1})+F(x))*f_{t-1}(mod x^{2^t})
]
然后时间复杂度依旧是(T(n)=T(frac{n}{2})+nlog n)，所以还是(O(nlog n))的
然后上个(ln)和求逆就好了
code
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const ll N=8e5+10,P=998244353;
ll n,m,r[N],a[N],b[N];
ll t1[N],t2[N],t3[N],t4[N],t5[N],t6[N];
ll power(ll x,ll b){
ll ans=1;
while(b){
if(b&1)ans=ans*x%P;
x=x*x%P;b>>=1;
}
return ans;
}
void GetL(ll len){
n=1;while(n<=len)n<<=1;
for(ll i=0;i<n;i++)
r[i]=(r[i>>1]>>1)|((i&1)?(n>>1):0);
return;
}
void NTT(ll *f,ll op){
for(ll i=0;i<n;i++)
if(i<r[i])swap(f[i],f[r[i]]);
for(ll p=2;p<=n;p<<=1){
ll len=p>>1,tmp=power(3,(P-1)/p);
if(op==-1)tmp=power(tmp,P-2);
for(ll k=0;k<n;k+=p){
ll buf=1;
for(ll i=k;i<k+len;i++){
ll tt=f[i+len]*buf%P;
f[i+len]=(f[i]-tt+P)%P;
f[i]=(f[i]+tt)%P;
buf=buf*tmp%P;
}
}
}
if(op==-1){
ll invn=power(n,P-2);
for(ll i=0;i<n;i++)
f[i]=f[i]*invn%P;
}
return;
}
void GetInv(ll *f,ll *g,ll m){
if(m==1){g[0]=power(f[0],P-2);return;}
GetInv(f,g,m>>1);GetL(m);
for(ll i=0;i<m;i++)t1[i]=f[i],t2[i]=g[i];
for(ll i=m;i<n;i++)t1[i]=t2[i]=0;
NTT(t1,1);NTT(t2,1);
for(ll i=0;i<n;i++)t1[i]=t1[i]*t2[i]%P*t2[i]%P;
NTT(t1,-1);
for(ll i=0;i<m;i++)g[i]=(2*g[i]-t1[i]+P)%P;
return;
}
void GetD(ll *f,ll *g,ll n){
for(ll i=0;i<n-1;i++)
g[i]=f[i+1]*(i+1)%P;
g[n-1]=0;return;
}
void GetJ(ll *f,ll *g,ll n){
for(ll i=1;i<n;i++)
g[i]=f[i-1]*power(i,P-2)%P;
g[0]=0;return;
}
void GetLn(ll *f,ll *g,ll m){
GetL(m);GetD(f,t3,n);GetInv(f,t4,n);
GetL(m);GetL(n);NTT(t3,1);NTT(t4,1);
for(ll i=0;i<n;i++)t3[i]=t3[i]*t4[i]%P;
NTT(t3,-1);GetJ(t3,g,n);
for(int i=0;i<n;i++)t3[i]=t4[i]=0;
return;
}
void GetExp(ll *f,ll *g,ll m){
if(m==1){g[0]=1;return;}
GetExp(f,g,m>>1);GetLn(g,t5,m);GetL(m);
for(ll i=0;i<m;i++)t6[i]=f[i];
for(ll i=m;i<n;i++)t5[i]=t6[i]=0;
NTT(t5,1);NTT(t6,1);NTT(g,1);
for(ll i=0;i<n;i++)
g[i]=g[i]*(1-t5[i]+t6[i]+P)%P;
NTT(g,-1);for(ll i=m;i<n;i++)g[i]=0;
return;
}
signed main()
{
scanf("%lld",&m);
for(ll i=0;i<m;i++)
scanf("%lld",&a[i]);
GetL(m);GetExp(a,b,n);
for(ll i=0;i<m;i++)
printf("%lld ",b[i]);
return 0;
}