defmodInv(a,m):#fermat little thmreturnmodExponent(a,m-2,m)defmodExponent(base,power,M):#quick mod O(log(power))result=1power=int(power)base=base%Mwhilepower>0:ifpower&1:result=(result*base)%Mbase=(base*base)%Mpower=power>>1returnresult
演算法
以下程式預設。
importnumpyasnp
這邊只要再從1到M-1中選出一個適當的alpha (Root Of Unity),其滿足"轉換公式"段落的(5)即可。
[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.
四月 01, 2023
數論轉換, 是一種計算摺積的快速演算法, 計算摺積的快速演算法中最常用的一種是使用快速傅里葉變換, 然而快速傅立葉變換具有一些實現上的缺點, 舉例來說, 資料向量必須乘上複數係數的矩陣加以處理, 而且每個複數係數的實部和虛部是一個正弦及餘弦函數, 因此大部分的係數都是浮點數, 也就是說, 必須做複數而且是浮點數的運算, 因此計算量會比較大, 而且浮點數運算產生的誤差會比較大, 而在中, 資料向量需要乘上的矩陣是一個整數的矩陣, 因此不需要作浮點數的乘法運算, 更進一步, 在模數為m, displaystyle, 的. 數論轉換是一種計算摺積的快速演算法 計算摺積的快速演算法中最常用的一種是使用快速傅里葉變換 然而快速傅立葉變換具有一些實現上的缺點 舉例來說 資料向量必須乘上複數係數的矩陣加以處理 而且每個複數係數的實部和虛部是一個正弦及餘弦函數 因此大部分的係數都是浮點數 也就是說 必須做複數而且是浮點數的運算 因此計算量會比較大 而且浮點數運算產生的誤差會比較大 而在數論轉換中 資料向量需要乘上的矩陣是一個整數的矩陣 因此不需要作浮點數的乘法運算 更進一步 在模數為M displaystyle M 的情況下 只有M 2 displaystyle M 2 種可能的加法與M 2 displaystyle M 2 種可能的乘法 這可以使用記憶體把這些有限數目的加法和乘法存起來 利用查表法來計算 使得數論轉換完全不須任何乘法與加法運算 不過需要較大的記憶體與消耗較大的存取記憶體量 雖然使用數論轉換可以降低計算摺積的複雜度 不過由於只進行整數的運算 所以只能用於對整數的信號計算摺積 而利用快速傅立葉變換可以對任何複數的信號計算摺積 這是降低運算複雜度所要付出的代價 目录 1 轉換公式 2 數論轉換的性質 3 快速數論轉換 4 N M alpha的選擇 5 特殊的數論轉換 5 1 梅森質數數論轉換 5 2 費馬數數論轉換 6 實作議題 7 演算法 8 參考資料轉換公式 编辑數論轉換的轉換式為F k n 0 N 1 f n a n k mod M k 0 1 2 N 1 displaystyle F k sum n 0 N 1 f n alpha nk pmod M quad k 0 1 2 ldots N 1 而數論轉換的反轉換式為f n N 1 k 0 N 1 F k a n k mod M n 0 1 2 N 1 displaystyle f n N 1 sum k 0 N 1 F k alpha nk pmod M quad n 0 1 2 ldots N 1 註解 1 M displaystyle M 是一個質數 2 mod M displaystyle pmod M 表示除以M的餘數 3 N displaystyle N 必須是M 1 displaystyle M 1 的因數 當N 1 displaystyle N neq 1 時N displaystyle N 和M displaystyle M 互質 4 N 1 displaystyle N 1 是N displaystyle N 對模數M displaystyle M 之模反元素 5 a displaystyle alpha 為模M的N階單位根 即a N 1 mod M displaystyle alpha N 1 pmod M 而且a k 1 mod M k 1 2 N 1 displaystyle alpha k neq 1 pmod M k 1 2 ldots N 1 若此時N M 1 displaystyle N M 1 我們稱a displaystyle alpha 為模M的原根舉一個例子 一個M 5 displaystyle M 5 的N 4 displaystyle N 4 點數論轉換與反轉換如下 取a 2 displaystyle alpha 2 注意此時N 1 4 displaystyle N 1 4 正轉換 F 0 F 1 F 2 F 3 1 1 1 1 1 2 4 3 1 4 1 4 1 3 4 2 f 0 f 1 f 2 f 3 mod M displaystyle begin pmatrix F 0 F 1 F 2 F 3 end pmatrix begin pmatrix 1 amp 1 amp 1 amp 1 1 amp 2 amp 4 amp 3 1 amp 4 amp 1 amp 4 1 amp 3 amp 4 amp 2 end pmatrix begin pmatrix f 0 f 1 f 2 f 3 end pmatrix pmod M 反轉換 f 0 f 1 f 2 f 3 4 4 4 4 4 2 1 3 4 1 4 1 4 3 1 2 F 0 F 1 F 2 F 3 mod M displaystyle begin pmatrix f 0 f 1 f 2 f 3 end pmatrix begin pmatrix 4 amp 4 amp 4 amp 4 4 amp 2 amp 1 amp 3 4 amp 1 amp 4 amp 1 4 amp 3 amp 1 amp 2 end pmatrix begin pmatrix F 0 F 1 F 2 F 3 end pmatrix pmod M 數論轉換的性質 编辑 1 正交性質數論轉換矩陣的每個列是互相正交的 即 n 0 N 1 a n l a n k N i f k l 0 i f k l mod M displaystyle sum n 0 N 1 alpha nl alpha nk begin cases N quad if k l 0 quad if k neq l end cases pmod M 2 對稱性若f n f N n displaystyle f n f N n 則f n displaystyle f n 的數論轉換F k displaystyle F k 也會有F k F N k displaystyle F k F N k 的特性 若f n f N n displaystyle f n f N n 則f n displaystyle f n 的數論轉換F k displaystyle F k 也會有F k F N k displaystyle F k F N k 的特性 3 反數論轉換可由正數論轉換實現f n N 1 k 0 N 1 F k a n k N 1 k 0 N 1 F k a n k displaystyle f n N 1 sum k 0 N 1 F k alpha nk N 1 sum k 0 N 1 F k alpha nk 即N 1 F k displaystyle N 1 F k 的數論轉換 步驟一 把F k displaystyle F k 改寫成F k displaystyle F k 若F k F 0 F 1 F 2 F N 1 displaystyle F k F 0 F 1 F 2 ldots F N 1 則F k F 0 F N 1 F 2 F 1 displaystyle F k F 0 F N 1 ldots F 2 F 1 步驟二 求F k displaystyle F k 的數論轉換 步驟三 乘上N 1 displaystyle N 1 4 線性性質若f n F k displaystyle f n leftrightarrow F k g n G k displaystyle g n leftrightarrow G k displaystyle leftrightarrow 表互為數論轉換對 則有a f n b g n a F k b g k mod M displaystyle af n bg n leftrightarrow aF k bg k pmod M 5 平移性質若f n F k displaystyle f n leftrightarrow F k 則f n l a l k F k mod M displaystyle f n l leftrightarrow alpha lk F k pmod M 而且f n a n l F k l displaystyle f n alpha nl leftrightarrow F k l 6 圓周摺積性質若f n F k displaystyle f n leftrightarrow F k g n G k displaystyle g n leftrightarrow G k 則f n g n F k G k mod M displaystyle f n otimes g n leftrightarrow F k G k pmod M 而且f n g n 1 N F k G k mod M displaystyle f n g n leftrightarrow frac 1 N F k otimes G k pmod M displaystyle otimes 代表圓形摺積 這個性質是數論轉換可以用來作為摺積的快速演算法的主要原因 7 帕塞瓦爾定理 Parseval s Theorem 若f n F k displaystyle f n leftrightarrow F k 則N n 0 N 1 f n f n k 0 N 1 F 2 k mod M displaystyle N sum n 0 N 1 f n f n sum k 0 N 1 F 2 k pmod M 而且N n 0 N 1 f 2 n k 0 N 1 F k F k mod M displaystyle N sum n 0 N 1 f 2 n sum k 0 N 1 F k F k pmod M 快速數論轉換 编辑如果轉換點數N是一個2的次方 則可以使用類似基數 2快速傅立葉變換的蝴蝶結構來有效率的實現數論轉換 同樣的互質因子算法也可以應用在數論轉換上 F k n 0 N 1 f n a n k r 0 N 2 1 f 2 r a 2 r k r 0 N 2 1 f 2 r 1 a 2 r 1 k r 0 N 2 1 f 2 r a 2 r k a k r 0 N 2 1 f 2 r 1 a 2 r k G k a k H k 0 k N 2 1 G k N 2 a k H k N 2 N 2 k N 1 displaystyle begin matrix F k amp amp sum n 0 N 1 f n alpha nk amp amp sum r 0 N 2 1 f 2r alpha 2rk sum r 0 N 2 1 f 2r 1 alpha 2r 1 k amp amp sum r 0 N 2 1 f 2r alpha 2 rk alpha k sum r 0 N 2 1 f 2r 1 alpha 2 rk amp amp begin cases G k alpha k H k quad 0 leq k leq frac N 2 1 G k frac N 2 alpha k H k frac N 2 quad frac N 2 leq k leq N 1 end cases end matrix 其中 G k r 0 N 2 1 f 2 r a 2 r k displaystyle G k sum r 0 N 2 1 f 2r alpha 2 rk H k r 0 N 2 1 f 2 r 1 a 2 r k displaystyle H k sum r 0 N 2 1 f 2r 1 alpha 2 rk 因此一個N displaystyle N 點數論轉換可以拆解成兩個N 2 displaystyle frac N 2 點的數論轉換 N M alpha的選擇 编辑由於數論轉換可以擁有類似快速傅立葉變換的快速演算法 因此通常N displaystyle N 會選擇適合使用快速演算的值 比如說取為2的冪次 或是可以分解成許多小質數相乘的數 在數論轉換中 需要大量地和a displaystyle alpha 的冪次做乘法 因此 如果可以取a displaystyle alpha 為2或2的冪次 則每一次的乘法在二進制中只會是一個移位的操作 可以省下大量的乘法運算 因為要做模M displaystyle M 的運算 所以M displaystyle M 的二進位表示法中 1的個數越少越好 同時M displaystyle M 的值不能取太小 這是因為數論轉換後的值都小於M displaystyle M 因此當真實的摺積的結果會大於M displaystyle M 時就會發生錯誤 所以必須謹慎選取M displaystyle M 的大小 特殊的數論轉換 编辑梅森質數數論轉換 编辑 梅森質數是指形如2 p 1 displaystyle 2 p 1 的質數記做M p displaystyle M p 其中p displaystyle p 是一個質數 在梅森質數數論轉換中 取M M p displaystyle M M p 可以證明a N displaystyle alpha N 可以如下選取 1 a 2 N p N 1 M p M p 1 p displaystyle alpha 2 N p N 1 M p frac M p 1 p 2 a 2 N 2 p N 1 M p M p 1 2 p displaystyle alpha 2 N 2p N 1 M p frac M p 1 2p 在這兩種選取方式下 由於a displaystyle alpha 是2的冪次 所以只需移位運算 但N displaystyle N 不是2的冪次 所以基數 2的蝴蝶結構不適用於此例中 同時N displaystyle N 為質數或一個質數的兩倍 並不是一個可以拆解成很多質因數乘積的數 因此也無法由互質因子演算法得到太大的好處 梅森質數數論轉換通常用於較短序列的摺積運算 費馬數數論轉換 编辑 費馬數是指形如2 2 t 1 displaystyle 2 2 t 1 的數記做F t displaystyle F t 在費馬數數論轉換中 取M F t displaystyle M F t 可以證明在0 lt t 4 displaystyle 0 lt t leq 4 之下a N displaystyle alpha N 可以如下選取 1 a 2 N 2 t 1 displaystyle alpha 2 N 2 t 1 2 a 3 N F t 1 displaystyle alpha 3 N F t 1 在這兩種選取方式下 N displaystyle N 是2的冪次 所以基數 2的蝴蝶結構適用於此例中 而若a displaystyle alpha 是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 displaystyle M q 0 N r 1 N q 1 r 1 r 2 displaystyle N q 1 r 1 r 2 r 1 q 2 r 2 r 3 displaystyle r 1 q 2 r 2 r 3 r 2 q 3 r 3 r 4 displaystyle r 2 q 3 r 3 r 4 整理一下M q 0 N r 1 displaystyle M q 0 N r 1 N q 1 r 1 r 2 displaystyle N q 1 r 1 r 2 r 1 q 2 r 2 r 3 displaystyle r 1 q 2 r 2 r 3 r 2 q 3 r 3 r 4 displaystyle 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 gt 1 print Inverse doesn t exist if t lt 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 displaystyle xM yN 1 gcd N M mod N 1 則存在x 1 N y 1 M m o d N 1 displaystyle x 1 N y 1 M mod 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 displaystyle x 1 N y 1 M mod 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 displaystyle x 1 q 1 y 1 y y 1 x displaystyle 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 displaystyle N M 1 mod M 1 N N M 2 m o d M m o d M 1 displaystyle N N M 2 mod M mod M 1 顯然N M 2 m o d M displaystyle N M 2 mod M 就是N的反元素 搭配快速mod可以容易的算出反元素 power是偶數時則可以用 a 2 m o d M p o w e r 2 displaystyle a 2 mod M power 2 power是基數時則可以用 a N m o d M p o w e r 1 displaystyle aN mod M power 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 gt 0 if power amp 1 result result base M base base base M power power gt gt 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 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 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 數論轉換 amp oldid 74418034, 维基百科,wiki,书籍,书籍,图书馆,