标签:efi ++ return 原函数 tor 写代码 size lan c++
前置知识\(\text{NTT}\)
所有操作均在对\(\text{998244353}\)取模下进行
所有的操作均没有经过任何卡常
所有的操作均用$\text{vector} $来实现,主要是为了理清思路,并且清零问题上会比较容易解决,同时对于每次计算完多项式的长度的要求会显得更加严格
相比经过卡常的代码来说,可读性应该强很多吧
代码总览
ll qpow(ll x,ll k) { // 朴素的快速乘
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
namespace Polynomial{
typedef vector <int> Poly;
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
void Show(Poly a){ for(int i:a) printf("%d ",i); puts(""); }
int rev[N];
int PreMake(int n){
int R=1,cc=-1;
while(R<n) R<<=1,cc++;
rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
return R;
}
void NTT(int n,Poly &a,int f){ //朴素的NTT模板
rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int w=qpow(f==1?3:(P+1)/3,(P-1)/i/2);
for(int l=0;l<n;l+=i*2) {
ll e=1;
for(int j=l;j<l+i;++j,e=e*w%P) {
int t=a[j+i]*e%P;
a[j+i]=a[j]-t,Mod2(a[j+i]);
a[j]+=t,Mod1(a[j]);
}
}
}
if(f==-1) {
ll base=qpow(n,P-2);
rep(i,0,n-1) a[i]=a[i]*base%P;
}
}
Poly operator * (Poly a,Poly b){
int n=a.size()+b.size()-1,R=PreMake(n);
a.resize(R),b.resize(R);
NTT(R,a,1),NTT(R,b,1);
rep(i,0,R-1) a[i]=1ll*a[i]*b[i]%P;
NTT(R,a,-1);
a.resize(n);
return a;
}
Poly operator + (Poly a,Poly b) { rep(i,0,a.size()-1) a[i]+=b[i],Mod1(a[i]); return a; }
Poly operator - (Poly a,Poly b) { rep(i,0,a.size()-1) a[i]-=b[i],Mod2(a[i]); return a; }
Poly Inv(Poly a) { // 多项式求逆
int n=a.size();
if(n==1) { Poly tmp; tmp.pb(qpow(a[0],P-2)); return tmp; }
Poly b=a; b.resize((n+1)/2); b=Inv(b);
int R=PreMake(n<<1);
a.resize(R),b.resize(R);
NTT(R,a,1),NTT(R,b,1);
rep(i,0,R-1) a[i]=(2-1ll*a[i]*b[i]%P+P)*b[i]%P;
NTT(R,a,-1);
a.resize(n);
return a;
}
Poly operator / (vector <int> a,vector <int> b){
reverse(a.begin(),a.end()),reverse(b.begin(),b.end());
int n=a.size(),m=b.size();
b.resize(n-m+1),b=Inv(b);
a=a*b,a.resize(n-m+1);
reverse(a.begin(),a.end());
return a;
}
Poly operator % (vector <int> a,vector <int> b) {
a=a-a/b*b;
a.resize(b.size()-1);
return a;
}
Poly Sqrt(vector <int> a){ // 开根号
int n=a.size();
if(n==1) { Poly tmp; tmp.pb(1); return tmp; }
Poly b=a; b.resize((n+1)/2),b=Sqrt(b),b.resize(n);
Poly c=Inv(b);
int R=PreMake(n*2);
a.resize(R),c.resize(R);
NTT(R,a,1),NTT(R,c,1);
rep(i,0,R-1) a[i]=1ll*a[i]*c[i]%P;
NTT(R,a,-1);
a.resize(n);
rep(i,0,n-1) a[i]=1ll*(P+1)/2*(a[i]+b[i])%P;
return a;
}
Poly Deri(Poly a){ // 导数
rep(i,1,a.size()-1) a[i-1]=1ll*i*a[i]%P;
a.pop_back();
return a;
}
int Mod_Inv[N];
Poly IDeri(Poly a) { // 原函数
Mod_Inv[0]=Mod_Inv[1]=1;
rep(i,2,a.size()+1) Mod_Inv[i]=1ll*(P-P/i)*Mod_Inv[P%i]%P;
a.pb(0);
drep(i,a.size()-1,1) a[i]=1ll*a[i-1]*Mod_Inv[i]%P;
a[0]=0;
return a;
}
Poly ln(Poly a){
int n=a.size();
a=Inv(a)*Deri(a);
a.resize(n-1);
return IDeri(a);
}
Poly Exp(Poly a){
int n=a.size();
if(n==1) { Poly tmp; tmp.pb(1); return tmp; }
Poly b=a; b.resize((n+1)/2),b=Exp(b),b.resize(n);
Poly c=ln(b);
int R=PreMake(n<<1);
a.resize(R),b.resize(R),c.resize(R);
NTT(R,b,1),NTT(R,a,1),NTT(R,c,1);
rep(i,0,R-1) a[i]=1ll*b[i]*(1-c[i]+a[i]+P)%P;
NTT(R,a,-1),a.resize(n);
return a;
}
Poly Pow(Poly a,int k){
a=ln(a);
rep(i,0,a.size()-1) a[i]=1ll*a[i]*k%P;
return exp(a);
}
#undef Mod1
#undef Mod2
}
求 \(G(x)\equiv \frac{1}{F(x)} (\mod x^n)\)
考虑递归求解,设已经求出了
其中递归边界是\(n=1\)时,\(G(x)=\frac{1}{F(0)} (\mod P)\)
则
注意通过平方可以扩大模数,这很常用
两边乘上\(F(x)\)
求\(G(x)^2\equiv F(x) (\mod x^n)\)
同样的,递归求解,设已经求出了,递归边界是\(n=1\)时,\(G(x)=1\)
\(H(x)\equiv G(x) (\mod x^{\lceil \frac{n}{2}\rceil})\)
\(G(x)\equiv \ln F(x) (\mod x^n)\)
\(G‘(x)\equiv F‘(x)\frac{1}{F(x)} (\mod x^n)\)
求出\(G‘(x)\),然后求原函数即可
这个推导其实非常有趣的
不说了,先上一下泰勒公式
其中\(f^{(i)}(x)\)表示\(f\)的\(i\)阶导数
即通过\(f(x)\)的导数在\(x_0\)处的函数值累加得到近似的多项式函数,累加到无穷之后认为相等
题目的问题是求
\(G(x)=e^{ F(x)},F(x)=\ln G(x)\)
\(\ln G(x)-F(x)\equiv 0 (\mod x^n)\)
那么我们把元\(x\)推广为多项式,把题目转化为
对于函数\(f(G)=\ln G-F\)
求出在\(\mod x^n\)意义下的零点
其中\(f(x)=\ln x-c\)
考虑迭代求解,设已经求出\(H(x)=e^{F(x)}( \mod x^{\frac{n}{2}})\)
带入\(f(G)\)在\(H\)上的泰勒展开式
和上面一样可以得到的是
\(H(x)\equiv G(x) (\mod x^{\frac{n}{2}})\)
\((H(x)-G(x))^2 \equiv 0 (\mod x^n)\)
所以只用考虑\(i=0,1\)的两项,后面的\((G-H)^i\)均为\(0\)
那么\(f(G)=f(H)+f‘(H) \cdot (G-H)\)
求导得到\(f‘(x)=\frac{1}{x}\)
所以\(f(G)=\ln H-F+\frac{G-H}{H}\equiv 0 (\mod x^n)\)
\(G=H(F-\ln H+1)\)
所以就可以好好得递归写代码啦
\(G(x)\equiv F(x)^k(\mod x^n)\)
\(\ln G(x)=k \ln F(x) (\mod x^n)\)
求出\(\ln G(x)\)之后,\(exp\)回来即可
很显然这个方法对于开根号也是适用的
标签:efi ++ return 原函数 tor 写代码 size lan c++
原文地址:https://www.cnblogs.com/chasedeath/p/12801809.html