Произведение чисел порядка 10^9 по модулю

#c #алгоритм #cpp

Столкнулся с задачей вычисления произведения двух разных чисел a,b по модулю n(не
обязательно простому):
10^9 < a,b <= n.
Собственно, простой вариант:
c=(a*b)%n;
(еще раз замечу, что всегда a,b <= n)
не подходит, ибо, при a,b > 4.29496729610^9, произведение ab >2^64.
Подскажите, пожалуйста, как реализовать подобную операцию.
Естественно, чем быстрее она будет выполняться, тем лучше.
Заранее спасибо всем откликнувшимся!    

Ответы

Ответ 1

Попробуйте так:
static int MulPeasant(int a, int b, int n)
{
    a = a % n; // не нужно, если гарантированно a < n
    b = b % n; // не нужно, если гарантированно b < n
    int acc = 0;
    while (b > 0)
    {
        if ((b & 1) != 0)
            acc = (acc + a) % n;
        a = (a * 2) % n;
        b >>= 1;
    }
    return acc;
}

Выражение (acc + a * b) % n — инвариант цикла.
Должно работать, если n < MAXINT/2.
(Если я не ошибаюсь, это так называемый «русский крестьянский метод»).

Для случая произвольного n подойдёт небольшая модификация:
static int MulPeasant(int a, int b, int n)
{
    a = a % n; // не нужно, если гарантированно a < n
    b = b % n; // не нужно, если гарантированно b < n
    int acc = 0;
    while (b > 0)
    {
        if ((b & 1) != 0)
            acc = SafeAdd(acc, a, n);
        a = SafeAdd(a, a, n);
        b >>= 1;
    }
    return acc;
}

static int SafeAdd(int x, int y, int n)
{
    if (x <= MAXINT - y) // проверка на переполнение сложения
    {
        int result = x + y; // x, y < n, x + y < 2n. (x + y) mod n равно или x + y,
        return (result >= n) ? (result - n) : result;  // или x + y - n
    }
    else // в этом случае x + y > MAXINT >= n, но x + y < n + n = 2n
    {    // значит (x + y) mod n = x + y - n
        return x - (n - y); // n - y > 0, вычитаем два положительных числа меньших n
    }
}


Можно ещё организовать поразрядное умножение:
const int halfbits = sizeof(int) * 8 / 2;
const int halfbase = (1 << halfbits);
const int lomask = halfbase - 1;

static int MulDigits(int a, int b, int n)
{
    int reducedbase = halfbase % n;
    int reducedbasesquare = (reducedbase * reducedbase) % n;
    int ah = a >> halfbits, al = a & lomask;
    int bh = b >> halfbits, bl = b & lomask;

    int ll = (al * bl) % n,
        lh = (((ah * bl) % n) * reducedbase) % n,
        hl = (((al * bh) % n) * reducedbase) % n,
        hh = (((ah * bh) % n) * reducedbasesquare) % n;
    return SafeAdd(
               SafeAdd(ll, lh, n),
               SafeAdd(hl, hh, n),
               n);
}

(Если ваше n одно и то же всё время, константы reducedbase и reducedbasesquare можно
предвычислять.)

А вообще, почитайте «Искусство программирования», главы 4.3.2 и 4.3.3.