Самый точный способ сделать операцию умножения и деления в 64-битной?
Как наиболее точно можно выполнить операцию умножения и деления для 64-разрядных целых чисел, которая работает как в 32-разрядных, так и в 64-разрядных программах (в Visual C++)? (В случае переполнения мне нужен результат mod 264.)
(Я ищу что-то вроде MulDiv64, за исключением того, что он использует встроенную сборку, которая работает только в 32-битных программах.)
Очевидно, приведение к double и обратно возможно, но мне интересно, есть ли более точный способ, который не слишком сложен. (т.е. я не ищу здесь арифметические библиотеки произвольной точности!)
12 ответов
Так как это помечено Visual C++, я дам решение, которое использует MSVC-специфичные свойства.
Этот пример довольно сложный. Это очень упрощенная версия того же алгоритма, который используется GMP и java.math.BigInteger для большого подразделения.
Хотя я имею в виду более простой алгоритм, он, вероятно, примерно в 30 раз медленнее.
Это решение имеет следующие ограничения / поведение:
- Требуется х64. Он не скомпилируется на x86.
- Коэффициент не равен нулю.
- Коэффициент насыщается, если он переполнен 64-битными.
Обратите внимание, что это для целого числа без знака. Тривиально создать оболочку для этого, чтобы она работала и для подписанных дел. Этот пример также должен давать правильно усеченные результаты.
Этот код не полностью протестирован. Тем не менее, он прошел все тесты, которые я бросил в него.
(Даже случаи, которые я специально построил, чтобы попытаться сломать алгоритм.)
#include <intrin.h>
uint64_t muldiv2(uint64_t a, uint64_t b, uint64_t c){
// Normalize divisor
unsigned long shift;
_BitScanReverse64(&shift,c);
shift = 63 - shift;
c <<= shift;
// Multiply
a = _umul128(a,b,&b);
if (((b << shift) >> shift) != b){
cout << "Overflow" << endl;
return 0xffffffffffffffff;
}
b = __shiftleft128(a,b,shift);
a <<= shift;
uint32_t div;
uint32_t q0,q1;
uint64_t t0,t1;
// 1st Reduction
div = (uint32_t)(c >> 32);
t0 = b / div;
if (t0 > 0xffffffff)
t0 = 0xffffffff;
q1 = (uint32_t)t0;
while (1){
t0 = _umul128(c,(uint64_t)q1 << 32,&t1);
if (t1 < b || (t1 == b && t0 <= a))
break;
q1--;
// cout << "correction 0" << endl;
}
b -= t1;
if (t0 > a) b--;
a -= t0;
if (b > 0xffffffff){
cout << "Overflow" << endl;
return 0xffffffffffffffff;
}
// 2nd reduction
t0 = ((b << 32) | (a >> 32)) / div;
if (t0 > 0xffffffff)
t0 = 0xffffffff;
q0 = (uint32_t)t0;
while (1){
t0 = _umul128(c,q0,&t1);
if (t1 < b || (t1 == b && t0 <= a))
break;
q0--;
// cout << "correction 1" << endl;
}
// // (a - t0) gives the modulus.
// a -= t0;
return ((uint64_t)q1 << 32) | q0;
}
Обратите внимание, что если вам не нужен идеально усеченный результат, вы можете полностью удалить последний цикл. Если вы сделаете это, ответ будет не более чем на 2 больше, чем правильный коэффициент.
Тестовые случаи:
cout << muldiv2(4984198405165151231,6132198419878046132,9156498145135109843) << endl;
cout << muldiv2(11540173641653250113, 10150593219136339683, 13592284235543989460) << endl;
cout << muldiv2(449033535071450778, 3155170653582908051, 4945421831474875872) << endl;
cout << muldiv2(303601908757, 829267376026, 659820219978) << endl;
cout << muldiv2(449033535071450778, 829267376026, 659820219978) << endl;
cout << muldiv2(1234568, 829267376026, 1) << endl;
cout << muldiv2(6991754535226557229, 7798003721120799096, 4923601287520449332) << endl;
cout << muldiv2(9223372036854775808, 2147483648, 18446744073709551615) << endl;
cout << muldiv2(9223372032559808512, 9223372036854775807, 9223372036854775807) << endl;
cout << muldiv2(9223372032559808512, 9223372036854775807, 12) << endl;
cout << muldiv2(18446744073709551615, 18446744073709551615, 9223372036854775808) << endl;
Выход:
3337967539561099935
8618095846487663363
286482625873293138
381569328444
564348969767547451
1023786965885666768
11073546515850664288
1073741824
9223372032559808512
Overflow
18446744073709551615
Overflow
18446744073709551615
Вам просто нужно 64-битные целые числа. Есть несколько избыточных операций, но это позволяет использовать 10 в качестве базы и шага в отладчике.
uint64_t const base = 1ULL<<32;
uint64_t const maxdiv = (base-1)*base + (base-1);
uint64_t multdiv(uint64_t a, uint64_t b, uint64_t c)
{
// First get the easy thing
uint64_t res = (a/c) * b + (a%c) * (b/c);
a %= c;
b %= c;
// Are we done?
if (a == 0 || b == 0)
return res;
// Is it easy to compute what remain to be added?
if (c < base)
return res + (a*b/c);
// Now 0 < a < c, 0 < b < c, c >= 1ULL
// Normalize
uint64_t norm = maxdiv/c;
c *= norm;
a *= norm;
// split into 2 digits
uint64_t ah = a / base, al = a % base;
uint64_t bh = b / base, bl = b % base;
uint64_t ch = c / base, cl = c % base;
// compute the product
uint64_t p0 = al*bl;
uint64_t p1 = p0 / base + al*bh;
p0 %= base;
uint64_t p2 = p1 / base + ah*bh;
p1 = (p1 % base) + ah * bl;
p2 += p1 / base;
p1 %= base;
// p2 holds 2 digits, p1 and p0 one
// first digit is easy, not null only in case of overflow
uint64_t q2 = p2 / c;
p2 = p2 % c;
// second digit, estimate
uint64_t q1 = p2 / ch;
// and now adjust
uint64_t rhat = p2 % ch;
// the loop can be unrolled, it will be executed at most twice for
// even bases -- three times for odd one -- due to the normalisation above
while (q1 >= base || (rhat < base && q1*cl > rhat*base+p1)) {
q1--;
rhat += ch;
}
// subtract
p1 = ((p2 % base) * base + p1) - q1 * cl;
p2 = (p2 / base * base + p1 / base) - q1 * ch;
p1 = p1 % base + (p2 % base) * base;
// now p1 hold 2 digits, p0 one and p2 is to be ignored
uint64_t q0 = p1 / ch;
rhat = p1 % ch;
while (q0 >= base || (rhat < base && q0*cl > rhat*base+p0)) {
q0--;
rhat += ch;
}
// we don't need to do the subtraction (needed only to get the remainder,
// in which case we have to divide it by norm)
return res + q0 + q1 * base; // + q2 *base*base
}
Это ответ сообщества вики, так как на самом деле это просто набор ссылок на другие статьи / ссылки (я не могу опубликовать соответствующий код).
Умножение двух 64-битных целочисленных значений на 128-битный результат довольно просто, если использовать простую технику карандаша и бумаги, которую каждый изучает в начальной школе.
Комментарий Грегса верен: Кнут описывает раздел "Искусство компьютерного программирования", второе издание, том 2 / Получисленные алгоритмы "в конце Раздела 4.3.1 Арифметика с множественной точностью / Классические алгоритмы (страницы 255 - 265 в моем экземпляре). Это не легко читать, по крайней мере, не для кого-то вроде меня, который забыл большинство математики за 7-й класс алгебры. Как раз перед этим Кнут рассматривает и сторону умножения.
Некоторые другие варианты для идей (эти примечания для алгоритмов деления, но большинство также обсуждают умножение):
- Джек Креншоу освещает алгоритмы деления Кнута более читабельно в серии статей из журнала Embedded System Programming 1997 (к сожалению, в моих заметках нет точных проблем). К сожалению, статьи из старых проблем ESP не так легко найти в Интернете. Если у вас есть доступ к университетской библиотеке, возможно, у вас есть какие-то проблемы или копия библиотеки ESP CD-ROM.
- Томас Родехеффер из Microsoft Research опубликовал статью о Software Integer Division: http://research.microsoft.com/pubs/70645/tr-2008-141.pdf
- Статья Карла Хасселстрёма "Быстрое деление больших целых чисел": http://www.treskal.com/kalle/exjobb/original-report.pdf
- "Язык ассемблера" Рэндалла Хайда (http://webster.cs.ucr.edu/AoA/Windows/HTML/AoATOC.html), в частности, том четвертый, раздел 4.2.5 (Расширенное прецизионное деление): http://webster.cs.ucr.edu/AoA/Windows/HTML/AdvancedArithmetica2.html это в варианте Hyde для языка ассемблера x86, но есть также некоторый псевдокод и достаточное объяснение, чтобы перенести алгоритм на C. Это тоже медленно - выполнение деление по крупицам...
Поскольку это помечено как Visual C++, вы можете использовать недавно доступные встроенные функции:
- Мул:
_mul128(),_umul128() - Div:
_div128(),_udiv128()
uint64_t muldiv_u64(uint64_t a, uint64_t b, uint64_t c)
{
uint64_t highProduct;
uint64_t lowProduct = _umul128(a, b, &highProduct);
uint64_t remainder;
return _udiv128(highProduct, lowProduct, c, &remainder);
}
Если вам нужен подписанный мультидив, просто используйте версию без
u
Смотрите также
Вам не нужна произвольная арифметика точности для этого. Вам нужна только 128-битная арифметика. Т.е. вам нужно умножение 64*64=128 и деление 128/64=64 (с правильным поведением переполнения). Это не так сложно реализовать вручную.
Ниже представлены исходники muldiv64.asm для Windows-x86-64 Visual Studio. Более подробная информация о компиляции/тестировании/исключениях доступна на https://github.com/LF994/MulDiv64 .
include listing.inc
_TEXT SEGMENT
PUBLIC asmMulDiv64
asmMulDiv64 PROC
mov rax,rdx ; rax <- b
mul rcx ; 128 bits of rax*rcx placed into rdx+rax (a was in rcx)
div r8 ; [rdx,rax] / r8: result->rax, reminder->rdx (c was in r8)
; Quotient already in rax, place here "mov rax,rdx" if you want to return remainder instead.
ret 0 ; 0: do not remove parameters from stack
asmMulDiv64 ENDP
_TEXT ENDS
END
Ну, вы можете разделить свои 64-битные операнды на 32-битные порции (низкие и высокие части). Чем делать операции на них, что вы хотите. Все промежуточные результаты будут менее 64 бит и поэтому могут быть сохранены в ваших типах данных.
Если вам нужно только поддерживать Windows 7 и новее, один хороший способ заключается в следующем:
#include <mfapi.h>
#include <assert.h>
#pragma comment( lib, "mfplat.lib" )
uint64_t mulDiv64( uint64_t a, uint64_t b, uint64_t c )
{
assert( a <= LLONG_MAX && b <= LLONG_MAX && c <= LLONG_MAX );
// https://docs.microsoft.com/en-us/windows/desktop/api/Mfapi/nf-mfapi-mfllmuldiv
return (uint64_t)MFllMulDiv( (__int64)a, (__int64)b, (__int64)c, (__int64)c / 2 );
}
Этот метод намного проще, чем в других ответах, он округляет результат вместо усечения и работает на всех платформах Windows, включая ARM.
Есть ли в вашем распоряжении тип COMP (64-битный целочисленный тип на основе x87) в VC++? Я использовал это время от времени в Delphi, когда мне была нужна 64-битная целочисленная математика. В течение многих лет это было намного быстрее, чем основанная на библиотеке 64-битная целочисленная математика - конечно, когда деление было вовлечено.
В Delphi 2007 (последняя версия, которую я установил - 32-битная версия) я бы реализовал MulDiv64 следующим образом:
function MulDiv64(const a1, a2, a3: int64): int64;
var
c1: comp absolute a1;
c2: comp absolute a2;
c3: comp absolute a3;
r: comp absolute result;
begin
r := c1*c2/c3;
end;
(Эти странные абсолютные операторы накладывают переменные comp поверх их 64-битных целочисленных счетчиков. Я бы использовал простые приведения типов, за исключением того, что компилятор Delphi не запутался в этом - возможно, потому что язык Delphi (или как они его сейчас называют) имеет нет четкого синтаксического различия между приведением типа (переинтерпретация) и преобразованием типа значения.)
В любом случае, Delphi 2007 отображает вышесказанное следующим образом:
0046129C 55 push ebp
0046129D 8BEC mov ebp,esp
0046129F 83C4F8 add esp,-$08
004612A2 DF6D18 fild qword ptr [ebp+$18]
004612A5 DF6D10 fild qword ptr [ebp+$10]
004612A8 DEC9 fmulp st(1)
004612AA DF6D08 fild qword ptr [ebp+$08]
004612AD DEF9 fdivp st(1)
004612AF DF7DF8 fistp qword ptr [ebp-$08]
004612B2 9B wait
004612B3 8B45F8 mov eax,[ebp-$08]
004612B6 8B55FC mov edx,[ebp-$04]
004612B9 59 pop ecx
004612BA 59 pop ecx
004612BB 5D pop ebp
004612BC C21800 ret $0018
Следующий оператор выдает 256204778801521550, что представляется правильным.
writeln(MulDiv64($aaaaaaaaaaaaaaa, $555555555555555, $1000000000000000));
Если вы хотите реализовать это как встроенную сборку VC++, возможно, вам понадобится немного изменить флаги округления по умолчанию, чтобы выполнить то же самое, я не знаю - у меня не было необходимости выяснять - еще:)
Для кода 64-битного режима вы можете реализовать умножение 64*64=128 аналогично реализации деления 128/64=64:64 здесь.
Для 32-битного кода это будет более сложным, потому что нет инструкции ЦП, которая бы делала умножение или деление таких длинных операндов в 32-битном режиме, и вам придется объединить несколько меньших умножений в большее и переопределить длинное деление.
Вы можете использовать код из этого ответа в качестве основы для построения длинного деления.
Конечно, если ваши делители всегда меньше чем 232 (или еще лучше 216), вы можете сделать намного более быстрое деление в 32-битном режиме путем объединения нескольких делений (делить старшие 64 (или 32) бита делимого на 32-разрядный (или 16-разрядный) делитель, затем объедините остаток от этого деления со следующими 64 (32) битами дивиденда, разделите его на делитель и продолжайте делать это, пока вы не израсходуете весь дивиденд). Кроме того, если делитель большой, но может быть разложен на достаточно малые числа, деление на его факторы с использованием этого цепочечного деления будет лучше, чем классическое решение с циклическим циклом.
Учтите, что вы хотите умножить a от b а затем разделить на d:
uint64_t LossyMulDiv64(uint64_t a, uint64_t b, uint64_t d)
{
long double f = long double(b)/d;
uint64_t highPart = uint64_t((a & ~0xffffffff) * f + 0.5);
uint64_t lowPart = uint64_t((a & 0xffffffff) * f + 0.5);
return highPart + lowPart;
}
Этот код разбивает значение a на более высокие и более низкие 32-битные части, затем умножает 32-битные части отдельно на 52-битное точное соотношение b в dокругляет частичные умножения и суммирует их обратно до целого числа. Некоторая точность все еще теряется, но результат более точен, чем просто return a * double(b) / d;,
Вот метод аппроксимации, который вы можете использовать: (полная точность, если a > 0x7FFFFFFF или b > 0x7FFFFFFF и c больше, чем a или b)
constexpr int64_t muldiv(int64_t a, int64_t b, int64_t c, unsigned n = 0) {
return (a < 0x7FFFFFFF && b < 0x7FFFFFFF) ? (a * b) / c : (n != 2) ? (c <= a) ? ((a / c) * b + muldiv(b, a % c, c, n + 1)) : muldiv(a, b, c / 2) / 2 : 0;
}
Модуль используется для определения потери точности, которая затем снова включается в функцию. Это похоже на классический алгоритм деления.
2 был выбран потому что (x % x) % x = x % x,