luogu P4726 多项式指数函数(模板题FFT多项式求逆多项式对数函数)
Posted suncongbo
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了luogu P4726 多项式指数函数(模板题FFT多项式求逆多项式对数函数)相关的知识,希望对你有一定的参考价值。
手动博客搬家: 本文发表于20181127 08:39:42, 原地址https://blog.csdn.net/suncongbo/article/details/84559818
题目链接: https://www.luogu.org/problem/show?pid=4726
题意: 给定(n)次多项式(A(x)) 求多项式(f(x))满足(f(x)equiv e^{A(x)} (mod x^n))
题解
这个比对数函数复杂一些。。
前铺知识
泰勒展开
对于一个函数,我们可以用以下方式用它的高阶导数进行无穷逼近:[f(x)=f(a)+f‘(a)(x-a)+f‘‘(a)frac{(x-a)^2}{2!}+f‘‘‘(a)frac{(x-a)^3}{3!}+...]
牛顿迭代
已知要求的多项式(f)满足(g(f(x))equiv 0(mod x^n)) (g)已知,那么可以通过如下的方式倍增求解:
假设求出了原问题(mod x^n)的答案(f_0(x)), 考虑转移到(f(x) mod x^{2n}).
根据泰勒展开公式: [0=g(f)=g(f_0)+g‘(f_0)(f-f_0)+g‘‘(f_0)frac{(f-f_0)^2}{2!}+...]
但是,由于一个显然的结论——(fequiv f_0(mod x^n)), 因此((f-f_0)^2equiv 0(mod x^{2n})), 因此公式里只需要计算前两项,后面的项都在(mod x^{2n})意义下为(0)!
即(g(f_0)+g‘(f_0)(f-f_0)=0)
(f=f_0-frac{g(f_0)}{g‘(f_0)})
如此即可求解。
(太神了啊太神了啊呜呜呜。。我大概永远也搞不出这么神奇的东西吧。。)
本题题解
根据牛顿迭代的法则,令(g(f)=ln f-A), 则(f=f_0-frac{ln f_0-A}{frac{1}{f_0}}=f_0(1-ln f_0+A))
递归求解即可。
值得注意的是FFT的大小
(A)应该带入(2n)次, (ln f_0)应该带入(2n)次, (f_0)应该带入(n)次, FFT乘法因为是(2n)次乘以(n)次,所以应该取到(4n)个单位根。
FFT这个地方太容易出错了!范围大了常数大好几倍,范围小了就会出错。
时间复杂度为(T(n)=T(frac{n}{2})+O(nlog n)), (T(n)=O(nlog n))
常数我算的应该是(48)倍。每次(n)到(2n)的迭代需要做一次(2n)的(ln)和三次(4n)的DFT/IDFT. 因此(18(2nlog n)+3(4nlog n)=48nlog n).
(飞了……)
代码
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define llong long long
#define ldouble long double
#define uint unsigned int
#define ullong unsigned long long
#define udouble unsigned double
#define uldouble unsigned long double
#define modinc(x) {if(x>=P) x-=P;}
#define pii pair<int,int>
#define piii pair<pair<int,int>,int>
#define piiii pair<pair<int,int>,pair<int,int> >
#define pli pair<llong,int>
#define pll pair<llong,llong>
#define Memset(a,x) {memset(a,x,sizeof(a));}
using namespace std;
const int N = 1<<19;
const int LGN = 19;
const int G = 3;
const int P = 998244353;
llong tmp1[N+3],tmp2[N+3],tmp3[N+3],tmp4[N+3]; //inv
llong tmp5[N+3],tmp6[N+3],tmp7[N+3],tmp8[N+3]; //ln
llong tmp9[N+3],tmp10[N+3],tmp11[N+3],tmp12[N+3]; //exp
llong a[N+3],b[N+3];
int id[N+3];
int n;
llong quickpow(llong x,llong y)
{
llong cur = x,ret = 1ll;
for(int i=0; y; i++)
{
if(y&(1ll<<i))
{
ret = ret*cur%P;
y-=(1ll<<i);
}
cur = cur*cur%P;
}
return ret;
}
llong mulinv(llong x) {return quickpow(x,P-2);}
void initid(int dgr)
{
int len = 0;
for(int i=0; i<=LGN; i++) if(dgr==(1<<i)) {len = i; break;}
id[0] = 0;
for(int i=1; i<dgr; i++) id[i] = (id[i>>1]>>1)|((i&1)<<(len-1));
}
void ntt(int dgr,int coe,llong poly[],llong ret[])
{
initid(dgr);
for(int i=0; i<dgr; i++) ret[i] = poly[i];
for(int i=0; i<dgr; i++) if(i<id[i]) swap(ret[i],ret[id[i]]);
for(int i=1; i<=(dgr>>1); i<<=1)
{
llong tmp = quickpow(G,(P-1)/(i<<1));
if(coe==-1) tmp = mulinv(tmp);
for(int j=0; j<dgr; j+=(i<<1))
{
llong expn = 1ll;
for(int k=0; k<i; k++)
{
llong x = ret[j+k],y = ret[j+k+i]*expn%P;
ret[j+k] = x+y; modinc(ret[j+k]);
ret[j+i+k] = x-y+P; modinc(ret[j+i+k]);
expn = (expn*tmp)%P;
}
}
}
if(coe==-1)
{
llong tmp = mulinv(dgr);
for(int j=0; j<dgr; j++) ret[j] = ret[j]*tmp%P;
}
}
void polyinv(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<dgr; i++) ret[i] = 0ll;
ret[0] = mulinv(poly[0]);
for(int i=1; i<=(dgr>>1); i<<=1)
{
for(int j=0; j<(i<<2); j++) tmp1[j] = j<i ? ret[j] : 0ll;
for(int j=0; j<(i<<2); j++) tmp2[j] = j<(i<<1) ? poly[j] : 0ll;
ntt((i<<2),1,tmp1,tmp3); ntt((i<<2),1,tmp2,tmp4);
for(int j=0; j<(i<<2); j++) tmp3[j] = tmp3[j]*tmp3[j]%P*tmp4[j]%P;
ntt((i<<2),-1,tmp3,tmp4);
for(int j=0; j<(i<<1); j++) ret[j] = (tmp1[j]+tmp1[j]-tmp4[j]+P)%P;
}
for(int i=dgr; i<(dgr<<1); i++) ret[i] = 0ll;
}
void polyder(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<dgr-1; i++) ret[i] = poly[i+1]*(i+1)%P;
}
void polyint(int dgr,llong poly[],llong ret[])
{
for(int i=1; i<dgr; i++) ret[i] = poly[i-1]*mulinv(i)%P;
}
void polyln(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<dgr; i++) ret[i] = 0ll;
polyder(dgr,poly,tmp5);
polyinv(dgr,poly,tmp6);
ntt((dgr<<1),1,tmp5,tmp7); ntt((dgr<<1),1,tmp6,tmp8);
for(int i=0; i<(dgr<<1); i++) tmp7[i] = tmp7[i]*tmp8[i]%P;
ntt((dgr<<1),-1,tmp7,tmp8);
polyint(dgr,tmp8,ret);
}
void polyexp(int dgr,llong poly[],llong ret[])
{
for(int i=0; i<dgr; i++) ret[i] = i==0;
for(int i=1; i<=(dgr>>1); i<<=1)
{
for(int j=0; j<(i<<2); j++) tmp9[j] = j>=(i<<1) ? 0ll : ((j==0)+poly[j])%P;
for(int j=0; j<(i<<2); j++) tmp10[j] = j>=i ? 0ll : ret[j];
polyln((i<<1),tmp10,tmp11);
for(int j=0; j<(i<<1); j++) {tmp9[j] = tmp9[j]-tmp11[j]+P; modinc(tmp9[j]);}
ntt((i<<2),1,tmp10,tmp12); ntt((i<<2),1,tmp9,tmp11);
for(int j=0; j<(i<<2); j++) tmp12[j] = tmp12[j]*tmp11[j]%P;
ntt((i<<2),-1,tmp12,tmp11);
for(int j=0; j<(i<<1); j++) ret[j] = tmp11[j];
}
}
int main()
{
scanf("%d",&n); int dgr = 1; while(dgr<=n) dgr<<=1;
for(int i=0; i<n; i++) scanf("%lld",&a[i]);
polyexp(dgr,a,b);
for(int i=0; i<n; i++) printf("%lld ",b[i]);
return 0;
}
以上是关于luogu P4726 多项式指数函数(模板题FFT多项式求逆多项式对数函数)的主要内容,如果未能解决你的问题,请参考以下文章