存一份MTT多项式求逆

https://www.luogu.org/problemnew/show/P4239

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define Rep(i,s,t) for(register int i=s;i<t;++i)
typedef double db;
typedef long long ll;
const int N=2333333,mod=1000000007,m=sqrt(mod);
const double pi=acos(-1);
inline int fp(int a,int b){
    int s=1;
    while(b){
        if(b&1)s=1ll*s*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return s;
}
struct cp{
    db a,b;
    friend inline cp operator+(cp a,cp b){return (cp){a.a+b.a,a.b+b.b};}
    friend inline cp operator-(cp a,cp b){return (cp){a.a-b.a,a.b-b.b};}
    friend inline cp operator*(cp a,cp b){return (cp){a.a*b.a-a.b*b.b,a.a*b.b+a.b*b.a};}
}W[N],A[N],B[N],C[N],D[N];
int r[N];
inline void FFT(cp *P,int n,int opt){
    Rep(i,0,n)
        if(i<r[i])
            swap(P[i],P[r[i]]);
    for(register int i=1;i<n;i<<=1)
        for(register int p=i<<1,j=0;j<n;j+=p)
            Rep(k,0,i){
                cp w=(cp){W[n/i*k].a,W[n/i*k].b*opt};
                cp X=P[j+k],Y=P[i+j+k]*w;
                P[j+k]=X+Y;P[i+j+k]=X-Y;
            }
    if(opt==-1)
        Rep(i,0,n)
            P[i].a/=1.0*n;
}
inline void Multi(int *a,int *b,int n,int *ret){
    Rep(i,0,n+n)
        A[i]=B[i]=C[i]=D[i]=(cp){0,0};
    Rep(i,0,n){
        a[i]%=mod;b[i]%=mod;
        A[i]=(cp){(a[i]/m)*1.0,0};
        B[i]=(cp){(a[i]%m)*1.0,0};
        C[i]=(cp){(b[i]/m)*1.0,0};
        D[i]=(cp){(b[i]%m)*1.0,0};
    }
    int N,l=0;
    for(N=1;N<=n;N<<=1)++l;
    Rep(i,0,N)
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(register int i=1;i<N;i<<=1)
        Rep(k,0,i)
            W[N/i*k]=(cp){cos(k*pi/i),sin(k*pi/i)};
    FFT(A,N,1);FFT(B,N,1);FFT(C,N,1);FFT(D,N,1);
    Rep(i,0,N){
        cp tmp=A[i]*C[i];
        C[i]=B[i]*C[i],B[i]=B[i]*D[i],D[i]=D[i]*A[i];
        A[i]=tmp;C[i]=C[i]+D[i];
    }
    FFT(A,N,-1);FFT(B,N,-1);FFT(C,N,-1);
    Rep(i,0,n){
        ret[i]=0;
        ret[i]=(ret[i]+(ll)(A[i].a+0.5)%mod*m%mod*m%mod)%mod;
        ret[i]=(ret[i]+(ll)(C[i].a+0.5)%mod*m%mod)%mod;
        ret[i]=(ret[i]+(ll)(B[i].a+0.5)%mod)%mod;
        ret[i]=(ret[i]+mod)%mod;
    }
}
int n,a[N],b[N],c[N],d[N];
inline void Inv(int *a,int *b,int n){
    if(n==1){
        b[0]=fp(a[0],mod-2);
        return;
    }
    Inv(a,b,n>>1);
    Multi(a,b,n,c);
    Multi(c,b,n,d);
    Rep(i,0,n)
        b[i]=(b[i]+b[i])%mod;
    Rep(i,0,n)
        b[i]=(b[i]+mod-d[i])%mod;
}
int main(){
    scanf("%d",&n);int N;
    Rep(i,0,n)
        scanf("%d",a+i);
    for(N=1;N<n;N<<=1);
    Inv(a,b,N);
    Rep(i,0,n)
        printf("%d ",b[i]);
    return 0;
}