數論轉換

數論轉換是一種計算摺積的快速演演算法。計算摺積的快速演演算法中最常用的一種是使用快速傅裡葉變換，然而快速傅立葉變換具有一些實現上的缺點，舉例來說，資料向量必須乘上複數係數的矩陣加以處理，而且每個複數係數的實部和虛部是一個正弦及餘弦函式，因此大部分的係數都是浮點數，也就是說，必須做複數而且是浮點數的運算，因此計算量會比較大，而且浮點數運算產生的誤差會比較大。而在數論轉換中，資料向量需要乘上的矩陣是一個整數的矩陣，因此不需要作浮點數的乘法運算，更進一步，在模數為formula_1的情況下，只有formula_2種可能的加法與formula_2種可能的乘法，這可以使用記憶體把這些有限數目的加法和乘法存起來，利用查表法來計算，使得數論轉換完全不須任何乘法與加法運算，不過需要較大的記憶體與消耗較大的存取記憶體量。雖然使用數論轉換可以降低計算摺積的複雜度，不過由於只進行整數的運算，所以只能用於對整數的訊號計算摺積，而利用快速傅立葉變換可以對任何複數的訊號計算摺積，這是降低運算複雜度所要付出的代價。轉換公式.數論轉換的轉換式為formula_4而數論轉換的反轉換式為formula_5註解：(1) formula_1是一個質數。(2) formula_7 表示除以M的餘數(4)formula_13是formula_8對模數formula_1之模反元素(5)formula_16為模M的N階單位根，即formula_17而且formula_18。若此時formula_19，我們稱formula_16為模M的原根舉一個例子：一個formula_21的formula_22點數論轉換與反轉換如下，取formula_23，注意此時formula_24正轉換formula_25反轉換formula_26數論轉換的性質.(1) 正交性質數論轉換矩陣的每個列是互相正交的，即formula_27(2) 對稱性若formula_28，則formula_29的數論轉換formula_30也會有formula_31的特性。若formula_32，則formula_29的數論轉換formula_30也會有formula_35的特性。(3) 反數論轉換可由正數論轉換實現formula_36，即formula_37的數論轉換。步驟一：把formula_30改寫成formula_39，若formula_40，則formula_41步驟二：求formula_39的數論轉換。步驟三：乘上formula_13。(4) 線性性質若formula_44，formula_45，(formula_46表互為數論轉換對)則有formula_47。(5) 平移性質若formula_44，則formula_49，而且formula_50。(6) 圓周摺積性質若formula_44，formula_45，則formula_53，而且formula_54。(formula_55代表圓形摺積。)這個性質是數論轉換可以用來作為摺積的快速演演算法的主要原因。若formula_44，則formula_57，而且formula_58快速數論轉換.如果轉換點數N是一個2的次方，則可以使用類似基數-2快速傅立葉變換的蝴蝶結構來有效率的實現數論轉換。同樣的互質因子演算法也可以應用在數論轉換上。formula_59其中，formula_60，formula_61。因此一個formula_8點數論轉換可以拆解成兩個formula_63點的數論轉換。N, M, alpha的選擇.由於數論轉換可以擁有類似快速傅立葉變換的快速演演算法，因此通常formula_8會選擇適合使用快速演算的值，比如說取為2的冪次，或是可以分解成許多小質數相乘的數。在數論轉換中，需要大量地和formula_16的冪次做乘法，因此，如果可以取formula_16為2或2的冪次，則每一次的乘法在二進制中只會是一個移位的操作，可以省下大量的乘法運算。因為要做模formula_1的運算，所以formula_1的二進位表示法中，1的個數越少越好，同時formula_1的值不能取太小，這是因為數論轉換後的值都小於formula_1，因此當真實的摺積的結果會大於formula_1時就會發生錯誤，所以必須謹慎選取formula_1的大小。特殊的數論轉換.梅森質數數論轉換.梅森質數是指形如formula_73的質數記做formula_74，其中formula_75是一個質數。在梅森質數數論轉換中，取formula_76，可以證明formula_77可以如下選取：(1) formula_78(2) formula_79在這兩種選取方式下，由於formula_16是2的冪次，所以只需移位運算，但formula_8不是2的冪次，所以基數-2的蝴蝶結構不適用於此例中，同時formula_8為質數或一個質數的兩倍，並不是一個可以拆解成很多質因數乘積的數，因此也無法由互質因子演演算法得到太大的好處。梅森質數數論轉換通常用於較短序列的摺積運算費馬數數論轉換.費馬數是指形如formula_83的數記做formula_84。在費馬數數論轉換中，取formula_85，可以證明在formula_86之下formula_77可以如下選取：(1) formula_88(2) formula_89在這兩種選取方式下，formula_8是2的冪次，所以基數-2的蝴蝶結構適用於此例中，而若formula_16是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為餘數，接上面的等式方程式化如下。formula_92formula_93formula_94formula_95整理一下formula_96formula_97formula_98formula_99最後的r 一定會變成1，所以我們只要不斷的將r乘上-q帶往下一個式子(像是r1*(-q1))，跟代往下下個式子(像是r3的左邊式子要帶入r1)即可求得最後我們想得到的 (1)，最後的N旁的係數就是反元素。def modInv(x, M): #by Euclidean algorithm - iterative versiont, new_t, r, new_r = 0, 1, M, xwhile new_r != 0:quotient = int(r / new_r)tmp_new_t = new_ttmp_t = ttmp_new_r = new_rtmp_r = rt, new_t = new_t, (t - quotient * new_t)r, new_r = new_r, (r % new_r)if r &gt; 1:print("Inverse doesn't exist")if t  0:if power &amp; 1:result = (result * base) % Mbase = (base * base) % Mpower = power » 1return result演演算法.以下程式預設。import numpy as np這邊只要再從1到M-1中選出一個適當的alpha (Root Of Unity)，其滿足"轉換公式"段落的(5)即可。def NTT(x,N,M):#TODO: RootOfUnityalpha = 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) % MA[j][i] = A[i][j]return np.matmul(A,x) % M, alphadef 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) % MB[j][i] = B[i][j]B = B * NInvreturn 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.