【密码学算法实现】数论基础

1.1 有限域上的四则运算

  在密码学当中，有限域通常特指 $GF(2^m)$ 域，因为这类域有较为良好的性质，便于复杂密码学算法的硬件实现。

1.1.1 有限域上的加减法

  $GF(2^m)$域实质上由二元域扩张形成，因此域中元素的加减法符合实数域中的模2加减法，即 $a+b(\mathrm{mod}2)$。在模2加减法中，$-1\equiv 1(\mathrm{mod}2)$，因此减法于加法的结果相同。同时，可以观察到 $0+0\equiv 0(\mathrm{mod}2)$，$1+0\equiv 1(\mathrm{mod}2)$，$1+1\equiv 0(\mathrm{mod}2)$，这与异或运算的结果相同。所以 $GF(2^m)$ 域上的加减法等同于两数异或运算，即 $a+b(\mathrm{mod}2)=a\oplus b$。

1.1.2 有限域上的乘法

  以AES算法常用的 $GF(2^8)$ 域为例，域中任意元素均可由下表所示的一组数异或组合而成，该组数可以看作这个域的一个基。

0x01(0b00000001)        0x10(0b00010000)
0x02(0b00000010)        0x20(0b00010000)
0x04(0b00000100)        0x40(0b01000000)
0x08(0b00001000)        0x80(0b10000000)

  取 $GF(2^8)$ 域中的元素0x78和0x39，设 $x=0x78$，$y=0x39$，不可约多项式 $p(x)=0x11b=x^8+x^4+x^3+x+1$ 那么有 $y=0b00111001=0x20\oplus 0x10\oplus 0x08\oplus 0x01$，由有限域内的分配律可得：
$$\begin{array}{rl}x\times y=& x\times (0x20\oplus 0x10\oplus 0x08\oplus 0x01)(\mathrm{mod}0x11b)\\ =& (x\times 0x20(\mathrm{mod}0x11b))\oplus \\ & (x\times 0x10(\mathrm{mod}0x11b))\oplus \\ & (x\times 0x08(\mathrm{mod}0x11b))\oplus \\ & (x\times 0x01(\mathrm{mod}0x11b))\end{array}$$
由此可知，若得到下表所示的一组关于 $x$ 的数，则可以计算出 $x\times y$ 的值。

x·0x01 (mod p(x))      x·0x10 (mod p(x))
x·0x02 (mod p(x))      x·0x20 (mod p(x))
x·0x04 (mod p(x))      x·0x40 (mod p(x))
x·0x08 (mod p(x))      x·0x80 (mod p(x))

上述过程的主要代码如下所示：

# 以GF(2^8)域为例
def get_multi(x, p):
    x <<= 1
    if x & 0x100 == 0x100:
        x ^= p
    return x & 0xff
def Multi(x, y, p):
    result = 0
    while y > 0:
        if y & 1 == 1:
            result ^= x
        
        # x乘以2，同时进行越界取模不可约多项式处理
        x = get_multi(x, p) 
        y >>= 1
    return result

1.1.3 有限域上的除法

  有限域上的除法运算实际上可看作对除法竖式的模拟。取 $GF(2^8)$ 域中的元素0x78和0x09，设被除数 $x=0x78$，除数 $y=0x09$，则有如下图模拟过程：

设被除数 $x$ 的二进制长度为 $le{n}_1$，除数 $y$ 的二进制长度为 $le{n}_2$，则每次将 $x$ 刷新为 $x=x-(y\ll (le{n}_1-le{n}_2))(\mathrm{mod}2)=x\oplus (y\ll (le{n}_1-le{n}_2))$，直到 $le{n}_1\le le{n}_2$ 为止。

上述过程的主要代码如下所示：

def Div(x, y):
    len1, len2 = len(bin(x)), len(bin(y))
    len_t = len1 - len2 + 1
    if len_t < 1:
        return 0, x
    elif len_t == 1:
        return 1, x ^ y
    else:
        ans = 0
        while len1 >= len2:
            x ^= y << (len1 - len2)
            ans ^= 1 << (len1 - len2)
            len1 = len(bin(x))
        return ans, x #ans为商，x为余数

注：上述代码中的方法名与后文一致

1.2 Eucilid 算法

   Eucilid算法又称辗转相除法，基于(1)式（Bezout等式），可以求出整数$a$，$b$的最大公因数、$a$模$b$的逆元$s$和$b$模$a$的逆元$t$。
$$sa+tb=gcd\text{\hspace{1em}(1)}$$
   在下面的一系列运算中可以看到Eucilid算法的流程，设$a\ge b$。第i轮带余除法的除数作为第i+1轮带余除法的被除数，第i轮带余除法的余数作为第i+1轮带余除法的除数，这就是辗转的核心。
$$\begin{array}{rlrl}a& =q_{1}b+r_1& 0<& r_1<b\\ b& =q_2{r}_1+r_2& 0<& r_2<r_1\\ & \cdots & & \cdots \\ r_{n-2}& =q_n{r}_{n-1}+r_n& 0<& r_n<r_{n-1}\\ r_{n-1}& =q_n{r}_n+0& & \\ gcd(a,b)& =r_n& & \end{array}$$
其主要代码如下所示：

#实数域
def euc_div(a, b):
    if b == 0:
        return 1, 0, a
    x, y, gcd = euc_div(b, a % b)
    x, y = y, (x - (a // b) * y)
    return x, y, gcd

#GF(2^m域)
def euc_div(a, b, p):
    if b == 0:
        return 1, 0, a
    div, rem = Div(a, b)
    x, y, gcd = euc_div(b, rem, p)
    x, y = y, x ^ Multi(div, y, p)
    return x, y, gcd

1.3 快速幂取模算法

   当幂和模数增大时，幂取模的运算复杂度呈指数倍增长，在幂和模数达到一定数量级（例如RSA算法中常用的1024bit）时，普通幂运算的运算效率远远低于所需效率，此时则需要更高效的幂取模运算。
   快速幂取模运算核心思想在于把指数转化为二进制序列，将大量复杂的取模运算转化为少量的取模运算和大量简便的乘幂运算，提升运算效率也适合硬件实现。举个例子：
$$\begin{array}{rl}base& =1234,power=17=(10001{)}_2,mod=2048\\ & \\ 1234^{17}& =1234^{1\times 2^0}\times 1234^{0\times 2^1}\times 1234^{0\times 2^2}\times 1234^{0\times 2^3}\times 1234^{1\times 2^4}\\ & =1234^{2^0}\times 1\times 1\times 1\times 1234^{2^4}\\ & =(1234^0{)}^2\times 1\times 1\times 1\times (1234^4{)}^2\\ & \\ 1234^{17}& (\mathrm{mod}2048)=(1234^0{)}^2(\mathrm{mod}2048)\times (1234^4{)}^2(\mathrm{mod}2048)\end{array}$$
其主要代码如下所示：

#实数域
def quick_mod(base, power, mod):
    ans = 1
    while power > 0:
        if power & 1 == 1:
            ans = (ans * base) % mod
        power >>= 1
        base = (base * base) % mod
    return ans

#GF(2^m域)
def quick_mod(base, power, mod):
    ans = 1
    while power > 0:
        if power & 1 == 1:
            ans = Multi(ans, base, mod)
        power >>= 1
        base = Multi(base, base, mod)
    return ans

1.4 中国剩余定理

  中国剩余定理（Chinese Remainder Theorem, CRT）可以将一个大整数M转化为若干个剩余类同构，如(2)式所示。
$$M=\prod _{i=1}^k{m}_i(m_i互素)\text{\hspace{1em}(2)}$$
  中国剩余定理最广泛的应用是解同余方程组，例如《孙子算经》中一个经典问题：“今有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二，问物几何？”这个问题可以表示成(3)式的方程组：
$$\{\begin{array}{l}x\equiv 2(\mathrm{mod}3)\\ x\equiv 3(\mathrm{mod}5)\\ x\equiv 2(\mathrm{mod}7)\end{array}\text{\hspace{1em}(3)}$$
式中 $m_1=3$、$m_2=5$、$m_3=7$，$M=3\times 5\times 7=105$。则有(4)式：
$$\{\begin{array}{l}{M}_1=\frac{M}{m_1}=35M_1^{-1}=M_1(\mathrm{mod}{m}_1)=2\\ M_2=\frac{M}{m_2}=21M_2^{-1}=M_2(\mathrm{mod}{m}_2)=1\\ M_3=\frac{M}{m_3}=15M_3^{-1}=M_3(\mathrm{mod}{m}_3)=1\end{array}\text{\hspace{1em}(4)}$$
由此可得方程式的解：
$$\begin{array}{rl}x& \equiv a_1{M}_1{M}_1^{-1}+a_2{M}_2{M}_2^{-1}+a_3{M}_3{M}_3^{-1}(\mathrm{mod}M)\\ & \equiv 2\times 35\times 2+3\times 21\times 1+2\times 15\times 1(\mathrm{mod}105)\\ & \equiv 233(\mathrm{mod}105)\\ & =23\end{array}$$
上述过程的主要代码如下所示：

def CRT(m, c):
    M_sum, x = 1, 0
    M_i, t_i = [], []
    for num in m:
        M_sum *= int(num)
    for num in m:
        M_i.append(M_sum//int(num))
        t_i.append(euc_div(M_sum//int(num), int(num))[0])
    for i in range(len(m)):
        x += int(c[i]) * t_i[i] * M_i[i]
    if x % M_sum > 0:
        return x % M_sum
    else:
        return M_sum