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!}+...$$

但是，由于一个显然的结论——$f\equiv f_0(\mod x^n)$, 因此$(f-f_0)^2\equiv 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(n\log n)$, $T(n)=O(n\log n)$

常数我算的应该是$48$倍。每次$n$到$2n$的迭代需要做一次$2n$的$\ln$和三次$4n$的DFT/IDFT. 因此$18(2n\log n)+3(4n\log n)=48n\log 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、多项式求逆、多项式对数函数）

题解

代码

相关文章