數論轉換

數論轉換是一種計算摺積的快速演算法。計算摺積的快速演算法中最常用的一種是使用快速傅里葉變換,然而快速傅利葉變換具有一些實現上的缺點,舉例來說,資料向量必須乘上複數系數的矩陣加以處理,而且每個複數系數的實部和虛部是一個正弦餘弦函數,因此大部分的系數都是浮點數,也就是說,必須做複數而且是浮點數的運算,因此計算量會比較大,而且浮點數運算產生的誤差會比較大。

而在數論轉換中,資料向量需要乘上的矩陣是一個整數的矩陣,因此不需要作浮點數的乘法運算,更進一步,在模數為的情況下,只有種可能的加法與種可能的乘法,這可以使用記憶體把這些有限數目的加法和乘法存起來,利用查表法來計算,使得數論轉換完全不須任何乘法與加法運算,不過需要較大的記憶體與消耗較大的存取記憶體量。

雖然使用數論轉換可以降低計算摺積的複雜度,不過由於只進行整數的運算,所以只能用於對整數的訊號計算摺積,而利用快速傅利葉變換可以對任何複數的訊號計算摺積,這是降低運算複雜度所要付出的代價。

轉換公式

數論轉換的轉換式為

 

而數論轉換的反轉換式為

 

註解:

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