前言
因为之前的 Ⅰ 写的太屑了，所以又写了篇更屑的Ⅱ，其实是不同时间的个人垃圾理解了..
牛顿迭代
- 要求解(f(x)=0)[frac{f(x_0)}{x_0-x}=f'(x_0)\ x=x_0-frac{f(x_0)}{f'(x_0)} ]不断往下迭代可以找到近似解。
多项式牛顿迭代
和普通的牛顿迭代其实一样
[F=F_0-frac{G(F_0)}{G'(F_0)}
]
如何理解近似解？
从(mod x^i)开始向上(mod x^{i+1}...)
- 多项式求逆[A*Bequiv 1pmod {x^n} ]构造[G(B)=A*B-1 ]牛顿迭代找根
[ egin{aligned}
B&=B_0-frac{AB_0-1}{A}\
&=B_0-B_0(AB_0-1)\
&=2B_0-AB_0^2
end{aligned}
]
向上倍增即可。
-
多项式开根
[B^2equiv Apmod {x^n} ]构造
[G(B)=A-B^2 ]牛顿迭代找根
[egin{aligned} B&=B_0-frac{A-B_0^2}{-2B_0}\ &=frac{A}{2B_0}+frac{B_0}{2} end{aligned} ]向上倍增即可
-
多项式exp
[Bequiv e^A pmod {x^n} ]并没有弄清楚这个咋定义的，因此还不会直接用。
[ln B=A ]构造
[G(B)=ln B-A ]搞
[egin{aligned} B&=B_0-frac{ln B_0-A}{frac{1}{B_0}}\ &=B_0(1-ln B_0+A) end{aligned} ]
多项式取ln
-
多项式积分和求导
[egin{aligned} F(x)&=sum_{i=0}^na_ix^i\ frac{mathrm dF}{mathrm dx}&=sum_{i=0}^{n-1}a_{i+1}(i+1)x^i\ int F(x)&=sum_{i=1}^nfrac{a_{i-1}}{i}x^i end{aligned} ] -
多项式取ln
[egin{aligned} B&=ln Apmod {x^n}\ B'&=ln'(A)A'pmod {x^n}\ &=frac{A'}{A}pmod {x^n} end{aligned} ] -
多项式快速幂
[F^k=e^{kln n}pmod {x^n} ]常数无敌的(nlog n)
拉格朗日反演和快速插值什么的，我还没学...
放一个exp的代码
提几个细节
- 不知道为什么递归版效率优于非递归版，而递归版又好写一点，因此以后就写递归版了
- 每个需要倍增的，每次要外面开两个自己的数组，不许改传进来的值，注意这两个数组一定需要清空
- 要注意一下每次NTT的长度什么的
- 不要封装卷积乘法，因为慢，可以把多个多项式乘起来全部先转换成点表达，然后一起乘，再转回系数，这样需要注意长度。
#include <cstdio>
#include <algorithm>
const int N=(1<<18)+10;
const int mod=998244353,Gi=332748118;
#define mul(a,b) (1ll*(a)*(b)%mod)
#define add(a,b) ((a+b)%mod)
int invA[N],invB[N],lnA[N],lnB[N],exA[N],exB[N],ans[N],f[N],turn[N],n;
int qp(int d,int k){int f=1;while(k){if(k&1)f=mul(f,d);d=mul(d,d),k>>=1;}return f;}
void NTT(int *a,int typ,int len)
{
int L=-1;for(int i=1;i<len;i<<=1) ++L;
for(int i=1;i<len;i++)
{
turn[i]=turn[i>>1]>>1|(i&1)<<L;
if(i<turn[i]) std::swap(a[i],a[turn[i]]);
}
for(int le=1;le<len;le<<=1)
{
int wn=qp(typ?3:Gi,(mod-1)/(le<<1));
for(int p=0;p<len;p+=le<<1)
{
int w=1;
for(int i=p;i<p+le;i++,w=mul(w,wn))
{
int tx=a[i],ty=mul(w,a[i+le]);
a[i]=add(tx,ty);
a[i+le]=add(tx,mod-ty);
}
}
}
if(!typ)
{
int inv=qp(len,mod-2);
for(int i=0;i<len;i++) a[i]=mul(a[i],inv);
}
}
void polyinv(int *a,int *b,int len)
{
if(len==1) {b[0]=qp(a[0],mod-2);return;}
polyinv(a,b,len>>1);
for(int i=0;i<len<<1;i++) invA[i]=invB[i]=0;//
for(int i=0;i<len;i++) invA[i]=a[i],invB[i]=b[i];
NTT(invA,1,len<<1),NTT(invB,1,len<<1);
for(int i=0;i<len<<1;i++) invA[i]=mul(invB[i],add(2,mod-mul(invA[i],invB[i])));
NTT(invA,0,len<<1);
for(int i=0;i<len;i++) b[i]=invA[i];
}
void direv(int *a,int n){for(int i=0;i<n;i++)a[i]=mul(a[i+1],i+1);a[n]=0;}
void inter(int *a,int n){for(int i=n;i;i--)a[i]=mul(a[i-1],qp(i,mod-2));a[0]=0;}
void polyln(int *a,int *b,int len)
{
for(int i=0;i<len;i++) lnA[i]=a[i],lnB[i]=b[i];
polyinv(lnA,lnB,len);
direv(lnA,len-1);
NTT(lnA,1,len<<1),NTT(lnB,1,len<<1);
for(int i=0;i<len<<1;i++) lnA[i]=mul(lnA[i],lnB[i]);
NTT(lnA,0,len<<1);
inter(lnA,len-1);
for(int i=0;i<len;i++) b[i]=lnA[i];
}
void polyexp(int *a,int *b,int len)
{
if(len==1){b[0]=1;return;}//默认a[0]=0
polyexp(a,b,len>>1);
for(int i=0;i<len<<1;i++) exA[i]=exB[i]=0;//
polyln(b,exA,len);
for(int i=0;i<len;i++) exA[i]=add(a[i]+(i==0),mod-exA[i]);
for(int i=0;i<len;i++) exB[i]=b[i];
NTT(exA,1,len<<1),NTT(exB,1,len<<1);
for(int i=0;i<len<<1;i++) exA[i]=mul(exA[i],exB[i]);
NTT(exA,0,len<<1);
for(int i=0;i<len;i++) b[i]=exA[i];
}
int main()
{
scanf("%d",&n);
for(int i=0;i<n;i++) scanf("%d",f+i);
int len=1;while(len<n) len<<=1;
polyexp(f,ans,len);
for(int i=0;i<n;i++) printf("%d ",ans[i]);
return 0;
}
2018.12.29