NTT
NTT
NTT是一種跑得比FFT快的東西(?)。
元素的冪
考慮有限群G,\(a \in G\)。元素的冪就是a的幾次方。
使得\(a^d=e\)的最小正整數(shù)d稱為a的階,記作\(d=ord(a)\)。
顯然,a的冪生成的集合S是G的子群。因此,\(a^{|G|}=e\)。
原根
有個(gè)結(jié)論:\(Z_n^*\)存在原根\(\Leftrightarrow\)\(n=2,4,p^\alpha,2p^\alpha\),p為奇素?cái)?shù)。
tip:469762049,998244353和1004535809都有原根3。
設(shè)對(duì)于n,我們找到了原根g。設(shè)\(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){ //區(qū)間半徑 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的思路。首先,應(yīng)該明確我們要寫的主要部分,是把系數(shù)表示轉(zhuǎn)換成點(diǎn)值表示。先推式子,對(duì)于一個(gè)系數(shù)表示,把它拆成偶數(shù)和奇數(shù)的兩個(gè)多項(xiàng)式。然后,把當(dāng)前多項(xiàng)式寫成兩個(gè)子多項(xiàng)式的式子。利用單位根的性質(zhì),使得分叉?zhèn)€數(shù)為4。注意數(shù)組大小要開4n。
為了方便,把ntt給封裝起來。
任意模數(shù)ntt
#include <cstdio> #include <algorithm> using namespace std;typedef long long LL; const int maxn=4e5+5, maxc=1e5+5; //maxn必須是4倍數(shù)組 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){ //區(qū)間半徑 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會(huì)爆! }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:注意對(duì)于998244352,1004535808,469762048的分解。\(998244352=2^{23}*x\),\(1004535808=2^{21}*x\),\(469762048=2^{26}*x\)。這是因?yàn)槿齻€(gè)質(zhì)數(shù)本來就是\(p=c*2^l+1\)的形式。
轉(zhuǎn)載于:https://www.cnblogs.com/MyNameIsPc/p/9592422.html
總結(jié)
- 上一篇: Python 33(1) UDP协议
- 下一篇: CMA内存管理子系统