luogu P4725 多项式对数函数 (模板题FFT多项式求逆求导和积分)

Posted 2021-02-05 suncongbo

手动博客搬家: 本文发表于20181125 13:25:03, 原地址https://blog.csdn.net/suncongbo/article/details/84487306

题目链接: https://www.luogu.org/problemnew/show/P4725

题目大意: 给定一个(n)次多项式(A(x)), 求一个(n)次多项式(B(x))满足(B(x)equiv ln A(x) (mod x^n))

题解: 神数学模板题……
数学真奇妙！
前驱知识
导数、积分相关
幂函数的求导
(f(x)=x^n, f‘(x)=nx^{n-1})
和的导数等于导数的和
((f+g)‘(x)=f‘(x)+g‘(x))
一般多项式的求导
(f(x)=sum^{n-1}_{i=0} a_ix^i, f‘(x)=sum^{n-2}_{i=0} (i+1)a_{i+1}x^i)
对数函数(ln)的求导
(f(x)=ln(x), f‘(x)=frac{1}{x})
复合函数求导——链式法则
(f(g(x))‘=f‘(g(x))g‘(x))
求导的逆运算——积分

本题解法
(g(x)equiv ln f(x) (mod x^n))
两边同时求导可得
(g‘(x)equiv frac{f‘(x)}{f(x)} (mod x^n))
结束！
多项式求逆算(frac{1}{f(x)})，再和(f‘(x))相乘即可得到(g‘(x))。
(多项式求逆见蒟蒻一篇博客 https://blog.csdn.net/suncongbo/article/details/84485718)
(g‘(x))积个分得到(g(x)). 常数项，直接为(0).
时间复杂度(O(nlog n))
常数，我写的大概(9)倍吧，求逆是(6)倍，再做个乘法就是(3)倍。
UPD: 仔细想了一下我这个常数好像是(18)倍（见我多项式求逆那篇博客）

代码

#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 P = 998244353;
const int LGN = 19;
const int G = 3;
llong a[N+3];
llong b[N+3];
llong tmp1[N+3],tmp2[N+3],tmp3[N+3],tmp4[N+3]; //inv
llong tmp7[N+3],tmp8[N+3],tmp9[N+3],tmp10[N+3]; //ln
int id[N+2];
int n;

void initid(int _len)
{
 id[0] = 0;
 for(int i=1; i<(1<<_len); i++) id[i] = (id[i>>1]>>1)|((i&1)<<(_len-1));
}

llong quickpow(llong x,llong y)
{
 llong cur = x,ret = 1ll;
 for(int i=0; y; i++)
 {
  if(y&(1ll<<i))
  {
   y-=(1ll<<i); ret = ret*cur%P;
  }
  cur = cur*cur%P;
 }
 return ret;
}
llong mulinv(llong x) {return quickpow(x,P-2);}

void ntt(int dgr,int coe,llong poly[],llong ret[])
{
 int len = 0; for(int i=0; i<=LGN; i++) if((1<<i)==dgr) {len = i; break;}
 initid(len); for(int i=0; i<dgr; i++) ret[i] = 0ll;
 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 = (expn*ret[j+i+k])%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 i=0; i<dgr; i++) ret[i] = ret[i]*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[])
{
 polyder(dgr,poly,tmp7);
 polyinv(dgr,poly,tmp8);
 ntt((dgr<<1),1,tmp8,tmp9); ntt((dgr<<1),1,tmp7,tmp10);
 for(int i=0; i<(dgr<<1); i++) tmp9[i] = tmp9[i]*tmp10[i]%P;
 ntt((dgr<<1),-1,tmp9,tmp10);
 polyint(dgr,tmp10,ret);
}

int main()
{
 scanf("%d",&n); int dgr = 1; while(dgr<=n) dgr<<=1;
 for(int i=0; i<n; i++) scanf("%lld",&a[i]);
 polyln(dgr,a,b);
 for(int i=0; i<n; i++) printf("%lld ",b[i]);
 return 0;
}