向量\(\mathbf{a}\)的DFT \(\mathbf{A}\)定义为：
\[A_j=\sum_{k=0}^{n-1}{a_k\omega^{jk}}\]
其中\(\omega\)是primitive root of unity of order \(n\)，复数时取\(exp(2i\pi/n)\)或exp(-2i/n)。
\[\sum_{k=0}^{n-1}{A_k\omega^{-jk}} = \sum_{k=0}^{n-1}{\sum_{l=0}^{n-1}{a_l\omega^{kl}}\omega^{-jk}} = \sum_{k=0}^{n-1}{\sum_{l=0}^{n-1}{a_l\omega^{(l-j)k}}} = na_j\]
因此
\[a_j=\frac{1}{n}\sum_{k=0}^{n-1}{A_k\omega^{jk}}\]
DFT可以用于复数以外的其他ring，常用于\(\mathbb{Z}/m\)
使用128 bits模数需要高效的uint64*uint64%uint64，其中模数是常数。
到AVX-512也没有提供把两个64 bits乘数的积放在一个128 bits寄存器的指令，GCC没有提供用乘法、移位等模拟的uint128除以uint64的常量除法。DIV指令性能很低。
当模数\(P<2^{63}\)时可以用64位mantissa浮点数计算uint64*uint64%uint64。
由等式\[(a\cdot b\% m) = a\cdot b - \lfloor\frac{\cdot b}{m}\rfloor\cdot m\]
两边模\(2^{64}\)，得 \[\begin{eqnarray*} (a\cdot b\% m)\% 2^{64} &=& (a\cdot b - \lfloor\frac{a\cdot b}{m}\rfloor\cdot m) \% 2^{64} \\ &=& ((a\cdot b)\%2^{64} - (\lfloor\frac{a\cdot b}{m}\rfloor\cdot m)\%2^{64}) \% 2^{64} \end{eqnarray*}\]即用uint64_t乘法计算\(a\cdot b\)的低64位，减去\(\lfloor\frac{\cdot b}{m}\rfloor\cdot m\)的低64位。其中\(\lfloor\frac{\cdot b}{m}\rfloor<m\)，可以用64位mantissa浮点数(Intel x87 80-bit double-extended precision)计算。
若\(m<2^{63}\)时，浮点数转uint64_t时取round时是否向上取整了，差的绝对值将会大于等于\(2^{63}\)，此时再加上\(P\)。若\(m\geq 2^{63}\)，则没有简单的办法判断是否向上取整了。
|
|
存储\(P\)的倒数\(Q=\frac{1}{P}\)，用(long double)a*b*Q代替(long double)a*b/P能快些。此时\(Q\)会引入额外的误差，Matters Computational说适用于\(m<2^{62}\)，原因不明。
\[(a\ast b)_j = \sum_{k=0}^{n-1}{a_kb_{(j-k)\;mod\;n}}\]
性质：\(\frac{1}{n}\sum_{k=0}^{n-1}{\omega^{jk}} = [j\;mod\;n=0]\)
\[\begin{eqnarray*} (a\ast b)_j &=& \sum_{p=0}^{n-1}{\sum_{q=0}^{n-1}{[(p+q)\;mod\;n=j]a_pb_q}}\\ &=& \sum_{p=0}^{n-1}{\sum_{q=0}^{n-1}{[(p+q-j)\;mod\;n=0]a_pb_q}}\\ &=& \sum_{p=0}^{n-1}{\sum_{q=0}^{n-1}{\frac{1}{n}\sum_{k=0}^{n-1}{\omega^{pk}\omega^{qk}\omega^{-jk}}a_pb_q}}\\ &=& \frac{1}{n}\sum_{k=0}^{n-1}{(\sum_{p=0}^{n-1}{a_p\omega^{kp}})(\sum_{q=0}^{n-1}{b_q\omega^{kq}})\omega^{-jk}} \\ &=& \frac{1}{n}\sum_{k=0}^{n-1}{(A_kB_k)\omega^{-jk}} \\ \end{eqnarray*}\]如果要计算linear convolution，比如多项式乘法，可以把\(\mathbf{a}\)长度补到\(2n-1\)，求cyclic convolution。
下面讨论用于计算整系数向量卷积的discrete Fourier transform。
使用complex<double>计算convolution，需要保证结果每一项的实数部分在\([-2^{53}-1, 2^{53}-1]\)(Number.MIN_SAFE_VALUE、Number.MAX_SAFE_INTEGER)范围内，\(2^{53}-1\)是double能精确表示的最大整数。采取round half to even规则，\(2^{53},2^{53}+1\)均表示为\(2^{53}\)，无法区分。
设每项系数的绝对值小于等于\(v\)，那么convolution结果每一项绝对值小于等于\(nv^2\)，若\(nv^2\leq 2^{53}-1\)则可放心使用complex<double>。
complex<double>还要受到浮点运算误差影响。根据Roundoff Error Analysis of the Fast Fourier Transform，没仔细看，relative error均值为 log2(n)浮点运算精度变换前系数最大值。
对于模\(P\)的number theoretic transform，\(v\leq P-1\)，若\(nv^2\leq P\)则可放心使用。
998244353, 897581057, 880803841小于\(2^{30}\)，两倍不被超过INT_MAX，且可表示为\(k*n+1\)，其中\(n\)为2的幂，适合用作number theoretic transform的模。
设\(\mathbf{a},\mathbf{b}\)系数取自\([0,v]\)的uniform distribution，则\(\mathbf{a}\ast\mathbf{b}\)系数均值为\(nv^2/4\)，方差为\(nv^4/9\)。若把系数平移至\([-v/2, v/2]\)，则\(\mathbf{a}\ast\mathbf{b}\)系数均值\(0\)，方差为\(nv^4/144\)。若\(\mathbf{a},\mathbf{b}\)其中之一independent and identically distributed，则方差会很小。Chebyshev’s inequality等可用于估计系数绝对值落在方差若干倍内的概率，上界\(nv^2\)可以减小。即使\(\mathbf{a},\mathbf{b}\)不是independent and identically distributed，也可以用\(\mathbf{a}\ast\mathbf{b}=\mathbf{a}\ast(\mathbf{b+c})-\mathbf{a}\ast\mathbf{c}\)来计算，\(\mathbf{c}\)是independent and identically distributed uniform \([-v/2,v/2]\)。
取\(M\)为接近\(\sqrt{v}\)的整数，分解\(\mathbf{a}=\mathbf{a0}+M\mathbf{a1}\)、\(\mathbf{a}=\mathbf{a0}+M\mathbf{a1}\)，则：
\[\mathbf{a}\ast\mathbf{b} = \mathbf{a0}\ast\mathbf{b0}+M(\mathbf{a0}\ast\mathbf{b1}+\mathbf{a1}\ast\mathbf{b0})+M^2(\mathbf{a1}\ast\mathbf{b1})\]
适当选择\(M\)可以使\(\mathbf{a0},\mathbf{a1}\)的系数小于等于\(\lfloor\sqrt{v}\rfloor\)，convolution结果系数最大值为\(n(\lfloor\sqrt{v}\rfloor)^2\approx nv\)，比原来的\(nv^2\)小。
求出\(DFT(\mathbf{a0}), DFT(\mathbf{a1}), DFT(\mathbf{b0}), DFT(\mathbf{b1})\)后，计算等式右边四个convolution，带权相加即得到原convolution。
需要4次长为\(2n\)的DFT、1次长为\(2n\)inverse DFT。
可以使用Toom-2 (Karatsuba)计算\(\mathbf{a0}\ast\mathbf{b0}, \mathbf{a1}\ast\mathbf{b1}, (\mathbf{a0}+\mathbf{a1})\ast(\mathbf{b0}+\mathbf{b1})\)，减少为3次DFT、1次inverse DFT。\(\mathbf{a0}\ast\mathbf{b1}+\mathbf{a1}\ast\mathbf{b0} = (\mathbf{a0}+\mathbf{a1})\ast(\mathbf{b0}+\mathbf{b1}) - \mathbf{a0}\ast\mathbf{b0} - \mathbf{a1}\ast\mathbf{b1}\)。
取\(S\)与\(\sqrt(P)\)接近且\(M=P-S*S%P\)尽可能小。
分解\(\mathbf{a}=\mathbf{a0}+S\mathbf{a1}\)、\(\mathbf{b}=\mathbf{b0}+S\mathbf{b1}\)。
\[\begin{eqnarray*} IDFT(DFT(a0+i\sqrt{S}a1)\cdot DFT(b0+i\sqrt{S}a1)) &=& (a0+i\sqrt{S}a1)\ast (b0+i\sqrt{S}b1) \\ &=& (\mathbf{a0}\ast \mathbf{b0})-S(\mathbf{a1}\ast\mathbf{b1})+i\sqrt{S}((\mathbf{a0}\ast\mathbf{b1})+(\mathbf{a1}\ast\mathbf{b0})) \\ \end{eqnarray*}\]右边同余于\((\mathbf{a0}\ast \mathbf{b0})+S^2(\mathbf{a1}\ast\mathbf{b1})+i\sqrt{S}(\mathbf{a0}\ast\mathbf{b1}+\mathbf{a1}\ast\mathbf{b0})\)，提取虚部与实部即可得到：
\[(\mathbf{a0}\ast \mathbf{b0})+S(\mathbf{a0}\ast\mathbf{b1}+\mathbf{a1}\ast\mathbf{b0})+S^2(\mathbf{a1}\ast\mathbf{b1})\]
需要2次长为\(2n\)DFT、1次长为\(2n\)inverse DFT。
\(\mathbf{a}\)的共轭的DFT可由\(\mathbf{a}\)的DFT求出：
\[\begin{eqnarray*} DFT(conj(\mathbf{a}))_j &=& \sum_{k=0}^{n-1}{conj(a_k)\omega^{jk}} \\ &=& \sum_{k=0}^{n-1}{conj(a_k\omega^{-jk})} \\ &=& conj(\sum_{k=0}^{n-1}{a_k\omega^{-jk}}) \\ &=& conj(DFT(\mathbf{a})_{-j\;mod\;n}) \\ &=& conj(rev(DFT(\mathbf{a}))_j) \\ \end{eqnarray*}\] \[\begin{eqnarray*} DFT(re(\mathbf{a})) &=& (DFT(\mathbf{a})+DFT(conj(\mathbf{a})))/2 = (DFT(\mathbf{a}) + conj(rev(DFT(\mathbf{a}))))/2 \\ DFT(im(\mathbf{a})) &=& (DFT(\mathbf{a})-DFT(conj(\mathbf{a})))/(2i) = (DFT(\mathbf{a}) - conj(rev(DFT(\mathbf{a}))))/(2i) \\ \end{eqnarray*}\]分解\(\mathbf{a}=\mathbf{a0}+M\mathbf{a1}\)、\(\mathbf{b}=\mathbf{b0}+M\mathbf{b1}\)后，根据上面的公式，用\(DFT(\mathbf{a0}+i\mathbf{a1})\)计算\(DFT(\mathbf{a0})\)和\(DFT(\mathbf{a1})\)，同法计算\(DFT(\mathbf{b0})\)和\(DFT(\mathbf{b1})\)。然后用\(IDFT(DFT(\mathbf{a0})\cdot DFT(\mathbf{b0}) + i DFT(\mathbf{a0})\cdot DFT(\mathbf{b1}))\)计算出\(\mathbf{a0}\ast\mathbf{b0}\)与\(\mathbf{a0}\ast\mathbf{b1}\)，同法计算出\(\mathbf{a1}\ast\mathbf{b0}\)与\(\mathbf{a1}\ast\mathbf{b1}\)。
需要2次长为\(2n\)DFT、2次长为\(2n\)inverse DFT。
该优化可以和其他方式叠加。
把\(\mathbf{a}\)看作多项式\(A(x)=\sum_{k=0}^{n-1}{a_kx^k}\)，同样地，\(\mathbf{b}\)看作多项式\(B(x)\)。
奇次项\(A0(x)=\sum_{0\leq k<n, k\;\text{is even}}{a_kx^k}\), 偶次项\(A1(x)=\sum_{0\leq k<n, k\;\text{is odd}}{a_kx^k}\)，同样地，定义\(B0(x)\)与\(B1(x)\)。\(A0(x)\)的系数为\(\mathbf{a0}，\)A1(x)\(的系数为\)\(，令其长为\)n$，高位用\(0\)填充。
\[\begin{eqnarray*} A(x)B(x) &=& (A0(x^2)+x A1(x^2))(B0(x^2)+x B1(x^2)) \\ &=& A0(x^2)B0(x^2)+x^2A1(x^2)B1(x^2) + x(A0(x^2)B1(x^2)+A1(x^2)B0(x^2)) \\ \end{eqnarray*}\]用正交计算两个实系数向量DFT的方式，用2次长度为\(n\)(之前都是 \(2n\))的DFT计算\(DFT(\mathbf{a0}), DFT(\mathbf{a1})\), \(DFT(\mathbf{b0}), DFT(\mathbf{b1})\)。\(\mathbf{a1}\)循环右移1位的DFT的第\(j\)项等于\(DFT(\mathbf{a1})_j\omega^j\)，因此根据\(A1(x^2)\)的DFT的系数可以得到\(x^2A1(x^2)\)的DFT的系数。
构造长为\(n\)的向量\(\mathbf{c}\)： \[c_j = (a0_j b0_j + a1_j b1_j \omega^j) + i(a0_j b1_j + a1_j b0_j)\]
\(DFT(\mathbf{c})\)的实部为结果的偶次项系数，虚部为结果的奇次项系数。
需要2次长为\(n\)的DFT、1次长为\(n\)的inverse DFT。
取\(t\)个可用于number theoretic transform的质数\(P_0,P_1,\ldots,P_{t-1}\)，使\(M = \prod_{0\leq j<t}{P_j}\geq nv^2\)，计算\(t\)个NTT，之后用Chinese remainder theorem合并。
求\(x \equiv v_j\quad(\text{mod}\;M/P_j),\quad 0\leq j<t\)有至少两种算法。
Gauss之前也有很多人提出。
对于每个\(P_j\)用Blankinship’s algorithm计算\(P_j p_j \equiv 1\quad(\text{mod}\;M/P_j)\)。
\[x = (\sum_{j=0}^{t-1}{v_jP_jp_j}) \% M\]
注意到\(M\geq nv^2\)可能超出机器single-precision表示范围。
\(x_j\)满足前\(j+1\)个方程，\(x=x_{t-1}\)满足所有方程。
稍加变形可用于求\(x\% P\)。
\[\begin{eqnarray*} y_0 &=& v_0 \% P_0 & x_0 &=& y_0 \% P \\ y_1 &=& (v_1-(y_0)\% P_1)c_1 \% P_1 & x_1 &=& (x_0 + y_1P_0) \% P \\ y_2 &=& (v_2-(y_0+y_1P_0)\% P_2)c_2 \% P_2 & x_2 &=& (x_1 + y_2P_0P_1) \% P \\ y_3 &=& (v_3-(y_0+y_1P_0+y_2P_0P_1)\% P_3)c_3 \% P_3 & x_3 &=& (x_2 + y_3P_0P_1P_2) \% P \\ \ldots \\ \end{eqnarray*}\]原来的每个\(x_j\)需要计算\(x_j\%P\)和\(x_j\%P_{j+1}\)两份，时间复杂度上升为\(O(t^2)\)。