如何计算高次幂模运算a^b%c的结果？

谢邀，看来有必要回答一下这个问题。

如何计算高次幂模运算a ^ b % c 这个的结果，例如：897987 ^ 23239 % 2323 这样的等式，其实这个问题，有人已经证明而且简化了。

蒙哥马利(Montgomery)幂模运算是快速计算a^b%k的一种算法，是RSA加密算法的核心之一。

模幂运算是RSA 的核心算法，最直接地决定了RSA 算法的性能。针对快速模幂运算这一课题，西方现代数学家提出了大量的解决方案，通常都是先将幂模运算转化为乘模运算。

例如求D=C**15 % N，由于：a*b % n = (a % n)*(b % n) % n，所以：

C1 =C*C % N =C**2 % NC2 =C1*C % N =C**3 % NC3 =C2*C2 % N =C**6 % NC4 =C3*C % N =C**7 % NC5 =C4*C4 % N =C**14 % NC6 =C5*C % N =C**15 % N

即：对于E=15的幂模运算可分解为6 个乘模运算，归纳分析以上方法可以发现对于任意E，都可采用以下算法计算D=C**E % N：

D=1WHILE E>=0IF E%2=0C=C*C % NE=E/2ELSED=D*C % NE=E-1RETURN D

继续分析会发现，要知道E 何时能整除 2，并不需要反复进行减一或除二的操作，只需验证E 的二进制各位是0 还是1 就可以了，从左至右或从右至左验证都可以，从左至右会更简洁，设E=Sum[i=0 to n](E*2**i)，0<=E<=1，则：

D=1FOR i=n TO 0D=D*D % NIF E[i]=1 D=D*C % NRETURN D

这样，模幂运算就转化成了一系列的模乘运算。

模乘运算

对于乘模运算 A*B%N，如果A、B都是1024位的大数，先计算A*B，再% N，就会产生2048位的中间结果，如果不采用动态内存分配技术就必须将大数定义中的数组空间增加一倍，这样会造成大量的浪费，因为在绝大多数情况下不会用到那额外的一倍空间，而采用动态内存分配技术会使大数存储失去连续性而使运算过程中的循环操作变得非常繁琐。所以模乘运算的首要原则就是要避免直接计算A*B。

设A=Sum[i=0 to k](A[i]*r**i)，r=0x10000000，0<=A[i]

可以用一个循环来处理：

C=0;FOR i=n to 0C=C*rC=C+A[i]*B

RETURN C

这样将一个多位乘法转换成了一系列单位乘法和加法，由于：

a*b %n = (a%n * b%n) %na+b %n = (a%n + b%n) %n

所以，对于求C=A*B %N，我们完全可以在求C的过程中完成：

C=0;FOR i=n to 0C=C*r %NC=C+A[i]*B %N RETURN C

这样产生的最大中间结果是A*B 或C*r，都不超过1056(1024+32)位，空间代价会小得多，但是时间代价却加大了，因为求模的过程由一次变成了多次。对于孤立的乘模运算而言这种时间换空间的交易还是值得的，但是对于反复循环的乘模运算，这种代价就无法承受，必须另寻出路。

蒙哥马利模乘要解决的问题：

｛注：蒙哥马利模乘实际上是解决了这样一个问题，即不使用除法（用移位操作）而求得模乘运算的结果。时刻注意是模运算而且是大数运算的基础上的。｝

｛例如：假设进制 R=10 一个数(大数表示)23 =2*10^1+3*10^0 。欲求23 mod 5，先求 23 *10**-K mod 5的值 ，不使用乘法，我们可以采取下面的办法 就是将 23+5*q 这时不影响模运算的结果 当23+5*q 是10的倍数时 就可以用移位操作除以10，一般k的值取23（大数）的位数（在进制R的基础上）一直移位，最后剩下一个大于模5 小于2*5的数，在减去一个5 就是最后的结果了，这样求出的是23 * 10^-k mod 5 的结果。我们想求23 mod 5的结果只需要先求出23*10^-k mod 5 = Z即可 在去求 z*10^k mod 5 即可 ，因为 23 % 5 = (23 * 10**K * 10 ** (-K))%5 = ((23 * 10**(-K)%5 *10 ** (-K)%5)%5。

这个例子举得不好，因为23+5q 永远也加不出10的倍数，但是由于RSA算法的特殊性，在蒙哥马利模乘时是可以加出进制倍数，而实现移位的，这个例子只是领会精神，也可以不看｝

·蒙哥马利约减表达式： Mon =（S+qM）/R = （S+qM）*1/R （S为 其中S表示被归约数，M表示模数，R = 2^n（二进制情况下），n表示指定0的个数。）。S+qM实际表示模M的S所在的剩余类的所有数

·看做剩余类的运算：

M(S)*M(1) = M(S)*M(R)*M(R^(-1))

所以Mon mod M = S*R^(-1) mod M （所以想求S*R^(-1) mod M的值可以用Mon mod M来计算，找到q使得S+qM是R的倍数，这样被除数就是整数了）

·Montgomery乘法：Z = X*Y*R^(-1) mod M

·(R与M要互素，当R是2^n时，M是奇数即可)｝

蒙哥马利模乘器可用硬件实现，在配合软件实现RSA运算，以加快运算速度。

蒙哥马利模乘(下列表述不是很容易理解，可以参考上面的注)

由于RSA 的核心算法是模幂运算，模幂运算又相当于模乘运算的循环，要提高RSA 算法的效率，首要问题在于提高模乘运算的效率。不难发现，模乘过程中复杂度最高的环节是求模运算，因为一次除法实际上包含了多次加法、减法和乘法，如果在算法中能够尽量减少除法甚至避免除法，则算法的效率会大大提高。

设A=Sum[i=0 to k](A[i]*2**i)，0<=A[i]<=1，则：C= A*B = Sum[i=0 to k](A*B*2**i) 可用循环处理为：

C=0FOR i FROM k TO 0C=C*2C=C+A[i]*B RETURN C

若令 C'= A*B *2**(-k)则

C'= Sum[i=0 to k](A[i]*B*2**(i-k))

用循环处理即：

C'=0FOR i FROM 0 TO kC'=C'+A[i]*BC'=C'/2RETURN C'

通过这一算法求A*B*2**(-k)是不精确的，因为在循环中每次除以2都可能有余数被舍弃了，但是可以通过这一算法求A*B*2**(-k) %N的精确值，方法是在对C'除2之前，让C'加上C'[0]*N。由于在RSA中N是两个素数的积，总是奇数，所以当C'是奇数时，C'[0]=1，C'+C'[0]*N 就是偶数，而当C'为偶数时C'[0]=0，C'+C'[0]*N还是偶数，这样C'/2 就不会有余数被舍弃。又因为C'+N %N = C' %N，所以在计算过程中加若干次N，并不会影响结果的正确性。可以将算法整理如下：

{A*B*2**(-k)%N 的含义是找到一个与A*B%N 同余的一个数H*2**K，再计算H%N}

C'=0FOR i FROM 0 TO kC'=C'+A*BC'=C'+C'[0]*NC'=C'/2IF C'>=N C'=C'-NRETURN C'

由于在RSA中A、B总是小于N，又0<=A,C'[0]<=1，所以：

C' = (C'+A*B+C'[0]*N)/2C' < (C'+2N)/22C' < C'+2NC' < 2N

既然C'总是小于2N，所以求C' %N 就可以很简单地在结束循环后用一次减法来完成，即在求A*B*2**(-k) %N的过程中不用反复求模，达到了我们避免做除法的目的。当然，这一算法求得的是A*B*2**(-k) %N，而不是我们最初需要的A*B %N。但是利用A*B*2**(-k)我们同样可以求得A**E %N。

设R=2**k %N，R'=2**(-k) %N，E=Sum[i=0 to n](E[i]*2**i)：

A'=A*R %N //这一步是怎么求的？X=A'FOR i FROM n TO 0X=X*X*R' %NIF E[i]=1 X=X*A'*R' %NX=X*1*R' %NRETURN X

最初：

X = A*R %N，

开始循环时：

X = X*X*R' %N= A*R*A*R*R' %N = A**2*R %N

反复循环之后：

X = A**E*R %N

最后：

X = X*1*R' %N= A**E*R*R' %N= A**E %N

如此，我们最终实现了不含除法的模幂算法，这就是著名的蒙哥马利算法，而X*Y*R' %N 则被称为“蒙哥马利模乘”。以上讨论的是蒙哥马利模乘最简单，最容易理解的二进制形式。蒙哥马利算法的核心思想在于将求A*B %N转化为不需要反复取模的A*B*R' %N，（移位即可，因为R是2^K,总之R是与进制相关的数），但是利用二进制算法求1024位的A*B*R' %N，需要循环1024次之多，我么必然希望找到更有效的计算A*B*R' %N的算法。

{A*B%N

=(A*B*2**(-K) * 2**K )%N

=(A*B*2**(-K)%N * 2**K%N)%N

=((H*2**K * 2**(-K))%N * 2**K%N)%N

=(H%N * 2**K%N)%N //这里把A*B*2**(-K) 通过移位转化成了 H%N 而H是一个小于2N的数，做H-N即可。

}

考虑将A表示为任意的r进制：

A = Sum[i=0 to k](A*r**i) 0<=A<=r

我们需要得到的蒙哥马利乘积为：

C'= A*B*R' %N R'=r**(-k)

则以下算法只能得到C'的近似值

C'=0FOR i FROM 0 TO kC'=C'+A*BC'=C'/rIF C'>=N C'=C'-NRETURN C'

因为在循环中每次C'=C'/r 时，都可能有余数被舍弃。假如我们能够找到一个系数 q，使得(C' + A*B + q*N) %r =0，并将算法修改为：

C'=0FOR i FROM 0 TO kC'=C'+A*B+q*NC'=C'/rIF C'>=N C'=C'-NRETURN C'

则C'的最终返回值就是A*B*R' %N的精确值，所以关键在于求q。由于：

(C' + A*B + q*N) %r =0==> (C' %r + A*B %r + q*N %r) %r =0==> (C'[0] + A*B[0] + q*N[0]) %r =0

若令N[0]*N[0]' %r =1，q=(C'[0]+A*B[0])*(r-N[0]') %r，则：

(C'[0] + A*B[0] + q*N[0]) %r= (C'[0]+A*B[0] - (C'[0]+A*B[0])*N[0]'*N[0]) %r) %r= 0

于是我们可以得出r为任何值的蒙哥马利算法：

m=r-N[0]'C'=0FOR i FROM 0 TO kq=(C'[0]+A*B[0])*m %rC'=(C'+A*B+q*N)/rIF C'>=N C'=C'-NRETURN C'

如果令 r=0x100000000，则 %r 和 /r 运算都会变得非常容易，在1024位的运算中，循环次数k 不大于32，整个运算过程中最大的中间变量C'=(C'+A*B+q*N)< 2*r*N < 1057位，算法效率就相当高了。唯一的额外负担是需要计算 N[0]'，使N[0]*N[0]' %r =1，而这一问题前面已经用欧几里德算法解决过了，而且在模幂运算转化成反复模乘运算时，N是固定值，所以N[0]'只需要计算一次，负担并不大。

素数测试方法

数论学家利用费马小定理研究出了多种素数测试方法，目前最快的算法是拉宾米勒测试算法，其过程如下：

（1）计算奇数M，使得N=(2**r)*M+1（2）选择随机数A

若N 通过一次测试，则N 不是素数的概率为 25%，若N 通过t 次测试，则N 不是素数的概率为1/4**t。事实上取t 为5 时，N 不是素数的概率为 1/128，N 为素数的概率已经大于99.99%。

在实际应用中，可首先用300—500个小素数对N 进行测试，以提高拉宾米勒测试通过的概率，从而提高整体测试速度。而在生成随机素数时，选取的随机数最好让r=0，则可省去步骤（3） 的测试，进一步提高测试速度。

素数测试是RSA 选取秘钥的第一步，奇妙的是其核心运算与加解密时所需的运算完全一致：都是模幂运算。而模幂运算过程中中所需求解的欧几里德方程又恰恰正是选取密钥第二步所用的运算。可见整个RSA 算法具有一种整体的和谐。