跳至內容

數論轉換

本頁使用了標題或全文手工轉換
維基百科，自由的百科全書

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

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

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

轉換公式

數論轉換的轉換式為

而數論轉換的反轉換式為

註解：

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