math - 模倒数和无符号整数

Question

模逆可以计算如下(来自 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;
}

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

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

对于 a 的值 [0, 2^64)和 b , 是否所有计算都不会溢出/下溢(或溢出而无害)？

怎么会(x1 < 0)看起来像在未签名的情况下？

Answer 1

首先这个算法是如何工作的？它基于 Extended Euclidean algorithm用于计算 GCD .简而言之，想法如下:如果我们能找到一些整数系数 m和 n以至于

a*m + b*n = 1

然后 m将是模逆问题的答案。很容易看到，因为

a*m + b*n = a*m (mod b)

幸运的是，扩展欧几里得算法正是这样做的:如果 a和 b是互质的，它找到这样的 m和 n .它的工作方式如下:对于每个迭代跟踪两个三元组 (ai, xai, yai)和 (bi, xbi, ybi)以至于在每一步

ai = a0*xai + b0*yai
bi = a0*xbi + b0*ybi

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

1 = GCD(a0,b0) = a0*xbi + b0*ybi

它使用更明确的方式来计算模数:如果

q = a / b
r = a % b

然后

r = a - q * b

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

您问题中的代码所做的是同一算法的简化版本:因为我们感兴趣的是 x[a/b]系数，它只跟踪它们并忽略 y[a/b]那些。使该代码适用于 uint64_t 的最简单方法是在单独的字段中明确跟踪标志，如下所示:

typedef struct tag_uint64AndSign {
    uint64_t  value;
    bool isNegative;
} uint64AndSign;


uint64_t mul_inv(uint64_t a, uint64_t b)
{
    if (b <= 1)
        return 0;

    uint64_t b0 = b;
    uint64AndSign x0 = { 0, false }; // b = 1*b + 0*a
    uint64AndSign x1 = { 1, false }; // a = 0*b + 1*a

    while (a > 1)
    {
        if (b == 0) // means original A and B were not co-prime so there is no answer
            return 0;
        uint64_t q = a / b;
        // (b, a) := (a % b, b)
        // which is the same as
        // (b, a) := (a - q * b, b)
        uint64_t t = b; b = a % b; a = t;

        // (x0, x1) := (x1 - q * x0, x0)
        uint64AndSign t2 = x0;
        uint64_t qx0 = q * x0.value;
        if (x0.isNegative != x1.isNegative)
        {
            x0.value = x1.value + qx0;
            x0.isNegative = x1.isNegative;
        }
        else
        {
            x0.value = (x1.value > qx0) ? x1.value - qx0 : qx0 - x1.value;
            x0.isNegative = (x1.value > qx0) ? x1.isNegative : !x0.isNegative;
        }
        x1 = t2;
    }

    return x1.isNegative ? (b0 - x1.value) : x1.value;
}

请注意，如果 a和 b不互质或当 b是 0 或 1，这个问题无解。在所有这些情况下，我的代码返回 0对于任何真正的解决方案来说，这是一个不可能的值。

还要注意，虽然计算出的值实际上是模逆，但简单的乘法并不总是产生 1，因为乘法超过 uint64_t 时会溢出。 .例如 a = 688231346938900684和 b = 2499104367272547425结果是 inv = 1080632715106266389

a * inv = 688231346938900684 * 1080632715106266389 = 
= 743725309063827045302080239318310076 = 
= 2499104367272547425 * 297596738576991899 + 1 =
= b * 297596738576991899 + 1

但是如果你对这些 a 做一个简单的乘法运算和 inv类型 uint64_t , 你会得到 4042520075082636476所以 (a*inv)%b将是 1543415707810089051而不是预期 1 .