跳转到内容

数论变换

本页使用了标题或全文手工转换
维基百科，自由的百科全书

数论变换是一种计算卷积的快速算法。计算卷积的快速算法中最常用的一种是使用快速傅里叶变换，然而快速傅里叶变换具有一些实现上的缺点，举例来说，资料向量必须乘上复数系数的矩阵加以处理，而且每个复数系数的实部和虚部是一个正弦及余弦函数，因此大部分的系数都是浮点数，也就是说，必须做复数而且是浮点数的运算，因此计算量会比较大，而且浮点数运算产生的误差会比较大。

而在数论变换中，资料向量需要乘上的矩阵是一个整数的矩阵，因此不需要作浮点数的乘法运算，更进一步，在模数为的情况下，只有种可能的加法与种可能的乘法，这可以使用记忆体把这些有限数目的加法和乘法存起来，利用查表法来计算，使得数论变换完全不须任何乘法与加法运算，不过需要较大的记忆体与消耗较大的存取记忆体量。

虽然使用数论变换可以降低计算卷积的复杂度，不过由于只进行整数的运算，所以只能用于对整数的信号计算卷积，而利用快速傅里叶变换可以对任何复数的信号计算卷积，这是降低运算复杂度所要付出的代价。

变换公式

数论变换的变换式为

而数论变换的逆变换式为

注解：

(1) 是一个质数。

(2) 表示除以M的余数

(3) 必须是的因数。(当时和互质)

(4)是对模数之模反元素

(5)为模M的N阶单位根，即而且。若此时，我们称为模M的原根

举一个例子：

一个的点数论变换与逆变换如下，取，注意此时

正变换

逆变换

数论变换的性质

(1) 正交性质

数论变换矩阵的每个列是互相正交的，即

(2) 对称性

若，则的数论变换也会有的特性。

若，则的数论变换也会有的特性。

(3) 反数论变换可由正数论变换实现

，即的数论变换。

步骤一：把改写成，若，则

步骤二：求的数论变换。

步骤三：乘上。

(4) 线性性质

若，，(表互为数论变换对)则有。

(5) 平移性质

若，则，而且。

(6) 循环卷积性质

若，，则，而且。(代表圆形卷积。)

这个性质是数论变换可以用来作为卷积的快速算法的主要原因。

(7) 帕塞瓦尔定理(Parseval's Theorem)

若，则，而且

快速数论变换

如果变换点数N是一个2的次方，则可以使用类似基数-2快速傅里叶变换的蝴蝶结构来有效率的实现数论变换。同样的互质因子算法也可以应用在数论变换上。

其中，，。 因此一个点数论变换可以拆解成两个点的数论变换。

N, M, alpha的选择

由于数论变换可以拥有类似快速傅里叶变换的快速算法，因此通常会选择适合使用快速演算的值，比如说取为2的幂次，或是可以分解成许多小质数相乘的数。

在数论变换中，需要大量地和的幂次做乘法，因此，如果可以取为2或2的幂次，则每一次的乘法在二进制中只会是一个移位的操作，可以省下大量的乘法运算。

因为要做模的运算，所以的二进制表示法中，1的个数越少越好，同时的值不能取太小，这是因为数论变换后的值都小于，因此当真实的卷积的结果会大于时就会发生错误，所以必须谨慎选取的大小。

特殊的数论变换

梅森质数数论变换

梅森质数是指形如的质数记做，其中是一个质数。

在梅森质数数论变换中，取，可以证明可以如下选取：

(1)

(2)

在这两种选取方式下，由于是2的幂次，所以只需移位运算，但不是2的幂次，所以基数-2的蝴蝶结构不适用于此例中，同时为质数或一个质数的两倍，并不是一个可以拆解成很多质因数乘积的数，因此也无法由互质因子算法得到太大的好处。梅森质数数论变换通常用于较短序列的卷积运算

费马数数论变换

费马数是指形如的数记做。

在费马数数论变换中，取，可以证明在之下可以如下选取：

(1)

(2)

在这两种选取方式下，是2的幂次，所以基数-2的蝴蝶结构适用于此例中，而若是2的幂次，只需移位运算。费马数数论变换通常用于较长序列的卷积运算。

实作议题

由于数论变换需运用到余数下的反元素，这边提供一些不同的算法。

(一) Euclidean algorithm - iterative version

假设M为质数的mod，N为我们当前的元素且符合M-1的因数，借由Euclidean Algorithm知道必然存在x, y的整数使得

xM + yN = 1 - (1)

由上式左右mod M 可以得到 yN mod M= 1，显然y就是我们这边想求的反元素。

我们已知最大公因数(gcd)有如下性质。

gcd(M, N) = gcd(N, M mod N) = gcd(M mod N, N mod (M mod N)) = ... = 1

这边设定q为商数，r为余数，接上面的等式方程式化如下。

整理一下

最后的r 一定会变成1，所以我们只要不断的将r乘上-q带往下一个式子(像是r1*(-q1))，跟代往下下个式子(像是r3的左边式子要带入r1)即可求得最后我们想得到的 (1)，最后的N旁的系数就是反元素。

def modInv(x, M): #by Euclidean algorithm - iterative version
    t, new_t, r, new_r = 0, 1, M, x

    while new_r != 0:
        quotient = int(r / new_r)
        tmp_new_t = new_t
        tmp_t = t
        tmp_new_r = new_r
        tmp_r = r
        t, new_t = new_t, (t - quotient * new_t)
        r, new_r = new_r, (r % new_r)
    if r > 1:
        print("Inverse doesn't exist")
    if t < 0:
        t = t + M
    return t

(二) Euclidean algorithm - recursion version

根据gcd的性质gcd(M, N) = gcd(N, M mod N) = gcd(M mod N, N mod (M mod N)) = ... = 1

可以看出递回的关系，借此我们可以从这边得到N的系数，也就是反元素。

gcd(M, N) = 1来看，我们知道存在 。

gcd(N, M mod N)=1，则存在

这边q就是相除的商数，比较M跟N的系数，这边可以得到一个递回关系，

可以借由下一层的系数来推出上一层的系数。

def egcd(a, b): # y = x1 - quotient * y1  , x = y1 
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)

def modInv(n, m): #by Euclidean algorithm - recursion version 
    g, y, x = egcd(n, m)
    if g != 1:
        print("Inverse doesn't exist")
    else:
        return y % m

(三) Fermat little theorem

当M是质数时，我们知道任何一个数字N，

显然就是N的反元素。

搭配快速mod可以容易的算出反元素，power是偶数时则可以用; power是基数时则可以用，让power变成偶数。反复直到power变成1。

def modInv(a, m): #fermat little thm
      return modExponent(a, m - 2, m)
      
def modExponent(base, power, M): #quick mod O(log(power))
    result = 1
    power = int(power)
    base = base % M
    while power > 0:
        if power & 1:
            result = (result * base) % M
        base = (base * base) % M
        power = power >> 1
    return result

算法

以下程式预设。

import numpy as np

这边只要再从1到M-1中选出一个适当的alpha (Root Of Unity)，其满足"变换公式"段落的(5)即可。

def NTT(x,N,M):
    #TODO: RootOfUnity 
    alpha = RootOfUnity(N, M)
    NInv = modInv(N, M)
    A = np.ones((N,N))
    for i in range(1,N):
        for j in range(1,i+1):
          A[i][j] = A[i][j-1]*modExponent(alpha,i,M) % M
          A[j][i] = A[i][j]
    return np.matmul(A,x) % M, alpha

def invNTT(x,alpha,N,M):
    alphaInv = modInv(alpha, M)
    NInv = modInv(N, M)
    B = np.ones((N,N))
    for i in range(1,N):
        for j in range(1,i+1):
          B[i][j] = B[i][j-1]*modExponent(alphaInv,i,M) % M
          B[j][i] = B[i][j]
    B = B * NInv 
    return np.matmul(B,x) % M

参考资料

[1] R.C. Agarval and C.S. Burrus,"Number Theoretic Transforms to Implement Fast Digital Convolution," Proc. IEEE, vol.63, no.4, pp. 550-560, Apr. 1975.

[2] I. Reed and T.K. Truong,"The use of finite fileds to compute convolution," IEEE Trans. Info. Theory, vol.21 ,pp.208-213, Mar. 1975.

检索自“https://zh.wikipedia.org/w/index.php?title=數論轉換&oldid=74418034”