COGS 2294 [HZOI 2015] 释迦（NTT mod any prime）

题目描述：
找不到传送门。。。
给两个次数界为n的多项式，求这两个多项式的乘积，输出x的0次项到n-1次项的系数 mod 23333333

分析：
%真正的dalao

NTT只能求在特定模数下 的值
对于任意模数，我们可以选择三个FFT模数分别做NTT，最后用CRT合并

一次卷积后每个数可以达到$10^{23}(n*mod^2)$左右，所以我们需要三个乘积大于$10^{23}$的FFT模数
合并的时候$10^{23}>2^{64}$，所以我们不能直接用CRT，需要一点小技巧

$$ans≡a_1~(mod~m_1)$$

$$ans≡a_2~(mod~m_2)$$

$$ans≡a_3~(mod~m_3)$$

先用CRT合并$a_1,a_2$，得到$ans≡A~~(mod~M)$
其中$A=a_1*m_2*m_2^{-1}+a_2*m_1*m_1^{-1},M=m_1*m_2$
$ans=k*M+A=x*m_3+a_3$

$$k*M≡a_3-A~~~(mod~m_3)$$

$$k≡(a_3-A)*M^{-1}~~~(mod~m_3)$$

求出k后带入$ans=k*M+A$，然后对23333333取模即可

中间做乘法的过程需要用快速乘，否则会爆ll

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm> 
#define ll long long

using namespace std;

const ll mod=23333333;
const int N=500005;
ll M[10],f[4][N],g[4][N],ans[N];
int n,m,fn;

ll KSM(ll a,ll b,ll p) {
    ll t=1;
    a%=p;          //防止爆ll 
    while (b) {
        if (b&1) t=(t*a)%p;
        b>>=1;
        a=(a*a)%p;
    }
    return t%p;
}

void NTT(int n,ll *a,int opt,ll p) {
    int i,j=0,k;
    for (i=0;i<n;i++) {
        if (i>j) swap(a[i],a[j]);
        for (int l=n>>1;(j^=l)<l;l>>=1);
    }
    for (i=2;i<=n;i<<=1) {
        int m=i>>1;
        ll wn=KSM(3,(p-1)/i,p);
        for (j=0;j<n;j+=i) {
            ll w=1;
            for (k=0;k<m;k++,w=(w*wn)%p) {
                ll z=(a[j+k+m]*w)%p;
                a[j+k+m]=((a[j+k]-z)%p+p)%p;
                a[j+k]=(a[j+k]+z)%p;
            }
        }
    }
    if (opt==-1) reverse(a+1,a+n);
} 

void solve(ll *a,ll *b,ll p) {
    NTT(fn,a,1,p);
    NTT(fn,b,1,p);
    for (int i=0;i<fn;i++) a[i]=(a[i]*b[i])%p;
    NTT(fn,a,-1,p);
    ll inv=KSM(fn,p-2,p);
    for (int i=0;i<fn;i++) a[i]=(a[i]*inv)%p;
}

ll mul(ll num,ll x,ll p) {     //num*x
    ll ans=0;
    num%=p;
    while (x) {
        if (x&1) ans=(ans+num)%p;
        num=(num+num)%p;
        x>>=1;
    }
    return ans%p;
}

ll china(ll a1,ll a2,ll a3) {
    ll MM=M[1]*M[2];
    ll x,y;
    x=KSM(M[2],M[1]-2,M[1]);     //CRT
    y=KSM(M[1],M[2]-2,M[2]);     
    ll A=( mul( a1*M[2]%MM, x%MM, MM )%MM + mul( a2*M[1]%MM, y%MM, MM)%MM )%MM;
    ll k=((a3-A)%M[3]+M[3])%M[3];
    k=((k*KSM(MM,M[3]-2,M[3]))%M[3]+M[3])%M[3];
    return (A%mod+(k%mod*MM%mod)%mod)%mod;
}

int main()
{
    scanf("%d",&n);
    M[1]=998244353; M[2]=1004535809; M[3]=469762049;
    for (int i=0;i<n;i++) scanf("%lld",&f[1][i]),f[2][i]=f[3][i]=f[1][i];
    for (int i=0;i<n;i++) scanf("%lld",&g[1][i]),g[2][i]=g[3][i]=g[1][i];
    fn=1;
    while (fn<=n*2) fn<<=1;
    for (int i=1;i<=3;i++) solve(f[i],g[i],M[i]);
    for (int i=0;i<n;i++) printf("%lld\n",china(f[1][i],f[2][i],f[3][i]));
    return 0;
}