Есть ли способ для сборки, например. (853467 * 21660421200929) % 100000000000007 без библиотек BigInteger (обратите внимание, что каждый номер соответствует 64-битовому целому числу, но результат умножения не соответствует)?
Это решение кажется неэффективным:
int64_t mulmod(int64_t a, int64_t b, int64_t m) {
if (b < a)
std::swap(a, b);
int64_t res = 0;
for (int64_t i = 0; i < a; i++) {
res += b;
res %= m;
}
return res;
}
Вы должны использовать Русское Крестьянское умножение. Он использует повторное удвоение для вычисления всех значений (b*2^i)%m и добавляет их, если установлен i -й бит a.
uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
int64_t res = 0;
while (a != 0) {
if (a & 1) res = (res + b) % m;
a >>= 1;
b = (b << 1) % m;
}
return res;
}
Он улучшает ваш алгоритм, потому что он занимает O(log(a)) время, а не O(a) время.
Предостережения: unsigned и работает только в том случае, если m равно 63 бит или меньше.
Ответ Keith Randall хорош, но, как он сказал, предостережение в том, что оно работает, только если m равно 63 бит.
Вот модификация, которая имеет два преимущества:
m - 64 бит.(Обратите внимание, что строки res -= m и temp_b -= m полагаются на 64-разрядное беззнаковое целочисленное переполнение для получения ожидаемых результатов. Это должно быть хорошо, так как непознанное целочисленное переполнение корректно определено на C и С++. причина важна для использования целочисленных типов без знака.)
uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
uint64_t res = 0;
uint64_t temp_b;
/* Only needed if b may be >= m */
if (b >= m) {
if (m > UINT64_MAX / 2u)
b -= m;
else
b %= m;
}
while (a != 0) {
if (a & 1) {
/* Add b to res, modulo m, without overflow */
if (b >= m - res) /* Equiv to if (res + b >= m), without overflow */
res -= m;
res += b;
}
a >>= 1;
/* Double b, modulo m */
temp_b = b;
if (b >= m - b) /* Equiv to if (2 * b >= m), without overflow */
temp_b -= m;
b += temp_b;
}
return res;
}
Улучшение алгоритма повторного удвоения состоит в том, чтобы проверить, сколько бит сразу можно вычислить без переполнения. Ранняя проверка выхода может быть выполнена для обоих аргументов - ускорение (маловероятное?) Событие N не является простым.
например. 100000000000007 == 0x00005af3107a4007, который позволяет рассчитывать 16 (или 17) бит за каждую итерацию. Фактическое количество итераций будет 3 с примером.
// just a conceptual routine
int get_leading_zeroes(uint64_t n)
{
int a=0;
while ((n & 0x8000000000000000) == 0) { a++; n<<=1; }
return a;
}
uint64_t mulmod(uint64_t a, uint64_t b, uint64_t n)
{
uint64_t result = 0;
int N = get_leading_zeroes(n);
uint64_t mask = (1<<N) - 1;
a %= n;
b %= n; // Make sure all values are originally in the proper range?
// n is not necessarily a prime -- so both a & b can end up being zero
while (a>0 && b>0)
{
result = (result + (b & mask) * a) % n; // no overflow
b>>=N;
a = (a << N) % n;
}
return result;
}
Вы можете попробовать что-то, что разбивает умножение на дополнения:
// compute (a * b) % m:
unsigned int multmod(unsigned int a, unsigned int b, unsigned int m)
{
unsigned int result = 0;
a %= m;
b %= m;
while (b)
{
if (b % 2 != 0)
{
result = (result + a) % m;
}
a = (a * 2) % m;
b /= 2;
}
return result;
}
Оба метода работают для меня. Первый из них такой же, как у вас, но я изменил ваши номера на принудительный ULL. Во-вторых, используется ассемблерная нотация, которая должна работать быстрее. Существуют также алгоритмы, используемые в криптографии (криптография на основе RSA и RSA, в основном, я думаю), как уже упоминалось об уменьшении Montgomery, но я думаю, что для их выполнения потребуется некоторое время.
#include <algorithm>
#include <iostream>
__uint64_t mulmod1(__uint64_t a, __uint64_t b, __uint64_t m) {
if (b < a)
std::swap(a, b);
__uint64_t res = 0;
for (__uint64_t i = 0; i < a; i++) {
res += b;
res %= m;
}
return res;
}
__uint64_t mulmod2(__uint64_t a, __uint64_t b, __uint64_t m) {
__uint64_t r;
__asm__
( "mulq %2\n\t"
"divq %3"
: "=&d" (r), "+%a" (a)
: "rm" (b), "rm" (m)
: "cc"
);
return r;
}
int main() {
using namespace std;
__uint64_t a = 853467ULL;
__uint64_t b = 21660421200929ULL;
__uint64_t c = 100000000000007ULL;
cout << mulmod1(a, b, c) << endl;
cout << mulmod2(a, b, c) << endl;
return 0;
}
Я могу предложить улучшение вашего алгоритма.
Фактически вы вычисляете a * b итеративно, добавляя каждый раз b, делая по модулю после каждой итерации. Лучше добавлять каждый раз b * x, тогда как x определяется так, что b * x не будет переполняться.
int64_t mulmod(int64_t a, int64_t b, int64_t m)
{
a %= m;
b %= m;
int64_t x = 1;
int64_t bx = b;
while (x < a)
{
int64_t bb = bx * 2;
if (bb <= bx)
break; // overflow
x *= 2;
bx = bb;
}
int64_t ans = 0;
for (; x < a; a -= x)
ans = (ans + bx) % m;
return (ans + a*b) % m;
}
a * b % m равно a * b - (a * b/m) * m
Используйте арифметику с плавающей точкой для приближения a * b/m. Приближение оставляет значение достаточно малым для обычных 64-битных целочисленных операций, от m до 63 бит.
Этот метод ограничен значением double, которое обычно составляет 52 бита.
uint64_t mod_mul_52(uint64_t a, uint64_t b, uint64_t m) {
uint64_t c = (double)a * b / m - 1;
uint64_t d = a * b - c * m;
return d % m;
}
Этот метод ограничен значением long double, которое обычно составляет 64 бита или больше. Целочисленная арифметика ограничена 63 битами.
uint64_t mod_mul_63(uint64_t a, uint64_t b, uint64_t m) {
uint64_t c = (long double)a * b / m - 1;
uint64_t d = a * b - c * m;
return d % m;
}
Эти методы требуют, чтобы a и b были меньше, чем m. Чтобы обработать произвольные a и b, добавьте эти строки перед вычислением c.
a = a % m;
b = b % m;
В обоих методах последняя операция % может быть сделана условной.
return d >= m ? d % m : d;