求解x=a^b(mod m)

http://www.narutoacm.com/archives/a-pow-b-mod-m/

本文致力于解决如下问题：求解x≡ab(mod m)，其中a,b,m都是正整数。

如果b足够小，则可直接用逐次平方法求解，如果你不知道逐次平方法，可以先看这里。所以这里假设b足够大（这不是说是一个64位整数，而是可以上百上千位的一个数），大到逐次平方法也已不足以快速出解。

用素数探路的思想，先假设m是素数，那么要么a与m互素，要么m|a。前者可利用费马小定理，令b = k(m-1) + b′，其中0≤b′<m-1，则有ab≡ak(m-1)+b′≡ab′(mod m)。之后可用逐次平方法快速求解。若m|a，结果显然为0。

现在考虑要求m可以是任意数的情况。同样，若a与m互素，由上边我们可联想到欧拉公式，利用欧拉公式求解。令b = kΦ(m) + b′，其中0≤b′<Φ(m)，则有ab≡akΦ(m)+b′≡ab′(mod m)，之后用逐次平方法快速求解。如果a与m不互素，即gcd(a,m)>1，这种情况下应该怎么做？

注意到b如此大，而模m的不同的数最多只有m个，显而易见的，ab一定和某个很小的指数ab′模m同余，如果找到这个小的指数，就可以利用逐次平方法求解。考虑如下序列：

a0, a1, ..., am(mod m)

由鸽巢原理可知，必有一个最小的r和一个最小的s，使得ar≡ar+s(mod m)，其中r+s≤m。若找到这样的r和s，那么显然有ax≡akr+x0≡akr+s+x0≡ax+s(mod m)，其中x是大于或等于r的任意数。要注意的是，对任意大于或等于r的x，不存在更小的数s1，使得ax≡ax+s1(mod m)。假如有更小的这样的数s1存在，不妨设x = r+ks，因为若不这样，可以令x1 = x + x′ = r + ks，同余式两边同时乘以ax′，就变成ax1≡ax1+s1(mod m)，而x1 = r+ks。所以有ax≡ar+ks≡ar≡ar+s1(mod m)，其中s1<s，这与我们找到的s是最小的相矛盾，所以这样的s1是不存在的。上面的论述给出以下的一个事实：序列

a0, a1, ..., ak, ...(mod m)

实际上是这样一个序列

a0, a1, ..., ar, ar+1, ..., ar+s-1, ar, ar+1, ...(mod m)

即a的0次、1次...幂模m的序列中，前r个数互不相同，从第r+1个数(注：指数为r)开始，每s个数就循环一次。我们把r称为a幂次模m的循环起始点，s称为循环长度。根据以上推导，如果我们能够找到r和s，那么大幂次b就能转换成一个小于m的幂次，然后用逐次平方法就可以求出问题的解。

不妨把a分解成素数乘积的形式：a = p1a1p2a2...pnan，那么a的b次幂模m就转换成每个素数的若干次幂模m的结果相乘再模m，而素数和m的关系只有两种，一种是互素，另一种是m的约数。第一种情况我们已经解决了，所以现在的问题转换成求素数p的b次幂模m，即pb(mod m)，其中p|m。

如前所述，我们可以找到幂模m的r和s，使得pr≡pr+s(mod m)。于是m|(pr+s - pr)，即m|pr(ps - 1)，令m = pr0m′，其中gcd(p,m′)=1。因为gcd(p,ps - 1)=1，所以r至少为r0，即r≥r0，由r的最小性可得r即为r0。于是有m′|(ps - 1)，即有ps≡1(mod m′)。由s的最小性可知，我们只需要找出最小的一个使上式成立的指数，该指数即为s。事实上，这个指数叫做p模m′的次数。

次数定义：对任意正整数a与m，其中gcd(a,m)=1，使得ae≡1(mod m)的最小指数e≥1叫做a模m的次数，记作em(a)。

现在证明如下次数性质：若有an≡1(mod m)，则次数em(a)整除n。特别的，次数总整除Φ(m)。

证明：令g = gcd(n,em(a))，则求得x, y使得nx + em(a)y = g，于是有anx≡1≡ag(mod m)，由次数的最小性可得g = em(a)，即有em(a)|n。又欧拉公式告诉我们aΦ(m)≡1(mod m)，结合上边推理可知次数总整除Φ(m)。证毕！

由次数定义极其性质可知，s即为em′(p)，并且s|Φ(m′)。由于m = pr0m′，且gcd(p,m′)=1，由欧拉Φ函数的积性可得Φ(m′)|Φ(m)，这说明s|Φ(m)。所以有px≡px+Φ(m)(mod m), x≥r=r0，于是pb≡pr+(b-r)(mod Φ(m))(mod m)。由于m = prm′，所以Φ(m)≥Φ(pr)=pr-1(p-1)≥r，其中最后一步可用数学归纳法证明。所以

pb≡pr(mod Φ(m))+Φ(m)+(b-r)(mod Φ(m))≡pΦ(m)+b(mod Φ(m))(mod m)

现在我们可以把a分解成素数乘积，然后对于乘积中每个素数，求出对应幂次模m的值，然后相乘再模m，就得出了解！

能否直接求出a的幂次模m的循环起始点r和循环长度s？！

假设素数p的幂次模m的循环起始点r0和循环长度s0已经求出了，那么pa的幂次模m的循环起始点和循环长度是多少？同样的分析方法，设循环起始点为r，循环长度为s，则有m|pra(psa - 1)，根据前面的论述，我们知道m = pr0m′，其中gcd(p,m′)=1，所以ra≥r0 → r≥r0/a → r≥ceil(r0/a)，由r的定义的最小性可得r = ceil(r0/a)。同样又有m′|(psa - 1)，即psa≡1(mod m′)，所以s = s0/gcd(s0,a)，即有s|s0|Φ(m)。

设数a0的幂次模m的循环起始点为r0，循环长度为s0，数a1的幂次模m的循环起始点为r1，循环长度为s1，其中gcd(a0,a1)=1。现在，我们求a0a1的幂次模m的循环起始点r和循环长度s。由已知可得m|a0r0(a0s0 - 1)，又显然gcd(a0,a0s0-1)=1，因为r0是循环起始点，所以m = a0r0m′，其中gcd(a0,m′)=1。同理也有m = a1r1m″，其中gcd(a1,m″)=1。又由于gcd(a0,a1)=1，所以有m = a0r0a1r1n，其中gcd(a0,n)=gcd(a1,n)=1。由m|(a0a1)r((a0a1)s-1)且gcd(a0a1,(a0a1)s-1)=1且gcd(a0,a1)=1可得，r≥r0且r≥r1。由r的最小性知r即为max(r0,r1)。于是n|((a0a1)s-1)，即(a0a1)s≡1(mod n)，所以s = en(a0a1)。又m′|(a0s0-1)，而n|m′，所以a0s0≡1(mod n)。同理有a1s1≡1(mod n)。所以(a0a1)lcm(s0,s1)≡a0lcm(s0,s1)a1lcm(s0,s1)≡1(mod n)，由次数性质可得s|lcm(s0,s1)。

当把a分解成素数乘积a = p1a1p2a2...pnan时，以上讨论告诉我们a的幂次模m的循环起始点r = max(ceil(ri/ai)),(1≤i≤n)，其中ri是m中包含pi的最大次数。a的幂次模m的循环长度s|lcm(si/gcd(si,ai)),(1≤i≤n) |Φ(m)，其中si为pi的幂次模m的循环长度。因为ri<=Φ(m)，所以r<=Φ(m)。因此有

ab≡ar+(b-r)(mod s)≡ar+(b-r)(mod Φ(m))≡ab(mod Φ(m))+Φ(m)(mod m)

问题解决！（所有公式都没用mathjax显示，看起来可能有点别扭~）