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

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

Ответ

Попробуйте так: 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.