模倒数和无符号整数

Eci*_*ana 5 math integer-overflow greatest-common-divisor modular-arithmetic

模逆可以计算如下(来自Rosetta Code):

#include <stdio.h>

int mul_inv(int a, int b)
{
    int b0 = b, t, q;
    int x0 = 0, x1 = 1;
    if (b == 1) return 1;
    while (a > 1) {
        q = a / b;
        t = b, b = a % b, a = t;
        t = x0, x0 = x1 - q * x0, x1 = t;
    }
    if (x1 < 0) x1 += b0;
    return x1;
}
Run Code Online (Sandbox Code Playgroud)

但是,ints如您所见,输入为 。上面的代码是否也适用于无符号整数(例如uint64_t)?我的意思是,可以用 全部int替换uint64_t吗?我可以尝试少量输入,但尝试所有 64 位组合是不可行的。

我对两个方面特别感兴趣:

  • 对于值 [0, 2^64)ab,所有计算都不会上溢/下溢(或上溢而无害)?

  • 如何将(x1 < 0)看起来像无符号的情况下?

Ser*_*gGr 5

首先这个算法是如何工作的?它基于用于计算GCD的扩展欧几里得算法。简而言之,想法如下:如果我们可以找到一些整数系数,并且使得mn

a*m + b*n = 1
Run Code Online (Sandbox Code Playgroud)

那么m就是模逆问题的答案。很容易看出,因为

a*m + b*n = a*m (mod b)
Run Code Online (Sandbox Code Playgroud)

幸运的是,扩展欧几里得算法正是这样做的:如果ab是互质的,它会找到这样的mn。它的工作方式如下:对于每次迭代,跟踪两个三元组(ai, xai, yai)(bi, xbi, ybi) 这样在每一步

ai = a0*xai + b0*yai
bi = a0*xbi + b0*ybi
Run Code Online (Sandbox Code Playgroud)

ai = 0所以当最终算法停止在和的状态时bi = GCD(a0,b0),那么

1 = GCD(a0,b0) = a0*xbi + b0*ybi
Run Code Online (Sandbox Code Playgroud)

这是使用更明确的方式来计算模来完成的:如果

q = a / b
r = a % b
Run Code Online (Sandbox Code Playgroud)

然后

r = a - q * b
Run Code Online (Sandbox Code Playgroud)

另一件重要的事情是,可以证明对于积极的ab在每一步|xai|,|xbi| <= b|yai|,|ybi| <= a。这意味着在计算这些系数期间不会出现溢出。不幸的是,负值也是可能的,而且,在每个方程中第一个之后的每一步中,一个为正,另一个为负。

您问题中的代码所做的是同一算法的简化版本:由于我们感兴趣的是系数x[a/b],因此它仅跟踪它们并忽略它们y[a/b]。使该代码起作用的最简单方法uint64_t是在单独的字段中显式跟踪符号,如下所示:

a*m + b*n = 1
Run Code Online (Sandbox Code Playgroud)

注意,如果ab不互质或者当b为 0 或 1 时,该问题无解。在所有这些情况下,我的代码返回的0值对于任何真正的解决方案来说都是不可能的。

另请注意,虽然计算值实际上是模逆,但简单的乘法并不总是产生 1,因为与 相乘时会溢出uint64_t。例如a = 688231346938900684b = 2499104367272547425结果是inv = 1080632715106266389

a * inv = 688231346938900684 * 1080632715106266389 = 
= 743725309063827045302080239318310076 = 
= 2499104367272547425 * 297596738576991899 + 1 =
= b * 297596738576991899 + 1
Run Code Online (Sandbox Code Playgroud)

a但如果你对这些和inv类型进行简单的乘法uint64_t,你会得到4042520075082636476这样的结果,而(a*inv)%b不是1543415707810089051预期的结果1