数论转换

数论转换是一种计算折积的快速算法。计算折积的快速算法中最常用的一种是使用快速傅里叶变换，然而快速傅立叶变换具有一些实现上的缺点，举例来说，资料向量必须乘上复数系数的矩阵加以处理，而且每个复数系数的实部和虚部是一个正弦及余弦函数，因此大部分的系数都是浮点数，也就是说，必须做复数而且是浮点数的运算，因此计算量会比较大，而且浮点数运算产生的误差会比较大。

而在数论转换中，资料向量需要乘上的矩阵是一个整数的矩阵，因此不需要作浮点数的乘法运算，更进一步，在模数为$M$的情况下，只有$M^2$种可能的加法与$M^2$种可能的乘法，这可以使用内存把这些有限数目的加法和乘法存起来，利用查表法来计算，使得数论转换完全不须任何乘法与加法运算，不过需要较大的内存与消耗较大的存取内存量。

虽然使用数论转换可以降低计算折积的复杂度，不过由于只进行整数的运算，所以只能用于对整数的信号计算折积，而利用快速傅立叶变换可以对任何复数的信号计算折积，这是降低运算复杂度所要付出的代价。

## 转换公式

数论转换的转换式为

$$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}\sum_{k=0}^{N-1} F[k]\alpha^{-nk}\pmod{M}\quad n=0,1,2,\ldots ,N-1$$

注解：

(1) $M$是一个质数。

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

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

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

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

举一个例子：

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

正转换$$\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}$$

反转换$$\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) 正交性质

数论转换矩阵的每个列是互相正交的，即$\sum_{n=0}^{N-1}\alpha^{nl}\alpha^{-nk}=\begin{cases} N \quad if\ k=l\\ 0 \quad if\ k\ne l\end{cases} \pmod{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]\alpha^{-nk}=N^{-1}\sum_{-k=0}^{N-1}F[-k]\alpha^{nk}$，即$N^{-1}F[-k]$的数论转换。

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

步骤二：求$F[-k]$的数论转换。

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

(4) 线性性质

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

(5) 平移性质

若$f[n]\leftrightarrow F[k]$，则$f[n+l]\leftrightarrow \alpha^{-lk}F[k]\pmod{M}$，而且$f[n]\alpha^{nl}\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]\pmod{M}$，而且$f[n]g[n]\leftrightarrow \frac{1}{N}F[k]\otimes G[k]\pmod{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]\pmod{M}$，而且$N\sum_{n=0}^{N-1}f^2[n]=\sum_{k=0}^{N-1}F[k]F[-k]\pmod{M}$

## 快速数论转换

如果转换点数N是一个2的次方，则可以使用类似基数-2快速傅立叶变换的蝴蝶结构来有效率的实现数论转换。同样的互质因子算法也可以应用在数论转换上。

$$\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^kH[k]\quad 0\le k\le \frac{N}{2}-1 \\ G[k-\frac{N}{2}]-\alpha^kH[k-\frac{N}{2}]\quad \frac{N}{2}\le k\le N-1 \end{cases}\end{matrix}$$

其中，$G[k]=\sum_{r=0}^{N/2-1}f[2r](\alpha^2)^{rk}$，$H[k]=\sum_{r=0}^{N/2-1}f[2r+1](\alpha^2)^{rk}$。

因此一个$N$点数论转换可以拆解成两个$\frac{N}{2}$点的数论转换。

## N, M, alpha的选择

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

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

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

## 特殊的数论转换

### 梅森质数数论转换

梅森质数是指形如$2^p-1$的质数记做$M_p$，其中$p$是一个质数。

在梅森质数数论转换中，取$M=M_p$，可以证明$\alpha , N$可以如下选取：

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

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

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

### 费马数数论转换

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

在费马数数论转换中，取$M=F_t$，可以证明在$0<t\le 4$之下$\alpha , N$可以如下选取：

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

(2) $\alpha = 3, N=F_t-1$

在这两种选取方式下，$N$是2的幂次，所以基数-2的蝴蝶结构适用于此例中，而若$\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},$$

$$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旁的系数就是反元素。

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

gcd(N, M mod N)=1，则存在$x_{1}N + y_{1}(M \, mod \,N) = 1$

$$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$$

$$y_{1} =x$$

可以借由下一层的系数来推出上一层的系数。

```python3
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} \, mod \, M = 1$

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

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

搭配快速mod可以容易的算出反元素，power是偶数时则可以用$(a^{2} \,mod\,M)^{power/2}$; power是基数时则可以用$(aN \,mod\,M)^{power-1}$，让power变成偶数。反复直到power变成1。

```python3
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
```

## 算法

以下程式预设。

```python3
import numpy as np
```

Root of Unity 是从1到M-1中选出一个适当的 alpha，其满足"转换公式"段落的(5)。

```python3
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
```

```python
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
```

```python3
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
```