扩展欧几里得算法
费马定理 (Fermat’s theorem )使我们能够通过二进制指数法(binary exponentiation )在O ( log n ) O(\log n) O ( log n ) 的操作内计算模数除法的乘法逆元(modular multiplicative inverses),但这只适用于对质数取模。它的一个推广是欧拉定理(Euler’s theorem ),此定理表明,如果 m 和 a 互质(即他们之间除了1,没有其他公约数),那么
a ϕ ( m ) ≡ 1 ( m o d m ) a^{\phi(m)} ≡ 1(\mod m)
a ϕ ( m ) ≡ 1 ( mod m )
其中,ϕ ( m ) \phi(m) ϕ ( m ) 是欧拉函数 ,定义为小于m m m 的正整数x x x 与m m m 互质的个数。当m m m 是一个质数时,所有小于m m m 的m − 1 m-1 m − 1 个正整数与m m m 都是互质的,所以ϕ ( m ) = m − 1 \phi(m)=m-1 ϕ ( m ) = m − 1 ,这就得出了费马定理。
这允许我们通过ϕ ( m ) \phi (m) ϕ ( m ) 来计算a a a 的模除乘法逆元,记作a ϕ ( m ) − 1 a^{\phi (m)−1} a ϕ ( m ) − 1 ,但反过来,计算ϕ ( m ) \phi (m) ϕ ( m ) 并不那么快速:通常你需要获得m m m 的因数分解 才能做到这一点。有一种更通用的方法,它是通过修改欧几里得算法 (Euclidean algorthm )实现的。
算法
扩展欧几里得算法 (Extended Euclidean algorithm),除了用来寻找最大公约数g = gcd ( a , b ) g = \gcd (a,b) g = g cd( a , b ) ,也可以用来寻找整数x x x 和y y y 满足:
a ⋅ x + b ⋅ y = g a \cdot x + b \cdot y = g
a ⋅ x + b ⋅ y = g
这可以通过将b b b 替换为m m m 和g g g 替换为1 1 1 来解决了寻找模运算逆元的问题:
a − 1 ⋅ a + k ⋅ m = 1 a^{−1}⋅a+k⋅m=1
a − 1 ⋅ a + k ⋅ m = 1
如果a a a 和m m m 不是互质的,那么就没有解,因为a a a 和m m m 的任何整数组合都不能得到一个不是他们最大公约数的倍数的数。
算法同样是递归的:它先计算出gcd ( b , a m o d b ) \gcd (b,a \mod b) g cd( b , a mod b ) 的系数x ′ x' x ′ 和y ′ y' y ′ ,然后恢复初始数字对的解。如果我们有一对( b , a m o d b ) (b,a \mod b) ( b , a mod b ) 的解( x ′ , y ′ ) (x',y') ( x ′ , y ′ ) :
b ⋅ x ′ + ( a m o d b ) ⋅ y ′ = g b⋅x′+(a \mod b)⋅y′=g
b ⋅ x ′ + ( a mod b ) ⋅ y ′ = g
那么,为了得到初始输入的解,我们可以将表达式( a m o d b ) (a \mod b) ( a mod b ) 重写为( a − ⌊ b a ⌋ ⋅ b ) (a−⌊\frac{b}{a}⌋⋅b) ( a − ⌊ a b ⌋ ⋅ b ) 并将它代入前述方程:
b ⋅ x ′ + ( a − ⌊ b a ⌋ ⋅ b ) ⋅ y ′ = g b⋅x′+(a−⌊\frac{b}{a}⌋⋅b)⋅y′=g
b ⋅ x ′ + ( a − ⌊ a b ⌋ ⋅ b ) ⋅ y ′ = g
现在我们按a a a 和b b b 进行分组,重新排列表达式可以得到:
a ⋅ y ′ ⏟ x + b ⋅ ( x ′ − ⌊ b a ⌋ ⋅ y ′ ) ⏟ y = g a \cdot \underbrace{y'}_x + b \cdot \underbrace{(x' - ⌊\frac{b}{a}⌋ \cdot y')}_y = g
a ⋅ x y ′ + b ⋅ y ( x ′ − ⌊ a b ⌋ ⋅ y ′ ) = g
将上式与初始表达式相比较,可以得出a a a 和b b b 的系数即为初始的x x x 和y y y 。
实现
我们将算法实现为递归函数。由于它的输出不是一个而是三个整数,我们通过引用将系数传递给它:
int gcd (int a, int b, int &x, int &y) { if (a == 0 ) { x = 0 ; y = 1 ; return b; } int x1, y1; int d = gcd (b % a, a, x1, y1); x = y1 - (b / a) * x1; y = x1; return d; }
为了计算模运算的乘法逆元,我们只需传递a a a 和m m m ,然后返回算法找到的x x x 系数。由于我们传递了两个正数,其中一个系数将是正的,另一个系数将是负的(具体哪一个取决于迭代次数是奇数还是偶数),所以我们需要检查x x x 是否为负,并加上m m m 以获得正确的余数。
int inverse (int a) { int x, y; gcd (a, M, x, y); if (x < 0 ) x += M; return x; }
这段代码的执行时间约为160ns,比用二进制指数法 求逆元快10ns。为了进一步优化它,我们同样可以用迭代方式实现,这样的执行时间为135ns:
int inverse (int a) { int b = M, x = 1 , y = 0 ; while (a != 1 ) { y -= b / a * x; b %= a; swap (a, b); swap (x, y); } return x < 0 ? x + M : x; }
需要注意的是,和二进制指数法不同,执行时间取决于a a a 的值。例如,对于特定的m m m 值( 1 0 9 + 7 ) (10^9+7) ( 1 0 9 + 7 ) ,性能表现最差的输入值为564400443 564400443 564400443 ,这种情况下,算法执行37次迭代,花费250ns。
练习 :可以尝试将相同的技术应用于二进制GCD (不过,除非你在优化方面比我强,否则这不会获得性能上的提升)。