luogu P4238 多项式求逆 (模板题、FFT)

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

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

题意: 给定(n)次多项式(A(x)), 求(n)次多项式(B(x))满足(B(x)A(x)equiv 1(mod x^n))

题解:
DFT，每个数对(998244353)求逆元。IDFT回来。
发现，错了。
为什么呢？
因为要对(x^n)取模。
例如，(1-x)在模(x^2)意义下的逆元是(1+x), 但是在实际上逆元是(1+x+x^2+x^3+...), 是无穷和式。
所以此路不通。

考虑求解多项式问题的常用方法——分治法。
设已求(B_0(x))满足(B_0(x)A(x)equiv 1(mod x^n)), 现要求(B(x))满足(B(x)A(x)equiv 1(mod x^{2n}))
显然有(B(x)-B_0(x)equiv 0(mod x^n))
为了凑出(x^{2n})两边平方得
(B^2(x)-2B_0(x)B(x)+B_0^2(x)equiv 0(mod x^{2n}))
如何求(B)呢？因为(A(x)B(x)equiv 0(mod x^{2n})), 我们将式子两边同乘(A(x))
(A(x)B(x)B(x)-2B_0(x)A(x)B(x)+A(x)B_0^2(x)equiv 0(mod x^{2n}))
(B(x)equiv 2B_0(x)-A(x)B_0^2(x) (mod x^{2n}))
右边的式子FFT计算即可。
时间复杂度？
(T(n)=T(frac{n}{2})+O(nlog n))
(T(n)=O(nlog n)).
常数？首先隐藏在递归复杂度背后有一个(2)倍常数。
然后我们把两个多项式相乘需要(3)次FFT, 三个就要(6)次。
因此常数为(12)倍。
如何优化？
(IDFT(DFT(IDFT(DFT(A) imes DFT(B_0))) imes DFT(B_0)))
变成(IDFT(DFT(A) imes DFT^2(B_0))
(3)次即可！
常数变为(6)倍。
UPD: 关于这里的常数问题: 因为我递归里DFT的范围是(2n)，最终的复杂度是(T(2n) = 6(2n)log (2n)), 因此我认为应为(12)倍常数。
空间？空间复杂度(O(n)), 但是要开(4d)的数组，其中(d)是(>n)的最小的(2)的幂。

代码
因为FFT数组清零等原因代码(对我来说)很难写。
贴一下我刚刚写的版本吧，还算是比较简单。
(话说怎么CSDN突然傲娇了啊。。缩进变成1格？？)

#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],tmp5[N+3],tmp6[N+3];
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;
            }
        }
    }
}

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++) tmp5[j] = tmp3[j]*tmp3[j]%P*tmp4[j]%P;
        ntt((i<<2),-1,tmp5,tmp6); llong tmp = mulinv(i<<2);
        for(int j=0; j<(i<<2); j++) tmp6[j] = tmp6[j]*tmp%P;
        for(int j=0; j<(i<<1); j++) ret[j] = (tmp1[j]+tmp1[j]-tmp6[j]+P)%P;
    }
}

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