快速数论变换(NTT)是快速傅立叶算法(FFT)的进一步算法,为的是解决FFT算法中的大量浮点数计算造成的误差。如果您还不知道什么是快速傅立叶算法,建议先看看 快速傅立叶算法之大数相乘这篇文章,因为NTT里面的很多东西都和FFT相似,如果不了解清楚,可能无法理解本文内容。

相关知识

质数大家应该都知道吧,就是只能被1和自己整除的数字,比如3,5,7,11等。

现在讲一个概念,叫做欧拉函数:φ(n)。这玩意是干嘛的呢?它是用来计算在小于n的数字中,和它互质数字的个数。比如数字5,在5以内和它互质的数字为1,2,3,4,一共有4个,因此φ(5)=4;数字7,在7以内和它互质的为1,2,3,4,5,6,共6个,因此φ(7)=6。我们很容易发现,如果一个数字p是质数,那么φ(p)=p-1。(当然了,如果p不质数,结果肯定不是p-1哈。)

另一方面,我们伟大的数学家欧拉,发现了一个有趣的定理,如果数字a和m互质,那么(a^m)%m=1。例如Math.pow(7,3)%3=1,Math.pow(7,6)%6=1

那么这两个知识点有什么关联吗?当然有,为的就是推出一个更重要的知识点:原根。原根的条件就比较苛刻了,如果(a^p)%m(1<=p<=φ(m))的结果两两不相同的话,那么a就是m的原根。例如3就是7的原根,首先φ(7)=6,然后3^1%7=3,3^2%7=2,3^3%7=6,3^4%7=4,3^5%7=5,3^6%7=1。我们可以发现,3的1,2,3,4,5,6次方各不相同,

原根的幂取模是一个循环的数列,如上面的3对于7来说是原根,会产生3,2,6,4,5,1,3,2,6,4,5,1。。这样的循环数列,可以看到不变的其实就是3,2,6,4,5,1这几个数字。同时,方程式y=a0+a1x+a2x^2.....anx^n,应该可以利用原根的这种性质来做些什么。

当然数学家确实是发现了这其中的奥妙,我们可以找到形如a*2^b+1的原根,如17=1*2^4+1,它的原根是3,它的循环数列为3,9,10,13,5,15,11,16,14,8,7,4,12,2,6,1,共16个数字。我们挑取出其中第2^0(1),2^1(2),2^2(4),2^3(8),2^4(16)个数字,分别是3,9,13,16,1。这几个数字有什么特殊呢?我们可以发现3*3%17=9,9*9%17=13,13*13%17=16,16*16%17=1,也就是前面的数字自相乘得到的是后面的数字,这个和我们的FFT中的单位复根比较类似。数学家已经证明了NTT可以采取和FFT类似的方法来计算,但是具体的细节可能有所不同。

NTT和逆NTT

还记得我们计算FFT的时候有个W(n,k)吗?在NTT计算的时候,我们用的是G(n,k)。和FFT不同,G(n,k)是事先计算好的,而且不是唯一的。就拿上面的3和17来说,3是17的原根,那么,G(4,0)=1,G(4,1)=16,G(4,2)=13,G(4,3)=9,G(4,4)=3,类似的,我们利用5,97产生的G(n,k)来简单演示一下如何跑NTT和逆NTT。

background Layer 1 15476 0x^7+0x^6+0x^5+x^4+5x^3+4x^2+7x+6 0 0 0 1 5 4 7 6 0 0 5 7 0 1 4 6 0 5 NTT (按照奇偶下标拆分) 自上而下 我们采用的是5和97,5是97的原根,得到 G(6,0)=1,G(6,1)=96,G(6,2)=22,G(6,3)=64,G(6,4)=8,G(6,5)=28,G(6,6)=5 自下而上 0 7 0 5 0 7 0 4 1 6 0 4 1 6 0 5 0 7 0 4 1 6 G(6,0) (6+1*G(6,0)%17)%17 (6-1*G(6,0)%17)+17%17 G(6,0) (4+0*G(6,0)%17)%17 (4-0*G(6,0)%17)+17%17 G(6,0) (7+0*G(6,0)%17)%17 (7-0*G(6,0)%17)+17%17 G(6,0) (5+0*G(6,0)%17)%17 (5-0*G(6,0)%17)+17%17 7 5 4 4 7 7 5 5 (7+4*G(6,0)%17)%17 G(6,0) 11 3 (5+4*G(6,2)%17)%17 G(6,2) 93 14 (7+5*G(6,0)%17)%17 G(6,0) 12 2 (7+5*G(6,2)%17)%17 G(6,2) 20 91 G(6,0) (11+12*G(6,0)%17)%17 23 96 G(6,3) (93+20*G(6,3)%17)%17 15 74 G(6,2) G(6,2)*G(6,3)%97 47 5 56 23 (7-4*G(6,0)%17)+17%1 (5-4*G(6,2)%17)%17 (7-5*G(6,0)%17)+17%1 (7-5*G(6,2)%17)%17 (3+2*G(6,2)%17)%17 (14+91*G(6,2.5)%17)%17 (11-12*G(6,0)%17)+17)%16 (93-20*G(6,3)%17)%17 (3-2*G(6,2)%17)+17)%16 (14-91*G(6,2.5)%17)%17 第一组 第二组 第三组

和FFT一样,我们首先采用分治,将形参划分,然后同样是蝴蝶操作,只不过这里注意的是,FFT乘的是W(n,k),这里乘以的是G(n,k),这里的G(n,k)的规律比较复杂,每组第一个相乘的一定是1,组内第二个是G(n,k),第三个是G(n,k-1),第四个是G(n,(2k-1)/2),第五个是G(n,k-2).....FFT中的原根在计算的时候固定的,在NTT里,类似的计算单元为G(n,k),之后的每个每个都是乘以它然后取模。

不过采用逆NTT,这里会有一个概念,叫做逆元,这个需要对最后的结果进行相乘取模操作,逆元=a^(φ(n)-1)%m。

background Layer 1 y0=23 y1=15 y2=47 y3=5 y4=96 y5=74 y6=56 y7=23 23 56 74 96 5 47 15 23 23 74 5 15 56 96 47 23 23 5 逆向NTT (按照奇偶下标拆分) 自上而下 我们采用的是5和97,5是97的原根,得到 G(6,0)=1,G(6,1)=96,G(6,2)=22,G(6,3)=64,G(6,4)=8,G(6,5)=28,G(6,6)=5 自下而上 74 15 23 5 74 15 56 47 96 23 56 47 96 23 23 5 74 15 56 47 96 23 G(6,0) G(6,0) G(6,0) G(6,0) 22 24 6 88 89 38 28 79 G(6,0) 28 16 G(6,2) 20 28 G(6,0) 20 61 G(6,2) 30 46 G(6,0) 48 8 G(6,3) 0 40 G(6,2) G(6,2)*G(6,3)%97 0 0 32 56 第一组 第二组 第三组 逆元_G=8^95%97=85 48 8 0 40 0 0 32 56 交换位置 48 8 56 0 32 40 0 0 6 1 7 0 4 5 0 0 *85%97 结果15476

逆NTT最后的时候,需要交换位置,交换位置的时候除了第一个不交换外,其余的都和n/2后的内容交换。

上面的逆NTT已经成功的将值转成了形参,说明NTT的形转值,值转形这两条路是互通的,那么我么就可以用来做类似FFT的大数相乘了。

NTT之大数相乘

现在我们就来走一遍大数相乘,这里我们依旧采用5和97来做取模操作,当然实际上计算不会用它们,原因后面会说明。

background Layer 1 13*26 1.找到两者中位数较长的*2,找到最接近这个数且不小于它的2^n。这里为4。 2.对13做NTT 0x^3+0x^2+1x+3 0 0 1 3 0 1 0 3 0 1 0 3 自上而下 0 1 0 3 (3+(0*G(6,0)%17))%17 G(6,0) 我们采用的是5和97,5是97的原根,得到 G(6,0)=1,G(6,1)=96,G(6,2)=22,G(6,3)=64,G(6,4)=8,G(6,5)=28,G(6,6)=5 (3-(0*G(6,0)%17)+17)%17 3 3 (1+(0*G(6,0)%17))%17 G(6,0) (1-(0*G(6,0)%17)+17)%17 1 1 (3+(1*G(6,0)%17))%17 (3-(1*G(6,0)%17)+17)%17 G(6,0) (3+(1*G(6,2)%17))%17 (3-(1*G(6,2)%17)+17)%17 G(6,2) 4 25 2 78 自下而上 3.对26做NTT 0x^3+0x^2+2x+6 0 0 2 6 0 2 0 6 0 2 0 6 自上而下 0 2 0 6 (3+(0*G(6,0)%17))%17 G(6,0) (3-(0*G(6,0)%17)+17)%17 6 6 (2+(0*G(6,0)%17))%17 G(6,0) (2-(0*G(6,0)%17)+17)%17 2 2 (6+(2*G(6,0)%17))%17 (6-(2*G(6,0)%17)+17)%17 G(6,0) (6+(2*G(6,2)%17))%17 (6-(2*G(6,2)%17)+17)%17 G(6,2) 8 50 4 59 自下而上 4.乘法 y0=4*8%97=32 y1=25*50%97=86 y2=2*4%97=8 y3=78*59%97=43 5.逆NTT 43 8 86 32 43 86 8 32 43 86 8 32 自上而下 43 86 8 32 (32+(8*G(6,0)%17))%17 G(6,0) (32-(8*G(6,0)%17)+17)%17 40 24 (86+(43*G(6,0)%17))%17 G(6,0) (86-(43*G(6,0)%17)+17)%17 32 43 (40+(32*G(6,0)%17))%17 (40-(32*G(6,0)%17)+17)%17 G(6,0) (24+(43*G(6,2)%17))%17 (24-(43*G(6,2)%17)+17)%17 G(6,2) 72 0 8 48 自下而上 逆元:4^95%97=73 48 8 0 72 元素交换 *13%17 0 2 12 18 0 8 48 72 0 3 3 8 进位

最后

因为我们采用的是a*2^b+1的形式产生的质数,因此,我们可以处理2^b内的数字,如17,只可以处理2^4=8位数之内的乘法(实际上更小),因此,想让范围更大一点,我们需要吧质数设置的更大一些,这里我设置的是23068673,原根是3,当然这不是唯一的答案。

最后给出测试demo:

最后再唠叨一句,稍微讲一下程序中求a^b%m的方法。

我们都知道a^b是一个很大的数字,例如3^96,如果b持续增大的话,js可能就无法准确计算出a^b,这样就计算不出a^b%m。

程序中有一个计算a^b%m的方法,代码如下:

原理其实就是逐步求模,例如3%17=3,那么3*3%17=9(3^2%17),9*9%17=13(3^4%17),13*13%17=16(3^8%17)。。。这种计算方法和二进制刚好不谋而合,例如3^3%17=9*3%17=10,其中3为3%17,9为3^2%17,那和二进制有什么关系呢?3^3中的幂为3,3的二进制为11,3%17相当于2^0,3^2%17相当于2^1。所以3^3%17就3^(2^0)%17和3^(2^1)%17之积%17。

background Layer 1 3^17%17 17的二进制为10001,包含2^0和2^4 3%17=3 3*3%17=9 9*9%17=13 13*13%17=16 16*16%17=1 相当于2^0 相当于2^1 相当于2^2 相当于2^3 相当于2^4 结果3*1%17=3 10001(包含) 1000 100 10 1(包含)

何如判定包含呢?将幂的二进制不断的右移,如果最后一位有1,那么就表示包含,这个就是代码中if(b&1){}的意思。

DEMO下载

点击下载 [0积分]一共下载0

其他文章

0
我要评论

评论

返回
×

我要评论

回复:

昵称:(昵称不超过20个字)

图片:

提交
还可以输入500个字