Модульные обратные и целые числа без знака - PullRequest
Новая
       69

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

0 голосов
/ 30 ноября 2018

Модульные инверсии могут быть вычислены следующим образом (из Код Розетты ): 

#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-битные комбинации.

Меня особенно интересуют два аспекта:

  • для значений [0, 2 ^ 64) как a, так и b, все вычисления не будут переполнены / переполнены (или переполнены без вреда)?

  • как будет выглядеть (x1 < 0)в неподписанном случае?

1 Ответ

0 голосов
/ 04 декабря 2018

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

...