快速数论变换(FNT)是环 $\mathbb{Z}/ m \mathbb{Z}$ 上的 Fourier 变换(FFT)。
至于 快速 Fourier 变换是怎样的有什么用处,这里就不多说了,可参考 这里。
NFT 的核心问题
无论是 FNT 还是 FFT 其本质其关键就是寻找一个 $w$ 使得 $w^{2^n} = 1$。在复数域中这个问题是显然的,而在一个环那就不那么简单了,这里我们考虑环 $R = \mathbb{Z}/ m \mathbb{Z}$。$R$ 为域当且仅当 $m$ 为素数。我们的问题是:
- 我们选取 $R$ 中找比较大(满足我们的需求)的 $n$ 使得 $w^{2^n} = 1$ ;
- 找出对应的“原根“;
- 类似 FFT 的处理
- 应用时各种可能出错的情形,最常见的是溢出,还有只适用于数据范围不超过 P 的非负整数。
NFT 问题解决
为使得分析问题更为简单,我们考虑在 $m = p$ 为素数的情形,此时,我们有 $2^n \mid p-1$ 即 $p = k \cdot 2^n + 1$ 为(Fermat)素数,例如:
- $p=479 \times 2^{21} +1 = 1004535809,g = 3$
- $p= 13 \times 2^{20} + 1 = 13631489,g = 15$
- $p= 17 \times 2^{27} + 1 = 2281701377,g=3$
更多常数选择可见 这里。
最终我们选择了 $FM = 1004535809$ 它的优势在于,它的两倍不超过 int 它的乘积不超过 long long 很有利于我们的运算,如果使用刚好不超过 long long 的数使用时很容易出现溢出并不方便。并且它恰好比较大。避免了做完 FFT 出现溢出。另外它可以取到的最大的 $N > 2e6$ 也很不错。例如现在如果我们要做 $2^k$, $k \leq 21$ 的 NFT。那么我们取 $w = g^{\frac{p-1}{2^k}}$ 即可。
HDU 1402
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
| const int N=132005; char sa[N>>1],sb[N>>1]; LL a[N],b[N]; LL pow_mod(LL x,LL n,LL p){ LL r=1; while(n){ if(n&1) r=r*x%p; n>>=1; x=x*x%p; } return r; } void change(LL *x,int len,int loglen){ for(int i=1;i<len;++i){ int t=i,k=0; for(int j=0;j<loglen;++j,t>>=1){ k=(k<<1)|(t&1); } if(k<i) swap(x[i],x[k]); } } const LL FM = 479<<21|1; void nft(LL *x,int len,int loglen,bool isInverse){ LL g = pow_mod(3,(FM-1)>>loglen,FM); if(isInverse){ g=inv(g,FM); LL invlen = pow_mod(len,FM-2,FM); for(int i=0;i<len;++i){ x[i]=x[i]*invlen%FM; } } change(x,len,loglen); for(int step=2;step<=len;step<<=1){ int half = step>>1; LL wn = pow_mod(g,len/step,FM); for(int i=0;i<len;i+=step){ LL w = 1; for(int j = i;j<i+half;++j){ LL t=(w*x[j+half])%FM; x[j+half]=(x[j]-t+FM)%FM; x[j]=(x[j]+t)%FM; w = w*wn%FM; } } } } int main(){ while(~scanf("%s%s",sa,sb)){ int alen=(int)strlen(sa); int blen=(int)strlen(sb); int len=1,loglen=0,tmp=alen+blen+3; while(len<tmp){ len<<=1;++loglen; } clr(a,0);clr(b,0); for(int i=0;i!=alen;++i) a[i]=sa[alen-i-1]-'0'; for(int i=0;i!=blen;++i) b[i]=sb[blen-i-1]-'0'; nft(a,len,loglen,0); nft(b,len,loglen,0); for(int i=0;i!=len;++i){ a[i] = a[i]*b[i]%FM; } nft(a,len,loglen,1); int cnt=0; while(cnt<len){ a[cnt+1]+=a[cnt]/10; a[cnt]%=10;++cnt; } cnt=alen+blen; while(cnt>1&&a[cnt-1]==0) --cnt; for(int i=cnt-1;i>=0;--i){ putchar((int)a[i]+'0'); } puts(""); } return 0; }
|