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\{ \begin{align} n,n\mid k \\ 0,n \nmid k \end{align} \right.\)， 可以用類似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\)的形式。

轉載于:https://www.cnblogs.com/MyNameIsPc/p/9592422.html

總結

以上是生活随笔為你收集整理的NTT的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。