多项式
多项式
形如
这里的多项式实际上指的是形式幂级数, 也就是 (mathrm R[[x]]), 因为我们只关注它每一项的系数.
牛顿迭代法
给定定义域为多项式的函数 (G(x)), 求多项式 (F(x)), 使得
考虑倍增求出. 如果现在求出了低 ({lceil frac n2 ceil}) 项的系数, 即
, 那么利用泰勒展开可以得到
求逆
求 (F(x)), 满足 (F(x)A(x) equiv 1 pmod{x^n}).
设 (A(x)) 为原函数, 那么 (G(F(x)) = A(x) - frac{1}{F(x)}).
带入 ((1)),
多项式有逆的充要条件是常数项 (a_0 ot = 0), 递归终点 (f_0 = a_0^{-1}).
还可以设 (G(F(x)) = F(x)A(x) - 1), 推导略有不同:
带入 ((1)), 即
由于(A(x)F_0(x) - 1) 的低 (frac n2) 位均为 (0), 那么在模 (x^n) 意义下, 它乘上的数只需考虑低 (frac n2) 位即可.
而 (A(x)F_0(x)) 的低 (frac n2) 位中, 只有常数项为 (1), 其他项为 (0). 那么
上面的式子就可以化为
倍增即可求解.
求对数
求 (F(x)), 满足 (F(x) equiv ln A(x) pmod{x^n}).
多项式可以求对数的充要条件是常数项 (a_0 ot = 0), 经过积分得到常数项为 (0).
容易发现在求 (ln) 过程中 (A(x)) 的常数项被消成了 (1) (分子分母约分).
这里的求对数其实是假设 (a_0 = 1), 因为在模意义下, (ln x) 在 (x ot = 1) 处无定义.
求指数
求 (F(x)), 满足 (F(x) equiv e^{A(x)} pmod{x^n}).
设
那么, 通过牛顿迭代
多项式可以求指数的充要条件是常数项 (a_0 = 0), 递归终点 (f_0 = 1).
快速幂
求 (F(x)), 满足 (F(x) equiv A^k(x) pmod{x^n}). (k) 可为 (模意义下有意义的) 任意有理数.
由于 (ln) 过程的限制, 这个过程其实假设了 (a_0 = 1).
如果 (a_0) 为 (> 1) 正整数, 那么设多项式 (exp) 的递归终点为 (f_0 = a_0^k) (模意义下; 如果 (k) 为分数或负数需要求高次剩余/逆元) 即可;
如果 (a_0 = 0), 需要化为 (A'(x) = frac{A(x)}{x^l}), 其中 (A'(x)) 常数项不为 (0). 答案则为 ((A'(x))^k cdot x^{lk}).
多点求值
过于难写... 推荐秦九韶算法, 时间复杂度 (O(n^2)).
快速插值
给定 ((x_i, y_i) , forall i in {1,2,dotsc n}), 求一个 (n-1) 次多项式 (F(x)), 满足 (F(x_i) equiv y_i pmod p).
同上... 不想写快速插值...
拉格朗日插值
求单点值的时间复杂度 (O(n^2)).
拉格朗日反演
设(F(x)), (G(x))是形式幂级数, 满足 (F(G(x))=x), 称 (F(x)) 为 (G(x)) 的复合逆.
根据形式幂级数的性质, 可以得到 (F(x)) 常数项为 (0), 并且 (G(F(x)) = x).
可以证明
如果有 (F(x) = frac{x}{H(x)}), 还有
这个过程实际上是求 (frac{F(x)}{x}) 的逆的 (n) 次幂. 多项式求逆, 快速幂即可.
另外, 给定另一个多项式 (A(x)), 还可以求多项式的复合:
实现的细节
- 长度为
n, 次数为n-1 - 数组长度为超过
2*n的2的幂 - 必要时要清空
- 不要忘记取模!
代码
多项式全家桶:
const int nsz=4e5+50;
const ll nmod=998244353,g=3,ginv=332748118;
//998244353 g=3
//1004535809 g=3
int n;
ll a[nsz],b[nsz],ans[nsz];
ll qp(ll a,ll b=nmod-2){
ll res=1;
for(;b;b>>=1,a=a*a%nmod)if(b&1)res=res*a%nmod;
return res;
}
namespace npoly{
//l means length of array
ll len,l2,rev[nsz];
il void cp(ll *a,ll *b,int l0,int l1){
memcpy(a,b,l0*sizeof(ll));
memset(a+l0,0,(l1-l0)*sizeof(ll));
}
il void fftinit(int l0){
l2=0,len=1;
while(len<l0)++l2,len<<=1;
rep(i,0,len-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l2-1));//
}
void dft(ll *a,int l,int fl){
rep(i,0,l-1)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int i=1;i<l;i<<=1){
ll wn=qp(fl==1?g:ginv,(nmod-1)/(i<<1));
for(int j=0,p=i<<1;j<l;j+=p){
ll w=1;
for(int k=0;k<i;++k,w=w*wn%nmod){
ll x=a[j+k],y=w*a[j+k+i]%nmod;
a[j+k]=(x+y)%nmod,a[j+k+i]=(x-y)%nmod;
}
}
}
if(fl==-1){
ll v=qp(l);
rep(i,0,l-1)a[i]=a[i]*v%nmod;
}
}
void mul(ll *a,int l1,ll *b,int l2,ll *c,int l3=-1){
if(l3==-1)l3=l1+l2-1;
static ll c1[nsz],c2[nsz];
fftinit(l1+l2-1);
cp(c1,a,l1,len),cp(c2,b,l2,len);
dft(c1,len,1),dft(c2,len,1);
rep(i,0,len-1)c1[i]=c1[i]*c2[i]%nmod;
dft(c1,len,-1);
cp(c,c1,l3,l3);
}
void inv(ll *a,int l,ll *b){
if(l==1){b[0]=qp(a[0]);return;}
int l0=(l+1)>>1;
inv(a,l0,b);
static ll c1[nsz],c2[nsz];
fftinit(l<<1);//需要两倍长度dft保证乘法正确
cp(c1,a,l,len),cp(c2,b,l0,len);
dft(c1,len,1),dft(c2,len,1);
rep(i,0,len-1)c1[i]=c2[i]*(2-c1[i]*c2[i]%nmod)%nmod;
dft(c1,len,-1);
cp(b,c1,l,l);
}
void deriv(ll *a,int l,ll *b){ //b could be eq to a
rep(i,0,l-2)b[i]=a[i+1]*(i+1)%nmod;
b[l-1]=0;
}
void integ(ll *a,int l,ll *b){
repdo(i,l-2,0)b[i+1]=a[i]*qp(i+1)%nmod;
b[0]=0;
}
//a[0] should be 1
void ln(ll *a,int l,ll *b){
static ll c1[nsz],c2[nsz];
inv(a,l,c1),deriv(a,l,c2);
mul(c1,l,c2,l,b,l);
integ(b,l,b);
}
//a[0] should be 0
ll expzero=1;
void exp(ll *a,int l,ll *b){
if(l==1){b[0]=expzero;return;}
int l0=(l+1)>>1;
exp(a,l0,b);
static ll c1[nsz],c2[nsz];
fftinit(l<<1);
rep(i,l0,l-1)b[i]=0;
ln(b,l,c1);
rep(i,0,l-1)c1[i]=(a[i]-c1[i]+(i==0))%nmod;
rep(i,l,len-1)c1[i]=0;
cp(c2,b,l0,len);
dft(c1,len,1),dft(c2,len,1);
rep(i,0,len-1)c1[i]=c1[i]*c2[i]%nmod;
dft(c1,len,-1);
cp(b,c1,l,l);
}
void pow(ll *a,int l,ll p,ll *b){//suppose a[0]!=0;
static ll c[nsz];
ln(a,l,c);
rep(i,0,l-1)c[i]=c[i]*p%nmod;
// expzero=qp(a[0],p);
exp(c,l,b);
}
void sqrt(ll *a,int l,ll *b){pow(a,l,qp(2),b);}
void dncfft(ll *a,int l0,ll *b){
static ll c1[nsz];
rep(i,0,l0-1)c1[i]=-a[i];
c1[0]=(c1[0]+1)%nmod;
inv(c1,l0,b);
}
//tests
void testfft(){
int n,m;
cin>>n>>m,++n,++m;
rep(i,0,n-1)cin>>a[i];
rep(i,0,m-1)cin>>b[i];
mul(a,n,b,m,ans);
rep(i,0,n+m-2)cout<<(ans[i]+nmod)%nmod<<' ';
cout<<'
';
}
//3 1 2 1
//1 -2 3
void testinv(){
int n;
cin>>n;
rep(i,0,n-1)cin>>a[i];
inv(a,n,ans);
rep(i,0,n-1)cout<<(ans[i]+nmod)%nmod<<' ';
cout<<'
';
}
void testdnc(){
int n;
cin>>n;
rep(i,1,n-1)cin>>a[i];
dncfft(a,n,ans);
rep(i,0,n-1)cout<<(ans[i]+nmod)%nmod<<' ';
cout<<'
';
}
//3 1 2 1
//0 2 -1
void testln(){
int n;
cin>>n;
rep(i,0,n-1)cin>>a[i];
ln(a,n,ans);
rep(i,0,n-1)cout<<(ans[i]+nmod)%nmod<<' ';
cout<<'
';
}
//3 0 1 2
//1 1 499122179(5/2)
void testexp(){
int n;
cin>>n;
rep(i,0,n-1)cin>>a[i];
exp(a,n,ans);
rep(i,0,n-1)cout<<(ans[i]+nmod)%nmod<<' ';
cout<<'
';
}
void testsq(){
int n;
cin>>n;
rep(i,0,n-1)cin>>a[i];
sqrt(a,n,ans);
rep(i,0,n-1)cout<<(ans[i]+nmod)%nmod<<' ';
cout<<'
';
}
void testpow(){
int n,k;
cin>>n>>k;
rep(i,0,n-1)cin>>a[i];
pow(a,n,k,ans);
rep(i,0,n-1)cout<<(ans[i]+nmod)%nmod<<' ';
cout<<'
';
}
}
int main(){
ios::sync_with_stdio(0),cin.tie(0);
npoly::testexp();
return 0;
}
拉格朗日插值:
pair<ll,ll> li[nsz];
ll lag(ll x){
sort(li+1,li+n+1);
n=unique(li+1,li+n+1)-li-1;
ll res=0;
rep(i,1,n){
ll tmp=1;
rep(j,1,n)if(i!=j)tmp=tmp*(li[i].first-li[j].first)%nmod;
tmp=inv(tmp)*li[i].second%nmod;
rep(j,1,n)if(i!=j)tmp=tmp*(x-li[j].first)%nmod;
res=(res+tmp)%nmod;
}
return res;
}