数论变换

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

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

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

变换公式

数论变换的变换式为

 

而数论变换的逆变换式为

 

注解:

(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.