NTT

NTT

NTT是一种跑得比FFT快的东西（？）。

元素的幂

考虑有限群G，(a in G)。元素的幂就是a的几次方。

使得(a^d=e)的最小正整数d称为a的阶，记作(d=ord(a))。

显然，a的幂生成的集合S是G的子群。因此，(a^{|G|}=e)。

原根

有个结论：(Z_n^*)存在原根(Leftrightarrow)(n=2,4,p^alpha,2p^alpha)，p为奇素数。

tip：469762049，998244353和1004535809都有原根3。

设对于n，我们找到了原根g。设(g_n=g^{frac{p-1}{n}})。那么：

• (g_n=g^{frac{p-1}{n}})

• (g^n_n=g^{p-1}=1)

• (g_{dn}^{dk}=(g^frac{p-1}{dn})^{dk}=(g^{frac{p-1}{n}})^k=g_n^k)（消去引理）

• ((g_n^k)^2=(g_n^{k+n/2})^2=(g^frac{p-1}{n/2})^k=g_{n/2}^k)（折半引理）

• 求和引理：(sum_{i=0}^{n-1}(g_n^k)^i =left{ egin{align} n,nmid k \ 0,n mid k end{align} ight.)， 可以用类似fft的方法证。

然后我们就可以用原根愉快的做fft了。

inline void plus(int &x, int y){ x=1ll*x*y%mod; }
inline void plus(LL &x, LL y){ x*=y; x%=mod; }
inline void pro(int &x){ if (x<0) x+=mod; }

int fpow(LL a, LL x){
    LL ans=1;
    for (LL base=a; x; x>>=1, plus(base, base))
        if (x&1) plus(ans, base);
    return ans;
}
int inv(int x){ return fpow(x, mod-2); }

void fft(int *a, int l, int flag){
    for (int i=0; i<l; ++i) 
        if (i<rev[i]) swap(a[i], a[rev[i]]);
    LL gn, g, x, y;
    for (int mid=1; mid<l; mid<<=1){  //区间半径 
        gn=fpow(G, (mod-1)/(mid<<1));
        if (flag==-1) gn=inv(gn);
        for (int j=0; j<l; j+=(mid<<1)){ g=1;
            for (int k=j; k<j+mid; ++k, plus(g, gn)){
                x=a[k]; y=g*a[k+mid]%mod;
                a[k]=(x+y)%mod; a[k+mid]=(x-y+mod)%mod; 
            }
        }
    }
}

void ntt(int *a, int *b, int &la, int &lb){  //a=a*b 
    int l=1, bits=0; while (l<=la+lb) l<<=1, ++bits;
    int linv=inv(l);
    for (int i=1; i<l; ++i)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bits-1));
    fft(a, l, 1); fft(b, l, 1);
    for (int i=0; i<l; ++i) a[i]=1ll*a[i]*b[i]%mod;
    fft(a, l, -1); la=la+lb;
    for (int i=0; i<=la; ++i) a[i]=1ll*a[i]*linv%mod;
}

归纳一下写fft的思路。首先，应该明确我们要写的主要部分，是把系数表示转换成点值表示。先推式子，对于一个系数表示，把它拆成偶数和奇数的两个多项式。然后，把当前多项式写成两个子多项式的式子。利用单位根的性质，使得分叉个数为4。注意数组大小要开4n。

为了方便，把ntt给封装起来。

任意模数ntt

#include <cstdio> 
#include <algorithm>
using namespace std;

typedef long long LL;
const int maxn=4e5+5, maxc=1e5+5;  //maxn必须是4倍数组 
const LL p1=998244353, p2=1004535809, p3=469762049, G=3;
int n, m, k, cntc, rev[maxn], mod, P;

inline void plus(int &x, int y){ x=1ll*x*y%mod; }
inline void plus(LL &x, LL y){ x*=y; x%=mod; }
inline void pro(int &x){ if (x<0) x+=mod; }

LL fmul(LL a, LL b, LL mod){
    LL ans=0;
    for (; b; b>>=1, a+=a, a%=mod)
        if (b&1) ans+=a, ans%=mod;
    return ans;
}
int fpow(LL a, LL x, LL mod){
	LL ans=1;
	for (LL base=a; x; x>>=1, (base*=base)%=mod)
		if (x&1) (ans*=base)%=mod;
	return ans;
}
int inv(int x){ return fpow(x, mod-2, mod); }

void fft(int *a, int l, int flag){
	for (int i=0; i<l; ++i) 
		if (i<rev[i]) swap(a[i], a[rev[i]]);
	LL gn, g, x, y;
	for (int mid=1; mid<l; mid<<=1){  //区间半径 
		gn=fpow(G, (mod-1)/(mid<<1), mod);
		if (flag==-1) gn=inv(gn);
		for (int j=0; j<l; j+=(mid<<1)){ g=1;
			for (int k=j; k<j+mid; ++k, plus(g, gn)){
				x=a[k]; y=g*a[k+mid]%mod;
				a[k]=(x+y)%mod; a[k+mid]=(x-y+mod)%mod; 
			}
		}
	}
}

void ntt(int *a, int *b, int &la, int &lb){  //a=a*b 
	int l=1, bits=0; while (l<=la+lb) l<<=1, ++bits;
	int linv=inv(l);
	for (int i=1; i<l; ++i)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(bits-1));
	fft(a, l, 1); fft(b, l, 1);
	for (int i=0; i<l; ++i) a[i]=1ll*a[i]*b[i]%mod;
	fft(a, l, -1);
	for (int i=0; i<=la+lb; ++i) a[i]=1ll*a[i]*linv%mod;
}
int A[maxn], B[maxn], C[3][maxn], D[3][maxn];

int crt(LL c1, LL c2, LL c3){
    static LL invp1=fpow(p1, p2-2, p2), invp2=fpow(p2, p1-2, p1);
    static LL p0=p1*p2, invp0=fpow(p0%p3, p3-2, p3);
    LL c0=fmul(p2*invp2, c1, p0)+fmul(p1*invp1, c2, p0); c0%=p0;
    LL k=(c3+p3-c0%p3)*invp0%p3;
    return (c0%mod+k*(p0%mod))%mod;  //k*p0会爆！
}

int main(){
	scanf("%d%d%d", &n, &m, &P);
	for (int i=0; i<=n; ++i) 
		scanf("%d", &A[i]), C[0][i]=C[1][i]=C[2][i]=A[i];
	for (int i=0; i<=m; ++i) 
		scanf("%d", &B[i]), D[0][i]=D[1][i]=D[2][i]=B[i];
	mod=p1; ntt(C[0], D[0], n, m);
	mod=p2; ntt(C[1], D[1], n, m);
	mod=p3; ntt(C[2], D[2], n, m); mod=P; 
	for (int i=0; i<=n+m; ++i)
		printf("%d ", crt(C[0][i], C[1][i], C[2][i]));
	return 0;
}

两天后的PS：注意对于998244352，1004535808，469762048的分解。(998244352=2^{23}*x)，(1004535808=2^{21}*x)，(469762048=2^{26}*x)。这是因为三个质数本来就是(p=c*2^l+1)的形式。