Модульные инверсии и целые числа без знака

Прежде всего, как работает этот алгоритм? Он основан на расширенном алгоритме Евклида для вычисления 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.