數論轉換

Template:NoteTA

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

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

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

轉換公式

數論轉換的轉換式為

$F [k] = \sum_{n = 0}^{N - 1} f [n] α^{n k} (\mod M) k = 0, 1, 2, \dots, N - 1$

而數論轉換的反轉換式為

$f [n] = N^{- 1} \sum_{k = 0}^{N - 1} F [k] α^{- n k} (\mod M) n = 0, 1, 2, \dots, N - 1$

註解：

(1) $M$ 是一個質數。

(2) $(\mod M)$ 表示除以M的餘數

(3) $N$ 必須是 $M - 1$ 的因數。(當 $N \neq 1$ 時 $N$ 和 $M$ 互質)

(4) $N^{- 1}$ 是 $N$ 對模數 $M$ 之模反元素

(5) $α$ 為模M的N階單位根，即 $α^{N} = 1 (\mod M)$ 而且 $α^{k} \neq 1 (\mod M) k = 1, 2, \dots, N - 1$ 。若此時 $N = M - 1$ ，我們稱 $α$ 為模M的原根

舉一個例子：

一個 $M = 5$ 的 $N = 4$ 點數論轉換與反轉換如下，取 $α = 2$ ，注意此時 $N^{- 1} = 4$

正轉換 $(\begin{matrix} F [0] \\ F [1] \\ F [2] \\ F [3] \end{matrix}) = (\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & 2 & 4 & 3 \\ 1 & 4 & 1 & 4 \\ 1 & 3 & 4 & 2 \end{matrix}) (\begin{matrix} f [0] \\ f [1] \\ f [2] \\ f [3] \end{matrix}) (\mod M)$

反轉換 $(\begin{matrix} f [0] \\ f [1] \\ f [2] \\ f [3] \end{matrix}) = (\begin{matrix} 4 & 4 & 4 & 4 \\ 4 & 2 & 1 & 3 \\ 4 & 1 & 4 & 1 \\ 4 & 3 & 1 & 2 \end{matrix}) (\begin{matrix} F [0] \\ F [1] \\ F [2] \\ F [3] \end{matrix}) (\mod M)$

數論轉換的性質

(1) 正交性質

數論轉換矩陣的每個列是互相正交的，即 $\sum_{n = 0}^{N - 1} α^{n l} α^{- n k} = {\begin{matrix} N i f k = l \\ 0 i f k \neq l \end{matrix} (\mod M)$

(2) 對稱性

若 $f [n] = f [N - n]$ ，則 $f [n]$ 的數論轉換 $F [k]$ 也會有 $F [k] = F [N - k]$ 的特性。

若 $f [n] = - f [N - n]$ ，則 $f [n]$ 的數論轉換 $F [k]$ 也會有 $F [k] = - F [N - k]$ 的特性。

(3) 反數論轉換可由正數論轉換實現

$f [n] = N^{- 1} \sum_{k = 0}^{N - 1} F [k] α^{- n k} = N^{- 1} \sum_{- k = 0}^{N - 1} F [- k] α^{n k}$ ，即 $N^{- 1} F [- k]$ 的數論轉換。

步驟一：把 $F [k]$ 改寫成 $F [- k]$ ，若 $F [k] = F_{0}, F_{1},, F_{2} \dots, F_{N - 1}$ ，則 $F [- k] = F_{0}, F_{N - 1}, \dots, F_{2}, F_{1}$

步驟二：求 $F [- k]$ 的數論轉換。

步驟三：乘上 $N^{- 1}$ 。

(4) 線性性質

若 $f [n] \leftrightarrow F [k]$ ， $g [n] \leftrightarrow G [k]$ ，( $\leftrightarrow$ 表互為數論轉換對)則有 $a f [n] + b g [n] \leftrightarrow a F [k] + b g [k] (\mod M)$ 。

(5) 平移性質

若 $f [n] \leftrightarrow F [k]$ ，則 $f [n + l] \leftrightarrow α^{- l k} F [k] (\mod M)$ ，而且 $f [n] α^{n l} \leftrightarrow F [k + l]$ 。

(6) 圓周摺積性質

若 $f [n] \leftrightarrow F [k]$ ， $g [n] \leftrightarrow G [k]$ ，則 $f [n] \otimes g [n] \leftrightarrow F [k] G [k] (\mod M)$ ，而且 $f [n] g [n] \leftrightarrow \frac{1}{N} F [k] \otimes G [k] (\mod M)$ 。( $\otimes$ 代表圓形摺積。)

這個性質是數論轉換可以用來作為摺積的快速演算法的主要原因。

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

若 $f [n] \leftrightarrow F [k]$ ，則 $N \sum_{n = 0}^{N - 1} f [n] f [- n] = \sum_{k = 0}^{N - 1} F^{2} [k] (\mod M)$ ，而且 $N \sum_{n = 0}^{N - 1} f^{2} [n] = \sum_{k = 0}^{N - 1} F [k] F [- k] (\mod M)$

快速數論轉換

如果轉換點數N是一個2的次方，則可以使用類似基數-2快速傅立葉變換的蝴蝶結構來有效率的實現數論轉換。同樣的互質因子算法也可以應用在數論轉換上。

$\begin{matrix} F [k] & = & \sum_{n = 0}^{N - 1} f [n] α^{n k} \\ = & \sum_{r = 0}^{N / 2 - 1} f [2 r] α^{2 r k} + \sum_{r = 0}^{N / 2 - 1} f [2 r + 1] α^{(2 r + 1) k} \\ = & \sum_{r = 0}^{N / 2 - 1} f [2 r] (α^{2})^{r k} + α^{k} \sum_{r = 0}^{N / 2 - 1} f [2 r + 1] (α^{2})^{r k} \\ = & {\begin{matrix} G [k] + α^{k} H [k] 0 \leq k \leq \frac{N}{2} - 1 \\ G [k - \frac{N}{2}] - α^{k} H [k - \frac{N}{2}] \frac{N}{2} \leq k \leq N - 1 \end{matrix} \end{matrix}$

其中， $G [k] = \sum_{r = 0}^{N / 2 - 1} f [2 r] (α^{2})^{r k}$ ， $H [k] = \sum_{r = 0}^{N / 2 - 1} f [2 r + 1] (α^{2})^{r k}$ 。因此一個 $N$ 點數論轉換可以拆解成兩個 $\frac{N}{2}$ 點的數論轉換。

N, M, alpha的選擇

由於數論轉換可以擁有類似快速傅立葉變換的快速演算法，因此通常 $N$ 會選擇適合使用快速演算的值，比如說取為2的冪次，或是可以分解成許多小質數相乘的數。

在數論轉換中，需要大量地和 $α$ 的冪次做乘法，因此，如果可以取 $α$ 為2或2的冪次，則每一次的乘法在二進制中只會是一個移位的操作，可以省下大量的乘法運算。

因為要做模 $M$ 的運算，所以 $M$ 的二進位表示法中，1的個數越少越好，同時 $M$ 的值不能取太小，這是因為數論轉換後的值都小於 $M$ ，因此當真實的摺積的結果會大於 $M$ 時就會發生錯誤，所以必須謹慎選取 $M$ 的大小。

特殊的數論轉換

梅森質數數論轉換

梅森質數是指形如 $2^{p} - 1$ 的質數記做 $M_{p}$ ，其中 $p$ 是一個質數。

在梅森質數數論轉換中，取 $M = M_{p}$ ，可以證明 $α, N$ 可以如下選取：

(1) $α = 2, N = p, N^{- 1} = M_{p} - \frac{M_{p} - 1}{p}$

(2) $α = - 2, N = 2 p, N^{- 1} = M_{p} - \frac{M_{p} - 1}{2 p}$

在這兩種選取方式下，由於 $α$ 是2的冪次，所以只需移位運算，但 $N$ 不是2的冪次，所以基數-2的蝴蝶結構不適用於此例中，同時 $N$ 為質數或一個質數的兩倍，並不是一個可以拆解成很多質因數乘積的數，因此也無法由互質因子演算法得到太大的好處。梅森質數數論轉換通常用於較短序列的摺積運算

費馬數數論轉換

費馬數是指形如 $2^{2^{t}} + 1$ 的數記做 $F_{t}$ 。

在費馬數數論轉換中，取 $M = F_{t}$ ，可以證明在 $0 < t \leq 4$ 之下 $α, N$ 可以如下選取：

(1) $α = 2, N = 2^{t + 1}$

(2) $α = 3, N = F_{t} - 1$

在這兩種選取方式下， $N$ 是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為餘數，接上面的等式方程式化如下。

$M = q_{0} N + r_{1},$

$N = q_{1} r_{1} + r_{2},$

$r_{1} = q_{2} r_{2} + r_{3},$

$r_{2} = q_{3} r_{3} + r_{4}, . . .$

整理一下

$M - q_{0} N = r_{1},$

$N - q_{1} r_{1} = r_{2},$

$r_{1} - q_{2} r_{2} = r_{3},$

$r_{2} - q_{3} r_{3} = r_{4}, . . .$

最後的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來看，我們知道存在 $x M + y N = 1$ 。

gcd(N, M mod N)=1，則存在 $x_{1} N + y_{1} (M m o d N) = 1$

$x_{1} N + y_{1} (M m o d N) = x_{1} N + y_{1} (M - q_{1} N) = y_{1} M + (x_{1} - q_{1} y_{1}) N = 1$

這邊q就是相除的商數，比較M跟N的係數，這邊可以得到一個遞迴關係，

$(x_{1} - q_{1} y_{1}) = y$

$y_{1} = x$

可以藉由下一層的係數來推出上一層的係數。

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^{M - 1} m o d M = 1$

$N (N^{M - 2} m o d M) m o d M = 1$

顯然 $N^{M - 2} m o d M$ 就是N的反元素。

搭配快速mod可以容易的算出反元素，power是偶數時則可以用

(a^{2} m o d M)^{p o w e r / 2}

; power是基數時則可以用

(a N m o d M)^{p o w e r - 1}

，讓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

Root of Unity 是從1到M-1中選出一個適當的 alpha，其滿足"轉換公式"段落的(5)。

def RootOfUnity(N, M):
    def mod_exp(base, exp, mod):
        result = 1
        base %= mod
        while exp > 0:
            if exp % 2 == 1:  # If exp is odd
                result = (result * base) % mod
            base = (base * base) % mod
            exp //= 2
        return result

    for r in range(1, M):  # Check integers from 1 to M-1
        if mod_exp(r, N, M) == 1:
            return r  # Return the first root of unity
    return None

def NTT(x,N,M):
    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.

數論轉換

目录

轉換公式

數論轉換的性質

快速數論轉換

N, M, alpha的選擇

特殊的數論轉換

梅森質數數論轉換

費馬數數論轉換

實作議題

演算法

參考資料

导航菜单

數論轉換

轉換公式

數論轉換的性質

快速數論轉換

N, M, alpha的選擇

特殊的數論轉換

梅森質數數論轉換

費馬數數論轉換

實作議題

演算法

參考資料

导航菜单

搜索