如何计算\(x\cdot y \mod N\)？
传统的模乘运算
在\(Z_{1024}\)中，将其中一个数\(x\)表示成\(x = \displaystyle\sum_{i=0}^{1023} a_i \cdot 2^i\)，于是，乘法运算可以写为
\[x\cdot y \mod N = (\displaystyle\sum_{i=0}^{1023}a_i \cdot 2^i)\mod N\\
\]
Input: x, y, N
Output: x*y mod N
z = 0
for i = 1023 downto 0 step -1
z *= 2
z += a[i] * y
endfor
return z % N
上述运算过程中，\(z\times 2\)可能会导致值超过\(1024\)比特，所以需要将模运算移到循环体中，修改后算法如下
Input: x, y, N
Output: x*y mod N
z = 0
for i = 1023 downto 0 step -1
z *= 2
z %= N // z = (z>=N) ? z-N : z
c += z[i] * y
z %= N // z = (z>=N) ? z-N : z
endfor
return z % N
蒙哥马利算法
用\(r\)进制表示乘数\(x\)，设\(k\)为\(N\)以\(r\)进制表示时的位数，我们考虑这样一种算法
\[x \cdot y \cdot r^{-k} \mod N\\
\]
其中
\[x = \displaystyle\sum_{i=0}^{k-1} a_i \cdot r^{i}\\
\]
这样原式可以表示为
\[x\cdot y\cdot r^{-k} \mod N = (\displaystyle\sum_{i=0}^{k-1} a_i \cdot r^{i} \cdot y \cdot r^{-k})\mod N\\
\]
Input: x,y,N
Output: x*y*r^{-k} mod N
z = 0
for i = 1023 downto 0 step -1
z *= r
z += a[i] * y
z /= r
endfor
return z % N
但是\(z\)可能不是\(r\)的倍数，\(\frac{z}{r}\)相当于\(z\)右移，可能失去最右边的\(1\)，因此需要给\(z\)加上一个值，使\(z\)可以被\(r\)整除，为了不影响取模的结果，这个值应当是\(N\)的整数倍。
Input: x,y,N
Output: x*y*r^{-k} mod N
z = 0
for i = 1023 downto 0 step -1
z *= r
z += a[i] * y + λ*N
z /= r
endfor
return z % N
其中\((z + a[i] \cdot y + \lambda \cdot N) \mod r = 0\)，得到\(\lambda = (z[0] + a[i]\cdot y[0]) * (-N[0]^{-1})\mod N\)
Input: x,y,N
Output: x*y*r^{-k} mod N
z = 0
z0 = 0
y0 = y % r
N0 = N % r
w = -N0^{-1} % r
for i = k-1 downto 0 step -1
z *= r
λ = ((z0 + a[i] * y0) * w) % r
z += a[i] * y + λ*N // 不会损失最低位
z /= r
z0 = z % r
endfor
return (z >= N) ? z - N : z
算法中的\(y[0],N[0],w\)都可以预计算得到，\(\% r\)运算可将转化为移位运算实现。这样就通过蒙哥马利算法得到了\(x\cdot y\cdot r^{-k}\mod N\)，接着
\[\begin{aligned}
a' &= x \cdot r^{2k} \cdot r^{-k} \mod N\\
b' &= y \cdot r^{2k} \cdot r^{-k} \mod N\\
z' &= a' \cdot b' \cdot r^{-k} \mod N\\
z &= z' \cdot 1 \cdot r^{-k} \mod N
\end{aligned}
\]
就可以得到\(x \cdot y \mod N\)
相较于传统的模乘运算方法，蒙哥马利算法将模运算转化为加减运算和移位运算，提高了模运算的执行效率。