关于算法:Montgomery算法原理及Python实现

72次阅读

共计 1464 个字符，预计需要花费 4 分钟才能阅读完成。

论文原文：Modular Multiplication Without Trial Division

用处：当 a，b，N 很大时，减速计算a*b(modN)

将 a，b 映射到 Montgomery 域 内：

$$
a’=aR(modN),b’=bR(modN)
$$

此处的 $R=2^k>N$，令 $X=a’b’=abR^2(modN)$

Montgomery约简：

$$
Input:X\\
Output:X/R
$$

约简 1 次：$X=abR^2(modN)$ -> $X/R=abR(modN)$

约简 2 次：$X=abR(modN)$ -> $X/R=ab(modN)$

X/R是不能够间接除的，因为 X 的低 k 位很可能不都为 0。因而心愿有一个数 m 使得

$$
X+mN=0(modR)
$$

由 拓展的欧几里得算法 能够解出以下方程的一组特解：

$$
RR’-NN’=gcd(R,N)
$$

对上述两个公式进行推导，就能够失去 m 的计算方法：

$$
XN’+mNN’=0(modR)\\
XN’+m(RR’-gcd(R,N))=0(modR)\\
XN’=m*gcd(R,N)(modR)\\
m=XN’/gcd(R,N)(modR)
$$

蒙哥马利算法限度 $gcd(R,N)=1$，接着

$$
X+mX=X+XNN'(modR)=X+(XN'(modR))N\\
(X+mX)/R=(X+(XN'(modR))N)/R
$$

以下运算将被大大简化：

$$
aR=a<<k\\
bR=b<<k\\
(X+mN)/R=(X+mN)>>k
$$

1 次蒙哥马利算法中，进行了 2 次对 N 取模，2 次对 R 取模。

申明：上面的代码来自 GitHub，原作者是 LivingLeopold

import math

def get_gcd(a,b):
    if a%b == 0:
        return b
    else :
        return get_gcd(b,a%b)

def get_(a, b):
    if b == 0:
        return 1, 0
    else:
        k = a // b
        remainder = a % b
        x1, y1 = get_(b, remainder)
        x, y = y1, x1 - k * y1
    return x, y

def Multiplicative_inverse_modulo(a , b):
    if b < 0:
        m = abs(b)
    else:
        m = b
    flag = get_gcd(a, b)

    if flag == 1:
        x, y = get_(a, b)
        x0 = x % m
        return x0
    else:
        print("Error!")
        exit()


print('To calculate a * b ( mod N),')
a , b , N = eval(input('Please input a , b , N (Like 68,57,109)：'))

R = 2**(int(math.log2(N)+1))

N1 = Multiplicative_inverse_modulo(N , R)
N2 = R-N1

a_ =(a * R) % N
b_ =(b * R) % N
c_ =(((a_ * b_) + (((a_ * b_) * N2 ) % R ) * N )/ R ) % N

c = int(((c_ + ( c_* N2) % R * N ))/R)

print(str(a) + '*' + str(b) + '(mod' + str(N) +') =' + str(c))

[1] 高效幂模算法探索：Montgomery 算法解析
[2] 蒙哥马利算法详解
[3] 蒙哥马利算法约简、模乘、幂模
[4] 欧几里得算法 & 辗转相除法，百度百科
[5] 拓展的欧几里得算法
[6] RSA 大数运算实现（1024 位 n）（5）蒙哥马利模幂
[7] Python 代码实现蒙哥马利算法

正文完