数论变换

数论变换的变换式为

${\displaystyle F[k]=\sum _{n=0}^{N-1}f[n]\alpha ^{nk}{\pmod {M}}\quad k=0,1,2,\ldots ,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) ${\displaystyle M}$ 是一个质数。

(2) ${\displaystyle {\pmod {M}}}$  表示除以M的余数

(3) ${\displaystyle N}$ 必须是${\displaystyle M-1}$ 的因数。(当${\displaystyle N\neq 1}$ 时${\displaystyle N}$ 和${\displaystyle M}$ 互质)

(4)${\displaystyle N^{-1}}$ 是${\displaystyle N}$ 对模数${\displaystyle M}$ 之模反元素

(5)${\displaystyle \alpha }$ 为模M的N阶单位根，即${\displaystyle \alpha ^{N}=1{\pmod {M}}}$ 而且${\displaystyle \alpha ^{k}\neq 1{\pmod {M}}\ k=1,2,\ldots ,N-1}$ 。若此时${\displaystyle N=M-1}$ ，我们称${\displaystyle \alpha }$ 为模M的原根

举一个例子：

一个${\displaystyle M=5}$ 的${\displaystyle N=4}$ 点数论变换与逆变换如下，取${\displaystyle \alpha =2}$ ，注意此时${\displaystyle N^{-1}=4}$ 

正变换 ${\displaystyle {\begin{pmatrix}F[0]\\F[1]\\F[2]\\F[3]\end{pmatrix}}={\begin{pmatrix}1&1&1&1\\1&2&4&3\\1&4&1&4\\1&3&4&2\end{pmatrix}}{\begin{pmatrix}f[0]\\f[1]\\f[2]\\f[3]\end{pmatrix}}{\pmod {M}}}$ 

逆变换 ${\displaystyle {\begin{pmatrix}f[0]\\f[1]\\f[2]\\f[3]\end{pmatrix}}={\begin{pmatrix}4&4&4&4\\4&2&1&3\\4&1&4&1\\4&3&1&2\end{pmatrix}}{\begin{pmatrix}F[0]\\F[1]\\F[2]\\F[3]\end{pmatrix}}{\pmod {M}}}$ 

(1) 正交性质

数论变换矩阵的每个列是互相正交的，即${\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) 对称性

若${\displaystyle f[n]=f[N-n]}$ ，则${\displaystyle f[n]}$ 的数论变换${\displaystyle F[k]}$ 也会有${\displaystyle F[k]=F[N-k]}$ 的特性。

若${\displaystyle f[n]=-f[N-n]}$ ，则${\displaystyle f[n]}$ 的数论变换${\displaystyle F[k]}$ 也会有${\displaystyle F[k]=-F[N-k]}$ 的特性。

(3) 反数论变换可由正数论变换实现

${\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}}$ ，即${\displaystyle N^{-1}F[-k]}$ 的数论变换。

步骤一：把${\displaystyle F[k]}$ 改写成${\displaystyle F[-k]}$ ，若${\displaystyle F[k]=F_{0},F_{1},,F_{2}\ldots ,F_{N-1}}$ ，则${\displaystyle F[-k]=F_{0},F_{N-1},\ldots ,F_{2},F_{1}}$ 

步骤二：求${\displaystyle F[-k]}$ 的数论变换。

步骤三：乘上${\displaystyle N^{-1}}$ 。

(4) 线性性质

若${\displaystyle f[n]\leftrightarrow F[k]}$ ，${\displaystyle g[n]\leftrightarrow G[k]}$ ，(${\displaystyle \leftrightarrow }$ 表互为数论变换对)则有${\displaystyle af[n]+bg[n]\leftrightarrow aF[k]+bg[k]{\pmod {M}}}$ 。

(5) 平移性质

若${\displaystyle f[n]\leftrightarrow F[k]}$ ，则${\displaystyle f[n+l]\leftrightarrow \alpha ^{-lk}F[k]{\pmod {M}}}$ ，而且${\displaystyle f[n]\alpha ^{nl}\leftrightarrow F[k+l]}$ 。

(6) 循环卷积性质

若${\displaystyle f[n]\leftrightarrow F[k]}$ ，${\displaystyle g[n]\leftrightarrow G[k]}$ ，则${\displaystyle f[n]\otimes g[n]\leftrightarrow F[k]G[k]{\pmod {M}}}$ ，而且${\displaystyle f[n]g[n]\leftrightarrow {\frac {1}{N}}F[k]\otimes G[k]{\pmod {M}}}$ 。(${\displaystyle \otimes }$ 代表圆形卷积。)

这个性质是数论变换可以用来作为卷积的快速算法的主要原因。

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

若${\displaystyle f[n]\leftrightarrow F[k]}$ ，则${\displaystyle N\sum _{n=0}^{N-1}f[n]f[-n]=\sum _{k=0}^{N-1}F^{2}[k]{\pmod {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快速傅里叶变换的蝴蝶结构来有效率的实现数论变换。同样的互质因子算法也可以应用在数论变换上。

${\displaystyle {\begin{matrix}F[k]&=&\sum _{n=0}^{N-1}f[n]\alpha ^{nk}\\\ &=&\sum _{r=0}^{N/2-1}f[2r]\alpha ^{2rk}+\sum _{r=0}^{N/2-1}f[2r+1]\alpha ^{(2r+1)k}\\\ &=&\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}\\\ &=&{\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}}}$ 

其中，${\displaystyle G[k]=\sum _{r=0}^{N/2-1}f[2r](\alpha ^{2})^{rk}}$ ，${\displaystyle H[k]=\sum _{r=0}^{N/2-1}f[2r+1](\alpha ^{2})^{rk}}$ 。 因此一个${\displaystyle N}$ 点数论变换可以拆解成两个${\displaystyle {\frac {N}{2}}}$ 点的数论变换。

由于数论变换可以拥有类似快速傅里叶变换的快速算法，因此通常${\displaystyle N}$ 会选择适合使用快速演算的值，比如说取为2的幂次，或是可以分解成许多小质数相乘的数。

在数论变换中，需要大量地和${\displaystyle \alpha }$ 的幂次做乘法，因此，如果可以取${\displaystyle \alpha }$ 为2或2的幂次，则每一次的乘法在二进制中只会是一个移位的操作，可以省下大量的乘法运算。

因为要做模${\displaystyle M}$ 的运算，所以${\displaystyle M}$ 的二进制表示法中，1的个数越少越好，同时${\displaystyle M}$ 的值不能取太小，这是因为数论变换后的值都小于${\displaystyle M}$ ，因此当真实的卷积的结果会大于${\displaystyle M}$ 时就会发生错误，所以必须谨慎选取${\displaystyle M}$ 的大小。

梅森质数数论变换

梅森质数是指形如${\displaystyle 2^{p}-1}$ 的质数记做${\displaystyle M_{p}}$ ，其中${\displaystyle p}$ 是一个质数。

在梅森质数数论变换中，取${\displaystyle M=M_{p}}$ ，可以证明${\displaystyle \alpha ,N}$ 可以如下选取：

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

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

在这两种选取方式下，由于${\displaystyle \alpha }$ 是2的幂次，所以只需移位运算，但${\displaystyle N}$ 不是2的幂次，所以基数-2的蝴蝶结构不适用于此例中，同时${\displaystyle N}$ 为质数或一个质数的两倍，并不是一个可以拆解成很多质因数乘积的数，因此也无法由互质因子算法得到太大的好处。梅森质数数论变换通常用于较短序列的卷积运算

费马数数论变换

费马数是指形如${\displaystyle 2^{2^{t}}+1}$ 的数记做${\displaystyle F_{t}}$ 。

在费马数数论变换中，取${\displaystyle M=F_{t}}$ ，可以证明在${\displaystyle 0<t\leq 4}$ 之下${\displaystyle \alpha ,N}$ 可以如下选取：

(1) ${\displaystyle \alpha =2,N=2^{t+1}}$ 

(2) ${\displaystyle \alpha =3,N=F_{t}-1}$ 

在这两种选取方式下，${\displaystyle N}$ 是2的幂次，所以基数-2的蝴蝶结构适用于此例中，而若${\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为余数，接上面的等式方程式化如下。

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

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

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

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

整理一下

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

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

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

${\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 > 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来看，我们知道存在${\displaystyle xM+yN=1}$  。

gcd(N, M mod N)=1，则存在${\displaystyle x_{1}N+y_{1}(M\,mod\,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的系数，这边可以得到一个递回关系，

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

${\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，${\displaystyle N^{M-1}\,mod\,M=1}$ 

${\displaystyle N(N^{M-2}\,mod\,M)\,mod\,M=1}$ 

显然${\displaystyle N^{M-2}\,mod\,M}$ 就是N的反元素。

搭配快速mod可以容易的算出反元素，power是偶数时则可以用${\displaystyle (a^{2}\,mod\,M)^{power/2}}$ ; power是基数时则可以用${\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 > 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.