2.6 蒙哥马利乘法

蒙哥马利空间乘法

  蒙哥马利算法的发明者不是二战时的英国陆军元帅蒙哥马利。外国嘛，同名同姓的人比较多，而这个蒙哥马利是Peter L. Montgomery,美国数学家，于1985年发明了蒙哥马利算法。
  蒙哥马利算法的重要概念是蒙哥马利空间，该空间的定义是三个数：x、r、n。n是要取模的数，x和r是两个乘数，$\overline{x}$是x在蒙哥马利空间的表示，按以下公式计算：
$$\overline{x}=x\cdot rmodn$$
  蒙哥马利空间的乘法，符号用$\ast$，比较特殊，按以下规则进行，需要注意的是本文中的$r^{-1}$不是传统数学里的$-1$次方，而是模乘逆元：
$$\begin{gathered} \overline{x}\ast \overline{y} \\ =\overline{x\cdot y} \\ =x\cdot y\cdot rmodn \\ =\overline{x}\cdot \overline{y}\cdot r^{-1}modn \\ =x\cdot r\cdot y\cdot r\cdot r^{-1}modn \end{gathered}$$
  举个例子：
$$\begin{gathered} n=997 \\ r=1024 \\ x=1024 \\ y=2048 \\ \overline{x}:=x\cdot rmodn=729 \\ \overline{y}:=y\cdot rmodn=461 \\ {r}^{-1}=517 \\ \overline{x}\ast \overline{y}:=\overline{x}\cdot \overline{y}\cdot r^{-1}modn=483 \end{gathered}$$
  需要注意的是蒙哥马利空间乘法不是蒙哥马利乘法。所以实际上$\overline{x}\ast \overline{y}$是$\overline{x}\cdot \overline{y}$的蒙哥马利空间表示。

蒙哥马利约化

  由蒙哥马利空间乘法，我们知道求两个数的乘积取模，或者叫模乘，有两种方法，第一种方法是直接相乘再取模，也就是使用公式$x\cdot y\cdot rmodn$，第二种方法是使用蒙哥马利空间乘法$\overline{x}\cdot \overline{y}\cdot r^{-1}modn$。
  从表面上看，蒙哥马利空间乘法没有任何优势啊，乘法次数一样多。但是有了蒙哥马利约化后就不一样了。蒙哥马利约化，就是一个由$x$求$\overline{x}$的算法，也就是由$\overline{x}\cdot \overline{y}$求$\overline{x\cdot y}$的算法。
  为了进行蒙哥马利约化Montgomery reduction，我们需要解一个方程。该方程如下：
$$r\cdot r^{-1}-n\cdot n^{\prime }=1$$
  当然，如果不知道$r^{-1}$时，可以将上述方程视为一个正号变成了负号的线性丢番图方程，可以通过扩展欧几里得法将$r^{-1}$和$n^{\prime }$求出来。
  有了这个就可以这样计算了：
$$\begin{gathered} x\cdot r^{-1} \\ =x\cdot r\cdot r^{-1}/r \\ =x\cdot (1-n\cdot n^{\prime })/r \\ \equiv x\cdot (1-n\cdot n^{\prime })/r+l\cdot n \\ =(x-x\cdot n\cdot n^{\prime })/r+l\cdot n \end{gathered}$$
  因为$x\cdot r^{-1}$同余于$(x-x\cdot n\cdot n^{\prime })/r+l\cdot n$，其中$l$是使得结果为小于n的整数的某个数。我们设置一个变量$q=x\cdot n^{\prime }$，上面的式子就可以改成：
$$\begin{gathered} x\cdot r^{-1} \\ \equiv (x-q\cdot n)/r+l\cdot n \end{gathered}$$
  又因为整个括号内是要除于r的，所以计算q的时候可以取模r的同余类，于是q可以是$q=(xmodr)\cdot n^{\prime }modr$。
  这样就给出了一个由$x$求$\overline{x}$的算法了，算法如下：
$$\begin{gathered} q=(xmodr)\cdot n^{\prime }modr \\ a=(x-q\cdot n)/r \\ whilea<0:a+=n \\ returna \end{gathered}$$

模幂乘优化

  蒙哥马利乘法的最大用途是模幂乘优化，也就是解决$a^{e}modn$这样的计算问题。首先对于幂运算，我们想到的是进行二进制快速幂算法，减少乘法次数。
  对于$a^{e}modn$，算法的步骤是首先计算a的蒙哥马利空间表示$\overline{a}$，然后计算$\overline{a^e}$，然后用蒙哥马利约简计算$\overline{a^e}\cdot r^{-1}$，就完了。不过需要注意的是两点：
  1、必须要满足互质才能解那个丢番图方程；
  2、在二进制快速幂过程中由于每次都乘以了r，所以需要利用蒙哥马利约化乘以$r^{-1}$以减少一个r。
  3、实际应用中的蒙哥马利，远比我说的上述过程复杂，因为用了更多的优化技术，但是我是讲蒙哥马利乘法原理，所以只简单地用了二进制快速幂。

Python实现

def extend_euclid(a: int, b: int):
    if a < b:
        x, y, g = extend_euclid(b, a)
        return y, x, g
    # 保存a/b的系数
    coefficients = []
    while b != 0:
        coefficients.append(a // b)
        a, b = b, a % b
    # 循环代入
    x, y = 1, 0
    coefficients.reverse()
    for i in coefficients:
        x, y = y, x - y * i
    return x, y, a


def diophantine(a, b, c):
    x_g, y_g, g = extend_euclid(a, b)
    if c % g != 0:
        raise 'No solutions'
    return x_g * c // g, y_g * c // g, b // g, -a // g


class Montgomery:

    def __init__(self, r, n):
        self.__r = r
        self.__n = n
        x0, y0, xk, yk = diophantine(r, n, 1)
        self.__r_ = x0
        self.__n_ = y0
        while self.__r_ < 0:
            if xk > 0:
                self.__r_ += xk
                self.__n_ += yk
                print(self.__r_)
            else:
                self.__r_ -= xk
                self.__n_ -= yk

                print(self.__r_, '-')

    def reduce(self, x):
        r, n, n_ = self.__r, self.__n, self.__n_
        q = (x % r) * n_ % r
        a = (x - q * n) // r
        while a < 0:
            a += n
        return a

    def present(self, x):
        return x * self.__r % self.__n

    @staticmethod
    def power(a, e, n):
        r = 2
        while r <= n:
            r <<= 1
        montgomery = Montgomery(r, n)
        # 先使用最笨的办法

        power = 1
        result = montgomery.present(1)
        p = montgomery.present(a)

        while power <= e:
            if power & e != 0:
                result = montgomery.reduce(result * p)
            p = montgomery.reduce(p * p)

            power = power << 1
        return montgomery.reduce(result)

Python测试代码

import unittest

import unittest

import modular_inverse as mi
from montgomery2 import Montgomery


class MyTestCase(unittest.TestCase):
    def test_montgomery_space(self):
        x, y, r, n = 1024, 2048, 1024, 997
        r_ = mi.modular_inverse(r, n)
        x_ = x * r % n
        print(f'r_={r_}')
        print(f"{x}*{r}%{n}={x_}")
        y_ = y * r % n
        print(f"{y}*{r}%{n}={y_}")
        print(f"{x_}*{y_}*{r_}%{n}={x_ * y_ * r_ % n}")

    def test_montgomery_pow(self):
        a, e, n = 3, 777, 17
        r1 = Montgomery.power(a, e, n)
        r2 = a ** e % n
        print(r1 == r2, r1, r2)

    def test_montgomery_reduce(self):
        r, n = 3, 4
        m = Montgomery(r, n)
        x, y = 1, 99
        x_, y_ = x * r % n, y * r % n
        xy_ = x * y * r % n

        print(m.reduce(x_ * y_), xy_)


if __name__ == '__main__':
    unittest.main()

  测试结果：

True 14 14