椭圆曲线和蒙哥马利算法的一些记录
0. 数论知识
以下有些地方的同余符号用了等号代替，或者省略了 \(\pmod{m}\) 的写法。另外，除非特别说明，代码都是用 Python 实现。
0.1. 同余（Modulo）
给定一个正整数 m，如果两个整数 a 和 b 分别除以 m 得到的余数相同，那么 a 和 b 对模 m 同余。记作 \(a \equiv b \pmod{m}\)
比如，7 和 29 对模 11 同余，\(7 \equiv 29 \pmod{11}\)
验算：
>>> 7 % 11
7
>>> 29 % 11
7
>>> (7 % 11) == (29 % 11)
True
0.2. 模逆（Modular Inverse）
对于非零整数 a，m，如果存在 b 使得 \(a \times b \equiv 1 \pmod{m}\)，那么就称 b 是 a 在模 m 下的逆，记作 \(a^{-1}\)
比如，对于模 11，在 (0, 11) 范围里，可以计算每一个数的逆：
>>> import itertools
>>>
>>> m: int = 11
>>> for a, b in itertools.product(range(1, m), repeat=2):
>>> if (a * b) % m == 1:
>>> print(f'{a}: {b}')
1: 1
2: 6
3: 4
4: 3
5: 9
6: 2
7: 8
8: 7
9: 5
10: 10
验算也很简单：\(1 \times 1 \equiv 1 \pmod{11}\)，\(2 \times 6 = 12 \equiv 1 \pmod{11}\)，……，\(10 \times 10 = 100 \equiv 1 \pmod{11}\)。
并且，在这个范围里，逆是唯一的。反证，如果 a 存在两个不相等的逆 \(b_{1}\) 和 \(b_{2}\)，那么 \(a \times b_{1} \equiv a \times b_{2} \pmod{11} \Longrightarrow a \times (b_{1} - b_{2}) \equiv 0 \pmod{11} \Longrightarrow b_{1} - b{2} \equiv 0 \pmod{11}\) 对于指定范围里的两个数，这里显然矛盾了。当然，在不指定范围时，如果存在逆，显然存在无穷多个。
0.3. 扩展欧几里德算法（Extended Euclidean）
扩展欧几里德算法用来求不定方程 \(ax + by = c\) 的解，其中 c 是 a 和 b 的最大公约数。
比如 \(7x+11y=1\)：
从 11 和 7 入手：
(1) 11 = 1 * 7 + 4 // 11 除以 7，余 4
(2) 7 = 1 * 4 + 3 // 7 除以 4，余 3
(3) 4 = 1 * 3 + 1 // 4 除以 3，余 1
(4) 3 = 3 * 1 // 3 除以 1，余 0，结束
从 (3) 开始：
1 = 4 - 1 * 3 // 由 (3) 移项得到
= (11 - 1 * 7) - 1 * (7 - 1 * 4) // 由 (1) 移项得到 4 = 11 - 1 * 7，并替换
// 由 (2) 移项得到 3 = 7 - 1 * 4，并替换
= 11 - 7 - 7 + 4 // 整理
= 11 - 2 * 7 + (11 - 1 * 7) // 替换
= 2 * 11 - 3 * 7 // 整理
所以 x0 = -3, y0 = 2。由于不定方程的解不唯一，可以继续求得通解公式 x = -3 + 11t, y = 2 - 7t
在通解公式基础上求逆：
x = -3 + 11t
inv(7) = (-3 + 11t) (mod 11) = (-3) % 11 = 8
验算：7 * 8 % 11 = 1, OK.
扩展欧几里德算法的实现：
def egcd(a: int, b: int) -> tuple[int, int, int]:
if a == 0:
return (b, 0, 1)
else:
g, y, x = egcd(b % a, a)
return (g, x - (b // a) * y, y)
def modinv(a: int, m: int) -> int:
a %= m
g, x, _ = egcd(a, m)
if g != 1:
raise Exception('modular inverse does not exist')
else:
return x % m
使用：
>>> modinv(5, 11)
9
>>> 5 * 9 % 11
1
>>> modinv(7, 11)
8
>>> 7 * 8 % 11
1
其实，在高版本的 Python 里，pow() 接口 已经内置了求逆的算法实现，比如：
>>> pow(5, -1, 11)
9
>>> pow(7, -1, 11)
8
1. 椭圆曲线（Elliptic Curve）
椭圆曲线是指方程为 \(y^2 = x^3 + ax + b\) （其中 a 和 b 为常数）的曲线，其中 \((x, y) \in \mathbb{R}^2\)。
椭圆曲线在平面直角坐标系里的图像，可以在线搜索到，比如维基 Elliptic curve。
1.1. 点加（Point Addition）：
取椭圆曲线上的两个点：\(P(x_1, y_1), Q(x_2, y_2)\)。连接 P、Q 两点的直线（如果 P == Q，取经过该点的切线），与椭圆曲线交于第三点 G，过 G 点作垂直于 X 轴的直线，将过椭圆曲线另一点 R（一般是关于 X 轴对称的点），R 点则被定义为 P+Q 的结果，即 \(P + Q = R\)。
1.2. 计算：
当 P != Q 时，直线的斜率为 \(\lambda = \frac{y_2-y_1}{x_2-x_1}\)；当 P = Q 时，直线（切线）的斜率为：\(\lambda = \frac{3x_1^2 + a}{2 y_1}\)，所以 R 的坐标：
2. 有限域上的椭圆曲线
取一个素数 p，将以上实数域 \(\mathbb{R}\) 改为有限素数域 \(\mathbb{F}_p: {0, 1, 2, \cdots, p-1}\)，等号改为关于模 p 的同余，除法改为模逆运算。
比如取 \(p = 11, a = 0, b = 1\)，就得到了一个椭圆曲线的方程 \(y^2 = x^3 + 1\)，考察 \(\mathbb{F}_{11}\) 里的每一个数，可以知道以下这些点在线上：\((0, 1), (0, 10), (2, 3), (2, 8), (5, 4), (5, 7), (7, 5), (7, 6), (9, 2), (9, 9), (10, 0)\)。
验算如下：比如 \((0, 1)\)：\(x^3+1=0^3+1=1\)，\(y^2=1^2=1\)，符合；再比如 \((5, 7)\)： \(x^3+1=5^3+1=126 \equiv 5 \pmod{11}\)，同时 \(y^2=7^2=49 \equiv 5 \pmod{11}\)，也符合；等等。
2.1. 点加运算的封闭性
跟实数域上一样，有限素域上的点加，也符合运算的封闭性，即点加的结果，仍然在曲线上。比如：
例 1，取两个点 \((0, 1), (2, 3)\)：
所以点加的结果是 \((0, 1) + (2, 3) = (10, 0)\)，它是曲线上的点。
例 2，取同一个点 \((7, 6)\) 做点加：
所以点加的结果是 \((7, 6) + (7, 6) = (2, 3)\)，它是曲线上的点。
例 3，再取 \((2, 3), (7, 6)\) 做加法：
所以点加 \((2, 3) + (7, 6) = (5, 4)\)，也是曲线上的点。
加法表：
添加零元素，构成集合 \(\{\mathbf{0}, (0, 1), (0, 10), (2, 3), (2, 8), (5, 4), (5, 7), (7, 5), (7, 6), (9, 2), (9, 9), (10, 0)\}\)，因为任何两点相加后的结果，仍然在集合里，满足运算的封闭性，可以生成如下的加法表：
++=======++=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+
|| ADD || 0 | (0,1) |(0,10) | (2,3) | (2,8) | (5,4) | (5,7) | (7,5) | (7,6) | (9,2) | (9,9) |(10,0) |
|========++=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+=======+
|| 0 || 0 | (0,1) |(0,10) | (2,3) | (2,8) | (5,4) | (5,7) | (7,5) | (7,6) | (9,2) | (9,9) |(10,0) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (0,1) || (0,1) |(0,10) | 0 |(10,0) | (2,3) | (9,9) | (7,6) | (5,4) | (9,2) | (5,7) | (7,5) | (2,8) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
||(0,10) ||(0,10) | 0 | (0,1) | (2,8) |(10,0) | (7,5) | (9,2) | (9,9) | (5,7) | (7,6) | (5,4) | (2,3) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (2,3) || (2,3) |(10,0) | (2,8) | (0,1) | 0 | (9,2) | (7,5) | (7,6) | (5,4) | (9,9) | (5,7) |(0,10) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (2,8) || (2,8) | (2,3) |(10,0) | 0 |(0,10) | (7,6) | (9,9) | (5,7) | (7,5) | (5,4) | (9,2) | (0,1) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (5,4) || (5,4) | (9,9) | (7,5) | (9,2) | (7,6) |(10,0) | 0 | (2,3) | (0,1) |(0,10) | (2,8) | (5,7) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (5,7) || (5,7) | (7,6) | (9,2) | (7,5) | (9,9) | 0 |(10,0) |(0,10) | (2,8) | (2,3) | (0,1) | (5,4) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (7,5) || (7,5) | (5,4) | (9,9) | (7,6) | (5,7) | (2,3) |(0,10) | (2,8) | 0 | (0,1) |(10,0) | (9,2) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (7,6) || (7,6) | (9,2) | (5,7) | (5,4) | (7,5) | (0,1) | (2,8) | 0 | (2,3) |(10,0) |(0,10) | (9,9) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (9,2) || (9,2) | (5,7) | (7,6) | (9,9) | (5,4) |(0,10) | (2,3) | (0,1) |(10,0) | (2,8) | 0 | (7,5) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
|| (9,9) || (9,9) | (7,5) | (5,4) | (5,7) | (9,2) | (2,8) | (0,1) |(10,0) |(0,10) | 0 | (2,3) | (7,6) |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
||(10,0) ||(10,0) | (2,8) | (2,3) |(0,10) | (0,1) | (5,7) | (5,4) | (9,2) | (9,9) | (7,5) | (7,6) | 0 |
|--------++-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+-------+
2.2. 乘法运算
一个正整数 k 以及一个点 P，定义：
通过加法表，可以整理出如下的乘法表：
P = 0: 0P=0, 1P=0, ...
P = (0,1): 0P=0, 1P=(0,1), 2P=(0,10), 3P=0, ...
P = (0,10): 0P=0, 1P=(0,10), 2P=(0,1), 3P=0, ...
P = (2,3): 0P=0, 1P=(2,3), 2P=(0,1), 3P=(10,0), 4P=(0,10), 5P=(2,8), 6P=0, ...
P = (2,8): 0P=0, 1P=(2,8), 2P=(0,10), 3P=(10,0), 4P=(0,1), 5P=(2,3), 6P=0, ...
P = (5,4): 0P=0, 1P=(5,4), 2P=(10,0), 3P=(5,7), 4P=0, ...
P = (5,7): 0P=0, 1P=(5,7), 2P=(10,0), 3P=(5,4), 4P=0, ...
P = (7,5): 0P=0, 1P=(7,5), 2P=(2,8), 3P=(5,7), 4P=(0,10), 5P=(9,9), 6P=(10,0),
7P=(9,2), 8P=(0,1), 9P=(5,4), 10P=(2,3), 11P=(7,6), 12P=0, ...
P = (7,6): 0P=0, 1P=(7,6), 2P=(2,3), 3P=(5,4), 4P=(0,1), 5P=(9,2), 6P=(10,0),
7P=(9,9), 8P=(0,10), 9P=(5,7), 10P=(2,8), 11P=(7,5), 12P=0, ...
P = (9,2): 0P=0, 1P=(9,2), 2P=(2,8), 3P=(5,4), 4P=(0,10), 5P=(7,6), 6P=(10,0),
7P=(7,5), 8P=(0,1), 9P=(5,7), 10P=(2,3), 11P=(9,9), 12P=0, ...
P = (9,9): 0P=0, 1P=(9,9), 2P=(2,3), 3P=(5,7), 4P=(0,1), 5P=(7,5), 6P=(10,0),
7P=(7,6), 8P=(0,10), 9P=(5,4), 10P=(2,8), 11P=(9,2), 12P=0, ...
P = (10,0): 0P=0, 1P=(10,0), 2P=0, ...
可以看到，对于每一个点 P，\(kP\) 也在曲线上。并且，每一个点 P，倍乘的结果个数是有限的、且循环出现的，比如 P = (2,3)，除了 0，只有 1, 2, 3, 4, 5 这几种倍数，即 \(kP = (k \pmod{6})P\)。再比如 P=(9,2)，只有十二种倍数，等等。
2.3. 计算难题
根据加法和乘法的定义，对于给定的正整数 k 以及一个点 P，可以计算 \(Q = kP\)。即使 k 很大，稍加优化，计算复杂度也很有限。比如 k=1234567890，考虑以下过程：
Q = 0
ks = f'{k:0b}' # k 转换成二进制形式 1001001100101100000001011010010
for i in ks:
Q = Q + Q
if i == '1':
Q += P
这样，经过几十次加法后，就可以得到 Q = kP，并不需要真的进行 k 次加法。
但反过来呢？已知 Q 和 P，求正整数 k，以满足 Q = kP。从目前的结论来看，这个过程并没有实际的快速算法。对于很小的有限素数域，那可以枚举出一条曲线上所有的点，以及任意两个点之间的 “ 倍数 ” 关系，从而通过查表即可得到 k；而对于几百位数的大素数，即使知道 kP 的结果只有有限个，但也远远超过上面举例乘法表的循环范围，所以也没有办法实现所有点的枚举，只能通过实时的反复迭代进行点加计算来求解。
这就是椭圆曲线加密算法采用的难题。这类似于 RSA 算法里的大素数分解难题，即给定两个素数 p 和 q 很容易相乘得到 n，而对 n 进行因式分解却相对困难。
3. 蒙哥马利乘法（Montgomery Multiplication）
在上面 \(\S1.2\) 计算过程中，计算坐标时的运算，特别是求余数时的除法运算，在几百位的大素数时，在计算 \(kP\) 时的反复迭代情况下，计算量是相当大的。
蒙哥马利算法用于降低大数模运算的复杂度，提高计算效率。
3.1. 定义：
其中：
- \(x\) 是有限素域上的数，\(\bar{x}\) 是 \(x\) 在蒙哥马利空间里对应的数；
- \(\otimes\) 是蒙哥马利模乘。
- R：大于 p 的最小的 2 的幂，这样对于与 R 的乘法，可以使用移位操作；
- RR： \(RR \equiv R^2 \pmod{p}\)
经过以下计算：
这样使用蒙哥马利模乘 \(\otimes\) 可以实现 \(x\) 和 \(\bar{x}\) 相互转化，即 \(\bar{x} = x \otimes RR\) 和 \(x = \bar{x} \otimes 1\)。
3.2. 乘法：
利用模乘计算 \(xy \pmod{p}\)：
- 把 x, y 转化为蒙哥马利空间上的数 \(\bar{x}, \bar{y}\)，利用公式：\(\bar{x} = x \otimes RR, \cdots\)
- 进行蒙哥马利模乘 \(\bar{x} \otimes \bar{y}\) 得到 \(\overline{xy}\)，利用公式：\(\bar{x} \otimes \bar{y} = \bar{x} \cdot \bar{y} \cdot R^{-1} = (x \cdot R) \cdot (y \cdot R) \cdot R^{-1} = xy \cdot R = \overline{xy}\)
- 得到结果 xy，利用公式：\(\overline{xy} \otimes 1 = xy\)
一共进行了四次模乘 \(\otimes\)。
再比如利用模乘计算 \(xyz \pmod{p}\)：
- 把 x, y, z 转化为蒙哥马利空间上的数 \(\bar{x}, \bar{y}, \bar{z}\)，利用公式：\(\bar{x} = x \otimes RR, \cdots\)
- 进行蒙哥马利模乘 \(\bar{x} \otimes \bar{y} \otimes \bar{z}\) 得到 \(\overline{xyz}\)，利用公式：\(\bar{x} \otimes \bar{y} \otimes \bar{z} = \overline{xy} \otimes \bar{z} = \overline{xyz}\)
- 得到结果 xyz，利用公式：\(\overline{xyz} \otimes 1 = xyz\)
一共进行了六次模乘 \(\otimes\)。
所以，对于计算 n 个数的模乘 \(x_1 x_2 \cdots x_n\)，也是类似的三个步骤，一共 2n 次模乘。
3.3. 蒙哥马利模乘的实现：
对于指定的椭圆曲线，\(p, a, b\) 都是已知的，所以可以预先计算得到下列参数：
p = 0x12ab655e9a2ca55660b44d1e5c37b00159aa76fed00000010a11800000000001 # 这是一个 253 位的大素数
Rbit = 256
R = 1 << Rbit # 0x10000000000000000000000000000000000000000000000000000000000000000
RR = (1 << (2 * Rbit)) % p # 0x0011fdae7eff1c939a7cc008fe5dc8593cc2c27b58860591f25d577bab861857b
assert p < R
Rmask = R - 1 # 0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
Rinv = modinv(R, p) # 0x07b301912290c02c838557e227b28e2fbd1eae9574fee8754122dd1a1beeec02
mp = (-modinv(p, R)) % R # 0x6992d0fa68b29556249765c34790a000452217cc900000010a117fffffffffff
算法实现：
def redc(t: int) -> int:
a = ((t & Rmask) * mp) & Rmask
u = (t + a * p) >> Rbit
return u if u < p else u - p
def mont_mul(a: int, b: int) -> int:
return redc(a * b)
def mont_add(a: int, b: int) -> int:
t = a + b
return t if t < p else t - p
def mont_to(a: int) -> int:
return mont_mul(a, RR)
def mont_from(a: int) -> int:
return mont_mul(a, 1)
mont_mul 实现了模乘，即 \(\bar{x} \otimes \bar{y}\) 的计算；具体实现在 redc 里。mont_to 和 mont_from 分别实现了 \(x\) 和 \(\bar{x}\) 的相互转换，而它们本身也是基于 mont_mul 模乘，即上面的 \(\bar{x} = x \otimes RR\) 和 \(x = \bar{x} \otimes 1\)。这段算法实现里，用到的计算包括乘法、加减法和位操作，而没有求余数的除法，这就是蒙哥马利算法的优势。
测试代码：
>>> a = 123456789012345678901234567890123456789012345678901234567890
>>> b = 234567890123456789012345678901234567890123456789012345678901
>>>
>>> ma = mont_to(a)
>>> ma
2607102700399835930425256892604740408878764494146331858108218222192440580623
>>>
>>> mb = mont_to(b)
>>> mb
752157449874033112950826507274103158123732489659345112413844821135028041396
>>>
>>> mret = mont_mul(ma, mb)
>>> mret
6837091055679821749041464481328698092587836782355091795225616558442288441760
>>>
>>> ret = mont_from(mret)
>>> ret
2526139555283354312652070910489431265518965079633632635491409854002678665553
>>>
>>> a * b
289589985200426886037189477355433626610714830548588629843609998483265828380439065691255154854446640
31397677651425088890
>>> (a * b) % p
2526139555283354312652070910489431265518965079633632635491409854002678665553
>>>
>>> (a * b) % p == ret
True
4. 参考资料
- Extended Euclidean algorithm https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
- Wikipedia: Elliptic curve https://en.wikipedia.org/wiki/Elliptic_curve
- Wikipedia: Elliptic curve point multiplication https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication
- https://github.com/supranational/blst