技术控

    今日:133| 主题:49390
收藏本版 (1)
最新软件应用技术尽在掌握

[其他] 整系数Discrete Fourier transform技巧

[复制链接]
房妹 发表于 2016-10-4 10:29:07
48 2

立即注册CoLaBug.com会员,免费获得投稿人的专业资料,享用更多功能,玩转个人品牌!

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
Discrete Fourier transform

  向量    \(\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)。  
  Inverse transform

  \[\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}}\]
  Number theoretic transform

  DFT可以用于复数以外的其他ring,常用于    \(\mathbb{Z}/m\)  
  使用128 bits模数需要高效的    uint64*uint64%uint64,其中模数是常数。  
  硬件除法指令

  到AVX-512也没有提供把两个64 bits乘数的积放在一个128 bits寄存器的指令,GCC没有提供用乘法、移位等模拟的uint128除以uint64的常量除法。    DIV指令性能很低。  
  64位mantissa浮点数

  当模数    \(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}\),则没有简单的办法判断是否向上取整了。  
  1. uint64_tmul_mod(uint64_ta,uint64_tb,longP)
  2. {
  3. uint64_tx = a*b, r = x - P*uint64_t((longdouble)a*b/P+0.5);
  4. returnint64_t(r) <0? r + P : r;
  5. }
复制代码
存储    \(P\)的倒数    \(Q=\frac{1}{P}\),用    (long double)a*b*Q代替    (long double)a*b/P能快些。此时    \(Q\)会引入额外的误差,Matters Computational说适用于    \(m<2^{62}\),原因不明。  
  Cyclic convolution

  \[(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]\)。  
  方案0:sqrt decomposition(FFT, NTT)

  取    \(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}\)。  
  优化0:imaginary unit用于进位(FFT)

  取    \(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。  
  优化1:正交计算两个实系数向量DFT(FFT)

      \(\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。  
  奇偶项优化(FFT)

  该优化可以和其他方式叠加。
  把    \(\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。  
  方案1:Chinese remainder theorem(NTT)

  取    \(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’s algorithm)

  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表示范围。  
  Garner’s algorithm

定义  \[\begin{eqnarray*} c_1 &=& inv(P_0, P_1) \\ c_2 &=& inv(P_0P_1, P_2) \\ c_3 &=& inv(P_0P_1P_2, P_3) \\ \ldots \\ \end{eqnarray*}\]  \[\begin{eqnarray*} y_0 &=& v_0 \% P_0 & x_0 &=& y_0 \\ y_1 &=& (v_1-x_0)c_1 \% P_1 & x_1 &=& x_0 + y_1P_0 \\ y_2 &=& (v_2-x_1)c_2 \% P_2 & x_2 &=& x_1 + y_2P_0P_1 \\ y_3 &=& (v_3-x_2)c_3 \% P_3 & x_3 &=& x_2 + y_3P_0P_1P_2 \\ \ldots \\ \end{eqnarray*}\]
      \(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)\)。  
  References

  感谢ftiasch老师教导。
  
       
  • fotile96, https://async.icpc-camp.org/d/408-fft   
  • Jörg Arndt, Matters Computational   
  • 毛啸,再探快速傅里叶变换,2016年信息学奥林匹克中国国家队候选队员论文集  
友荐云推荐




上一篇:国庆节,我们用代码来画个国旗~
下一篇:Python科学计算与数据分析库/包大全
酷辣虫提示酷辣虫禁止发表任何与中华人民共和国法律有抵触的内容!所有内容由用户发布,并不代表酷辣虫的观点,酷辣虫无法对用户发布内容真实性提供任何的保证,请自行验证并承担风险与后果。如您有版权、违规等问题,请通过"联系我们"或"违规举报"告知我们处理。

那年的那年 发表于 2016-10-6 05:59:41
听说顶贴有福利,真的假的?
回复 支持 反对

使用道具 举报

酷u 发表于 2016-10-25 14:45:09
房妹是男的还是女的?
回复 支持 反对

使用道具 举报

*滑动验证:
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

我要投稿

推荐阅读

扫码访问 @iTTTTT瑞翔 的微博
回页顶回复上一篇下一篇回列表手机版
手机版/CoLaBug.com ( 粤ICP备05003221号 | 文网文[2010]257号 )|网站地图 酷辣虫

© 2001-2016 Comsenz Inc. Design: Dean. DiscuzFans.

返回顶部 返回列表