分治FFT&&多项式求逆学习笔记

分治FFT：

解决的是形似以下的问题：

给定n次多项式$g(x)$，求多项式$f(x)$，其中$f$的第$i$项系数的表达式为$\sum f_j\times g_{i-j}$。

解法：

不难发现式子也是卷积的形式，但是与普通多项式乘法不一样的是，每一项的系数依赖前面的项的系数，使得普通的FFT无法起作用。

考虑分治，将区间$[l,r]$分为$[l,mid],[mid+1,r]$两个区间计算，计算完$[l,mid]$中的多项式的系数之后，可以很方便的将$[l,mid]$对$[mid+1,r]$的贡献给算出来，贡献的计算同上面的递推式。

由于贡献的计算也需要依赖前面的系数，所以分治时的进程是：
divide(l,mid);
calc([l,mid]->[mid+1,r]);
divide(mid+1,r);

不难发现每个点之前的部分会被拆成若干个区间累加到这个位置上，所以正确性显而易见。

/*=======================================
 * Author : ylsoi
 * Time : 2019.1.29
 * Problem : luogu4721
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
    freopen("luogu4721.in","r",stdin);
    freopen("luogu4721.out","w",stdout);
}

template<typename T>void read(T &_){
    _=0; T fl=1; char ch=getchar();
    for(;!isdigit(ch);ch=getchar())if(ch=='-')fl=-1;
    for(;isdigit(ch);ch=getchar())_=(_<<1)+(_<<3)+(ch^'0');
    _*=fl;
}

const int maxn=1e5+10;
const int mod=998244353;

ll qpow(ll x,ll y){
    ll ret=1; x%=mod;
    while(y){
        if(y&1)ret=ret*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return ret;
}

struct ksslbh{
    int lim,cnt,dn[maxn<<2];
    ll g[maxn<<2],ig[maxn<<2];
    void ntt(ll *A,int ty){
        REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
        for(int len=1;len<lim;len<<=1){
            ll w= ty==1 ? g[len<<1] : ig[len<<1];
            for(int L=0;L<lim;L+=len<<1){
                ll wk=1;
                REP(i,L,L+len-1){
                    ll u=A[i],v=A[i+len]*wk%mod;
                    A[i]=(u+v)%mod;
                    A[i+len]=(u-v)%mod;
                    wk=wk*w%mod;
                }
            }
        }
    }
    void work(int n,int m,ll *a,ll *b){
        lim=1,cnt=0;
        while(lim<=n+m)lim<<=1,++cnt;
        if(!cnt)cnt=1;
        REP(i,0,lim){
            dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
            if(i>n)a[i]=0;
            if(i>m)b[i]=0;
        }
        g[lim]=qpow(3,(mod-1)/lim);
        ig[lim]=qpow(g[lim],mod-2);
        for(int i=lim>>1;i;i>>=1){
            g[i]=g[i<<1]*g[i<<1]%mod;
            ig[i]=ig[i<<1]*ig[i<<1]%mod;
        }
        ntt(a,1),ntt(b,1);
        REP(i,0,lim-1)a[i]=a[i]*b[i]%mod;
        ntt(a,-1);
        ll inv=qpow(lim,mod-2);
        REP(i,0,lim-1)a[i]=a[i]*inv%mod;
    }
}T;

int n;
ll f[maxn<<2],g[maxn<<2],a[maxn<<2],b[maxn<<2];

void divide(int l,int r){
    if(l==r)return;
    int mid=(l+r)>>1;
    divide(l,mid);
    int nn=mid-l,mm=r-l-1;
    //nn means the pows of f/a
    //mm means the pows of g/b
    REP(i,0,nn)a[i]=f[l+i];
    REP(i,0,mm)b[i]=g[i+1];
    T.work(nn,mm,a,b);
    REP(i,mid+1,r)f[i]=(f[i]+a[i-1-l])%mod;
    divide(mid+1,r);
}

int main(){
    //File();
    read(n);
    REP(i,1,n-1)read(g[i]);
    f[0]=1;
    divide(0,n);
    REP(i,0,n-1)printf("%lld ",(f[i]+mod)%mod);
    return 0;
}

多项式求逆：

对于一个多项式$A(x)$，找到一个多项式$A^{-1}(x)$，使得$A(x)\times A^{-1}(x)$只有常数项的系数为1，其余项的系数全部为0，那么我们称$A^{-1}(x)$为多项式$A(x)$的逆。

当常数项存在逆元的时候，利用递推的方法去求解多项式的逆，不难发现所得的多项式是一个无限项的级数，所以我们一般只取前面$n-1$次，即对$x^n$取模。

暴力求法：

从低次依次向高次递推，每次求解$n$次项的系数$a_n$时，将前$n$项带入，所得的和需为$0$，直接解出来即可。

优化：

观察上面的递推过程，不难发现后面的项求解需要用到前面的项的信息，利用分治FFT可以优化到$n{\log }^{2}n$。

更加优秀的解法：

直接一项一项地递推显然太浪费时间，考虑利用倍增的思想来利用前面的项推出后面的项。

假设我们现在得到了前$\frac{n}{2}$项的系数，那么我们就要想办法得出前$n$项的系数。

假设前$\frac{n}{2}$项的多项式为$B(x)$，最终的多项式为$A^{-1}(x)$，那么我们可以得到
$$A^{-1}(x)-B(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\frac{n}{2}}$$
为了使模数提升至$x^n$，可以将等式两边平方得：
$$(A^{-1}(x)-B(x){)}^2\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$$
考虑证明这一步，任意一个$i(i<n)$次项的系数的拆分必定有一个数$<\frac{n}{2}$，而任意$<\frac{n}{2}$次项的系数都为0，所以$i(i<n)$次项的系数都为$0$。

同时注意这种变化只有在等号右边等于$0$时才成立，这也是为什么要拿两式相减的理由。

然后将平方项拆开：
$$A^{-2}(x)-2\times A^{-1}(x)\times B(x)+B^2(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$$
两边同时乘以$A(x)$：
$$\begin{gathered} A^{-1}(x)-2\times B(x)+A(x)\times B^2(x)\equiv 0\mathrm{m}\mathrm{o}\mathrm{d}{x}^n \\ {A}^{-1}(x)\equiv 2\times B(x)-A(x)\times B^2(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n \end{gathered}$$
然后我们就成功地得到了在模$x^n$意义下的多项式$A^{-1}(x)$。

具体实现：

考虑递归实现上述过程。

实现过程中，理论上$A(x)$和$B(x)$项数都很多，但为了保证复杂度，$A(x)$每次取对计算有贡献的前$n$项，而$B(x)$原本是在模$x^{\frac{n}{2}}$意义下的多项式，只取前$\frac{n}{2}$项即可，所以最后乘出来的多项式为$2n$次的，时间复杂度满足$T(n)=n\log n+T(\frac{n}{2})$。

所以最终时间复杂度为$n\log n$。

 /*=======================================
 * Author : ylsoi
 * Time : 2019.2.2
 * Problem : luogu4238
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
    freopen("luogu4238.in","r",stdin);
    freopen("luogu4238.out","w",stdout);
}

template<typename T>void read(T &_){
	_=0; T fl=1; char ch=getchar();
	for(;!isdigit(ch);ch=getchar())if(ch=='-')fl=-1;
	for(;isdigit(ch);ch=getchar())_=(_<<1)+(_<<3)+(ch^'0');
	_*=fl;
}

const int maxn=1e5+10;
const int mod=998244353;

int lim,cnt,dn[maxn<<2];
ll g[maxn<<2],ig[maxn<<2];

ll qpow(ll x,ll y){
	ll ret=1; x%=mod;
	while(y){
		if(y&1)ret=ret*x%mod;
		x=x*x%mod;
		y>>=1;
	}
	return ret;
}

void ntt(ll *A,int ty){
	REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
	for(int len=1;len<lim;len<<=1){
		ll w= ty==1 ? g[len<<1] : ig[len<<1];
		for(int L=0;L<lim;L+=len<<1){
			ll wk=1;
			REP(i,L,L+len-1){
				ll u=A[i],v=A[i+len]*wk%mod;
				A[i]=(u+v)%mod;
				A[i+len]=(u-v)%mod;
				wk=wk*w%mod;
			}
		}
	}
}

int n;
ll f[maxn<<2],a[maxn<<2],b[maxn<<2];

void init(){
	read(n);
	REP(i,0,n-1)read(f[i]);
	lim=1;
	while(lim<=n+n)lim<<=1;
	g[lim]=qpow(3,(mod-1)/lim);
	ig[lim]=qpow(g[lim],mod-2);
	for(int i=lim>>1;i;i>>=1){
		g[i]=g[i<<1]*g[i<<1]%mod;
		ig[i]=ig[i<<1]*ig[i<<1]%mod;
	}
}

void poly_inv(int k){
	if(k==1){
		b[0]=qpow(f[0],mod-2);
		return;
	}
	poly_inv((k+1)>>1);
	lim=1,cnt=0;
	while(lim<=k+k)lim<<=1,++cnt;
	if(!cnt)cnt=1;
	REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
	REP(i,0,k-1)a[i]=f[i];
	REP(i,k,lim-1)a[i]=0;
	ntt(a,1),ntt(b,1);
	REP(i,0,lim-1)b[i]=b[i]*(2ll-a[i]*b[i]%mod)%mod;
	ntt(b,-1);
	ll inv=qpow(lim,mod-2);
	REP(i,0,k-1)b[i]=b[i]*inv%mod;
	REP(i,k,lim)b[i]=0;
}

int main(){
    File();
	init();
	poly_inv(n);
	REP(i,0,n-1)printf("%lld ",(b[i]%mod+mod)%mod);
    return 0;
}