码了个100行结构体多项式代码,把除了复合和复合逆以外的多项式操作都写了。顺便把多项式各种操作的原理写一下。
多项式乘法
这个。。NTT
class Array{
private:
vector<int>a;
public:
Array(const int size,const int f):a(size,f){}
void push(int n){a.push_back(n);}
Array(int* l=NULL,int* r=NULL){while(l!=r)push(*l),++l;}
inline int size(){return a.size();}
inline int& operator [] (const int x){return a[x];}
void resize(int n){a.resize(n);}
void clear(){a.clear();}
void swap(){reverse(a.begin(),a.end());}
};
int* NTT(const int len,Array& a,const bool Ty,int* r=NULL){
if(!r){
r=new int[len];
r[0]=0;int L=log2(len);
for(int i=0;i<len;i++)
r[i]=(r[i>>1]>>1)|((i&1)<<L-1);
}
for(int i=0;i<len;i++)
if(i<r[i])swap(a[i],a[r[i]]);
for(int i=1;i<len;i<<=1){
int T=poww(Ty?g:332748118,(mod-1)/(i<<1),mod);
for(int W=i<<1,j=0;j<len;j+=W){
ll omega=1;
for(int k=0;k<i;++k,omega=omega*T%mod){
ll x(a[j+k]),y(omega*a[i+j+k]%mod);
a[j+k]=x+y;(a[j+k]>mod)&&(a[j+k]-=mod);
a[i+j+k]=x-y+mod;(a[i+j+k]>mod)&&(a[i+j+k]-=mod);
}
}
}
return r;
}
Array Mul(Array x,Array y){
int n=x.size()-1,m=y.size()-1;
int limit=1;
while(limit<=n+m)limit<<=1;
Array ans;
x.resize(limit+1);
y.resize(limit+1);
int *r;
r=NTT(limit,x,1);
NTT(limit,y,1,r);
for(int i=0;i<=limit;i++) x[i]=1ll*x[i]*y[i]%mod;
NTT(limit,x,0,r);
int tem=poww(limit,mod-2,mod);
for(int i=0;i<=n+m;i++) x[i]=(ll)x[i]*tem%mod;
x.resize(n+m+1);
return x;
}
牛顿迭代
就是一个在实数域和复数域上近似求解方程的方法。
公式为:
已
知
方
程
f
(
x
)
=
0
,
x
n
+
1
=
x
n
−
f
(
x
n
)
f
′
(
x
n
)
已知方程f(x)=0, x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}
已知方程f(x)=0,xn+1=xn−f′(xn)f(xn)
这个
x
n
x_n
xn不断迭代就会以平方收敛,靠近方程的解,也就是说牛顿迭代每次迭代都会使有效数字翻一倍
那么,把这个方法拓展到多项式上,就可以做到乘法逆元、exp的运算,迭代logN次就能得到正确答案。
乘法逆元
设关于 f ( x ) f(x) f(x)的方程 f ( x ) g ( x ) = 1 f(x)g(x)=1 f(x)g(x)=1,即 f ( x ) g ( x ) − 1 = 0 f(x)g(x)-1=0 f(x)g(x)−1=0
把这个式子套牛顿迭代,把 f n ( x ) f_n(x) fn(x)替代 x n x_n xn:
f n + 1 ( x ) = f n ( x ) − f n ( x ) g ( x ) − 1 g ( x ) f_{n+1}(x)=f_n(x)-\frac{f_n(x)g(x)-1}{g(x)} fn+1(x)=fn(x)−g(x)fn(x)g(x)−1
化简得:
f n + 1 ( x ) = 2 f n ( x ) − g ( x ) f n 2 ( x ) f_{n+1}(x)=2f_n(x)-g(x)f_n^2(x) fn+1(x)=2fn(x)−g(x)fn2(x)
然后递归就好了,递归logN次,同时注意递归边界,即 f ( x ) f(x) f(x) 为0次多项式时的解,就是 g ( x ) g(x) g(x) 0次项系数的模意义下的乘法逆元。
设这个复杂度为 T ( N ) T(N) T(N),那么 T ( N ) = T ( N / 2 ) + O ( N l o g N ) T(N)=T(N/2)+O(NlogN) T(N)=T(N/2)+O(NlogN),可以证明(我不会证), T ( N ) T(N) T(N)其实就是 O ( N l o g N ) O(NlogN) O(NlogN)的,但是常数更大。
void Rev(Array &x,Array y){
int n=x.size()-1,m=y.size()-1;
int limit=1;
while(limit<=n+m)limit<<=1;
Array ans;
x.resize(limit+1);
y.resize(limit+1);
int *r;
r=NTT(limit,x,1);
NTT(limit,y,1,r);
for(int i=0;i<=limit;i++) x[i]=(ll)(2ll-(ll)x[i]*y[i]%mod+mod)%mod*y[i]%mod;
NTT(limit,x,0,r);
int tem=poww(limit,mod-2,mod);
for(int i=0;i<=n+m;i++) x[i]=(ll)x[i]*tem%mod;
x.resize(n+m+1);
}
Array Inv(Array a){
int N=a.size();
if(N==1)return Array(1,poww(a[0],mod-2,mod));
Array b=a;b.resize(N+1>>1);
b=Inv(b);b.resize(N);
Rev(a,b);
a.resize(N);
return a;
}
ln
其实这个。。就是一个式子:
l
n
(
f
(
x
)
)
=
∫
f
′
(
x
)
f
(
x
)
d
x
ln(f(x))=\int\frac{f'(x)}{f(x)}dx
ln(f(x))=∫f(x)f′(x)dx
但是它需要逆元运算,所以复杂度是O(NlogN)
Array Drv(Array a){ //求导
for(int i=0;i<a.size()-1;i++) a[i]=1ll*a[i+1]*(i+1)%mod;
a[a.size()-1]=0;
return a;
}
Array Idf(Array a){ //不定积分
for(int i=a.size()-1;i>0;i--) a[i]=1ll*a[i-1]*poww(i,mod-2,mod)%mod;
a[0]=0;
return a;
}
Array ln(Array a){
return Idf(Mul(Drv(a),Inv(a)));
}
exp
同样是牛顿迭代法,推式子即可:
f
n
+
1
(
x
)
=
f
n
(
x
)
−
ln
(
f
n
(
x
)
)
−
g
(
x
)
1
f
n
(
x
)
f_{n+1}(x)=f_n(x)-\frac{\ln(f_n(x))-g(x)}{\frac{1}{f_n(x)}}
fn+1(x)=fn(x)−fn(x)1ln(fn(x))−g(x)
f
n
+
1
(
x
)
=
f
n
(
x
)
(
1
−
ln
(
f
n
(
x
)
+
g
(
x
)
)
f_{n+1}(x)=f_n(x)(1-\ln(f_n(x)+g(x))
fn+1(x)=fn(x)(1−ln(fn(x)+g(x))
然后做法就和之前一样了,这个复杂度也是O(NlogN),但是是套上逆元的O(NlogN),也就是常数是逆元的常数的平方,常数很大。
Array Exp(Array a){
int N=a.size();
if(N==1) return Array(1,1);
Array b=a;b.resize(N+1>>1);
b=Exp(b);b.resize(N);
Array b1=ln(b);
for(int i=0;i<a.size();i++) b1[i]=1ll*(a[i]-b1[i]+mod)%mod;
b1[0]=1ll*(b1[0]+1)%mod;
a=Mul(b,b1);
return a;
}
平方根
我一开始想法是求ln,然后除2再exp回去,但是T了(略略略),所以还是用牛顿迭代。
需要注意的是,这里的边界不太好处理,要处理模意义的平方根,所以用了bsgs算法。
int bsgs(int y,int z,int p){
int m=sqrt(p)+1;
map<int,int> mp;
for(int i=0,t=z;i<m;i++,t=1ll*y*t%p) mp[t]=i;
for(int i=1,tt=poww(y,m,p),t=tt;i<=m+1;i++,t=1ll*t*tt%p){
int k=mp[t];if(k==0) continue;
return i*m-k;
}
}
int sqrtt(int x){
return min(poww(3,bsgs(3,x,mod)/2,mod),mod-poww(3,bsgs(3,x,mod)/2,mod));
}
Array Sqrt(Array a){
int N=a.size();
if(N==1) return Array(1,sqrtt(a[0]));
Array b=a;b.resize(N+1>>1);
b=Sqrt(b);b.resize(N);
a=Mul(a,Inv(b));
a.resize(N);
int inv2=poww(2,mod-2,mod);
for(int i=0;i<a.size();i++) a[i]=1ll*(a[i]+b[i])*inv2%mod;
return a;
}