Самый точный способ выполнить комбинированную операцию умножения и деления в 64-битном режиме?
Каков наиболее точный способ выполнения операции умножения и деления для 64-разрядных целых чисел, которые работают как в 32-разрядных, так и в 64-разрядных программах (в Visual С++)? (В случае переполнения мне нужен результат mod 2 64.)
(Я ищу что-то вроде MulDiv64, за исключением того, что этот использует встроенную сборку, которая работает только в 32-битных программах.)
Очевидно, что при нажатии на double и обратно возможно, но мне интересно, есть ли более точный способ, который не слишком сложный. (т.е. я не ищу здесь арифметические библиотеки произвольной точности!)
Ответы
Ответ 1
Так как это помечено Visual С++, я дам решение, которое нарушает специфические для MSVC встроенные функции.
Этот пример довольно сложный. Это очень упрощенная версия того же алгоритма, который используется GMP и java.math.BigInteger для большого деления.
Хотя у меня есть более простой алгоритм, он, вероятно, примерно на 30 раз медленнее.
Это решение имеет следующие ограничения/поведение:
- Для этого требуется x64. Он не будет компилироваться на 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
Ответ 2
Вам просто нужны 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
}
Ответ 3
Это ответ вики сообщества, поскольку это действительно просто куча указателей на другие документы/ссылки (я не могу опубликовать соответствующий код).
Умножение двух 64-битных ints на 128-битный результат довольно легко, используя прямое применение карандаша и бумажной техники, которые каждый изучает в начальной школе.
Комментарий GregS верен: в разделе "Искусство компьютерного программирования, второе издание, том 2/" Семинумерные алгоритмы "в конце раздела 4.3.1" Множество прецизионных арифметических/классических алгоритмов "(стр. 255 - 265) копия). Это нелегко прочитать, по крайней мере, не для кого-то вроде меня, который забыл большинство математик за пределами алгебры 7-го класса. Как раз перед, Кнут также охватывает сторону умножения вещей.
Некоторые другие варианты идей (эти примечания предназначены для алгоритмов деления, но большинство также обсуждает умножение):
- Джек Креншоу раскрывает алгоритмы деления Кнута более читаемым образом в серии статей из журнала Embedded System Programming 1997 (к сожалению, в моих заметках нет точных проблем). К сожалению, статьи из старых вопросов ESP нелегко найти в Интернете. Если у вас есть доступ к университетской библиотеке, возможно, вам понадобятся некоторые проблемы с обратной связью или копия библиотеки CD-ROM ESP.
- Томас Родеффер из исследования Microsoft опубликовал статью о подразделении Software Integer: http://research.microsoft.com/pubs/70645/tr-2008-141.pdf
- Статья Карла Хассельстрема "Быстрое разделение больших целых чисел": http://www.treskal.com/kalle/exjobb/original-report.pdf
- Randall Hyde "Искусство языка ассемблера" (http://webster.cs.ucr.edu/AoA/Windows/HTML/AoATOC.html), в частности раздел четвертый раздел 4.2.5 (расширенный прецизионный отдел): http://webster.cs.ucr.edu/AoA/Windows/HTML/AdvancedArithmetica2.html#998729, это вариант Hyde для ассемблера x86, но также есть псевдокод и достаточно объяснений для переноса алгоритма на C. Это тоже медленное - выполнение бит-по-бит...
Ответ 4
Для этого вам не нужна арифметика произвольной точности. Вам нужна только 128-разрядная арифметика. То есть вам нужно 64 * 64 = 128 умножения и 128/64 = 64 деления (с надлежащим поведением переполнения). Это не так сложно реализовать вручную.
Ответ 5
Хорошо, вы можете нарезать 64-разрядные операнды на 32-битные куски (низкая и высокая часть). Затем сделайте операцию, которую вы хотите. Все промежуточные результаты будут меньше 64 бит и поэтому могут храниться в типах данных, которые у вас есть.
Ответ 6
У вас есть тип COMP (64-разрядный целочисленный тип на основе x87) в вашем распоряжении в VС++? Я использовал его иногда в 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 (или что-то, что он называет теперь) не имеет четкого синтаксического различия между типом casting (reinterpret) и преобразованием типа значения.)
В любом случае, 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));
Если вы хотите реализовать это как встроенную сборку VС++, возможно, вам понадобится выполнить некоторую настройку флажков округления по умолчанию, чтобы выполнить одно и то же, я не знаю - у меня не было необходимости узнайте - пока:)
Ответ 7
Для 64-битного режима кода вы можете реализовать умножение 64 * 64 = 128 аналогично реализации 128/64 = 64: 64 раздела здесь.
Для 32-битного кода он будет более сложным, потому что нет инструкции по процессору, которая будет делать умножение или деление таких длинных операндов в 32-битном режиме, и вам придется объединить несколько меньших умножений в более крупный и переопределить длинное разделение.
Вы можете использовать код этого ответа в качестве основы для построения длинного разделения.
Конечно, если ваши разделители всегда меньше 2 32 (или еще лучше 2 16), вы можете сделать гораздо более быстрое деление в 32-битном режиме путем цепочки несколько делений (разделите самые значительные 64 (или 32) бита дивиденда на 32-разрядный (или 16-разрядный) делитель, затем объедините остаток от этого деления со следующими 64 (32) битами дивиденда и разделите это по делителю и продолжайте делать это, пока не будете использовать весь дивиденд). Кроме того, если делитель большой, но может быть учтен в достаточно малые числа, деление на его коэффициенты с использованием этого цепного деления будет лучше, чем классическое решение петли.
Ответ 8
Здесь метод аппроксимации, который вы можете использовать: (полная точность, если 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.
Ответ 9
В C нет 32x32→64 или 64x64→128 умножений/делений. Результат всегда имеет тот же размер, что и самый большой из множителей и множитель. Это означает, что только int32 x int32 → int32 и int64 x int64 → int64. В случае переполнения результатом являются младшие биты. То же самое для делений. Таким образом, на самом деле вам нужно объявить переменные как __int64 или передать один из множителей/делителей в __int64, и результат будет ниже 64 бит по мере необходимости.
__int64 a, b, c;
c = a * b;
__int32 d, e;
b = __int64(e) / d;
Ответ 10
Считайте, что вы хотите умножить 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;.
