Я работаю над собственным классом C ++ / CLI, который выполняет целочисленную арифметику со значениями с множественной точностью.Отдельные целые числа представлены массивами 64-битных целых чисел без знака.Знак представлен логическим значением, отрицательные значения хранятся вместе с абсолютными значениями, а не как два дополнения.Это значительно облегчает решение проблем со знаком.В настоящее время я оптимизирую операцию умножения.Я уже провел несколько раундов оптимизации, но все же моя функция требует в два раза больше времени * оператора * двух значений .NET BigInteger , что показывает, что еще есть значительный потенциал для дальнейшей оптимизации.
Прежде чем обратиться за помощью, позвольте мне показать вам, что я уже пробовал.Моей первой попыткой был наивный подход: умножить пары всех 64-битных элементов, используя элементарное умножение с 64 на 128 бит, и сдвинуть / сложить результаты.Я не показываю код здесь, потому что он был очень медленным.Следующей попыткой был рекурсивный алгоритм "разделяй и властвуй" , который оказался намного лучше.В моей реализации оба операнда рекурсивно разделяются посередине, пока не останутся два 64-битных значения.Они умножаются, давая 128-битный результат.Собранные элементарные результаты сдвигаются / складываются по всем уровням рекурсии для получения окончательного результата.Этот алгоритм, вероятно, выигрывает от того факта, что нужно вычислять гораздо меньше 64-128-битных элементарных продуктов, что является основным узким местом.
Итак, вот мой код.Первый фрагмент показывает точку входа верхнего уровня:
// ----------------------------------------------------------------------------
// Multi-precision multiplication, using a recursive divide-and-conquer plan:
// Left split: (a*2^k + b)i = ai*2^k + bi
// Right split: a(i*2^k + j) = ai*2^k + aj
public: static UINT64* Mul (UINT64* pu8Factor1,
UINT64* pu8Factor2,
UINT64 u8Length1,
UINT64 u8Length2,
UINT64& u8Product)
{
UINT64* pu8Product;
if ((u8Length1 > 0) && (u8Length2 > 0))
{
pu8Product = _SnlMemory::Unsigned ((u8Length1 * u8Length2) << 1);
u8Product = Mul (pu8Product, pu8Factor1, 0, u8Length1,
pu8Factor2, 0, u8Length2);
}
else
{
pu8Product = _SnlMemory::Unsigned (0);
u8Product = 0;
}
return pu8Product;
}
Коэффициенты передаются в виде UINT64* указателей массива, а длины указываются отдельно как количество элементов UINT64 в соответствующем массиве.Функция выделяет блок памяти, достаточно большой, чтобы содержать значение максимальной ожидаемой длины, которое также используется в качестве блокнота для временных подчиненных результатов.Функция вызывает другую функцию Mul, которая выполняет рекурсивную оценку и возвращает количество элементов UINT64, фактически использованных в конечном результате.
Это рекурсивная часть «разделяй и разделяй»алгоритм:
// ----------------------------------------------------------------------------
// Recursively expand the arbitrary-precision multiplication to the sum of a
// series of elementary 64-to-128-bit multiplications.
private: static UINT64 Mul (UINT64* pu8Product,
UINT64* pu8Factor1,
UINT64 u8Offset1,
UINT64 u8Length1,
UINT64* pu8Factor2,
UINT64 u8Offset2,
UINT64 u8Length2)
{
UINT64 *pu8Lower, u8Lower, *pu8Upper, u8Upper, u8Split;
UINT64 u8Product = 0;
if (u8Length1 > 1)
{
// left split: (a*2^k + b)i = ai*2^k + bi
u8Split = u8Length1 >> 1;
u8Lower = Mul (pu8Lower = pu8Product,
pu8Factor1, u8Offset1, u8Split, // bi
pu8Factor2, u8Offset2, u8Length2);
u8Upper = Mul (pu8Upper = pu8Product + ((u8Split * u8Length2) << 1),
pu8Factor1, u8Offset1 + u8Split, // ai
u8Length1 - u8Split,
pu8Factor2, u8Offset2, u8Length2);
u8Product = Mul (u8Split, pu8Lower, u8Lower, pu8Upper, u8Upper);
}
else if (u8Length2 > 1)
{
// right split: a(i*2^k + j) = ai*2^k + aj
u8Split = u8Length2 >> 1;
u8Lower = Mul (pu8Lower = pu8Product,
pu8Factor1, u8Offset1, u8Length1, // aj
pu8Factor2, u8Offset2, u8Split);
u8Upper = Mul (pu8Upper = pu8Product + ((u8Length1 * u8Split) << 1),
pu8Factor1, u8Offset1, u8Length1, // ai
pu8Factor2, u8Offset2 + u8Split,
u8Length2 - u8Split);
u8Product = Mul (u8Split, pu8Lower, u8Lower, pu8Upper, u8Upper);
}
else // recursion base: 64-to-128-bit multiplication
{
AsmMul1 (pu8Factor1 [u8Offset1],
pu8Factor2 [u8Offset2],
u8Lower, u8Upper);
if (u8Upper > 0)
{
pu8Product [u8Product++] = u8Lower;
pu8Product [u8Product++] = u8Upper;
}
else if (u8Lower > 0)
{
pu8Product [u8Product++] = u8Lower;
}
}
return u8Product;
}
В первой условной ветви левый операнд разбивается.Во втором правый операнд разделен.Третья ветвь - это база рекурсии, которая вызывает элементарную процедуру умножения:
; -----------------------------------------------------------------------------
; 64-bit to 128-bit multiplication, using the x64 MUL instruction
AsmMul1 proc ; ?AsmMul1@@$$FYAX_K0AEA_K1@Z
; ecx : Factor1
; edx : Factor2
; [r8] : ProductL
; [r9] : ProductH
mov rax, rcx ; rax = Factor1
mul rdx ; rdx:rax = Factor1 * Factor2
mov qword ptr [r8], rax ; [r8] = ProductL
mov qword ptr [r9], rdx ; [r9] = ProductH
ret
AsmMul1 endp
Это простой ASM PROC, который использует инструкцию CPU MUL для умножения с 64 на 128 бит.Я пробовал несколько других кандидатов в ASM и C ++, и этот был наиболее эффективным.
Последняя часть - это «завоевать» часть алгоритма «разделяй и властвуй»:
// ----------------------------------------------------------------------------
// Shift-add recombination of the results of two partial multiplications.
private: static UINT64 Mul (UINT64 u8Split,
UINT64* pu8Lower,
UINT64 u8Lower,
UINT64* pu8Upper,
UINT64 u8Upper)
{
FLAG fCarry;
UINT64 u8Count, u8Lower1, u8Upper1;
UINT64 u8Product = u8Lower;
if (u8Upper > 0)
{
u8Count = u8Split + u8Upper;
fCarry = false;
for (u8Product = u8Split; u8Product < u8Count; u8Product++)
{
u8Lower1 = u8Product < u8Lower ? pu8Lower [u8Product] : 0;
u8Upper1 = pu8Upper [u8Product - u8Split];
if (fCarry)
{
pu8Lower [u8Product] = u8Lower1 + u8Upper1 + 1;
fCarry = u8Lower1 >= MAX_UINT64 - u8Upper1;
}
else
{
pu8Lower [u8Product] = u8Lower1 + u8Upper1;
fCarry = u8Lower1 > MAX_UINT64 - u8Upper1;
}
}
if (fCarry)
{
pu8Lower [u8Product++] = 1;
}
}
return u8Product;
}
Здесь добавляются два частичных результата, причем второй операнд смещен влево на «коэффициент разделения» соответствующего шага рекурсии.
Я потратил значительное время на оптимизацию кода для скорости, причем довольно многоуспех, но теперь я достиг точки, где я не вижу дальнейших возможностей, кроме использования совершенно другого алгоритма.Однако, поскольку я не специалист по числовым трюкам, я застрял здесь.
Надеясь на некоторые интересные идеи о том, как улучшить этот расчет ...
РЕДАКТИРОВАТЬ 2019-03-26: Ну, иногда кажется, что лучше не пытаться быть умным ... после нескольких дополнительных попыток оптимизации, некоторые из которых даже умеренно успешны, я попытался написать настоящую глупую версию умножения, котораяпросто использует вычислительные возможности встроенных компиляторов _umul128 и _addcarry_u64.Код предельно прост:
public: static UINT64* Mul (UINT64* pu8Factor1,
UINT64* pu8Factor2,
UINT64 u8Length1,
UINT64 u8Length2,
UINT64& u8Product)
{
u8Product = u8Length1 + u8Length2;
CHAR c1Carry1, c1Carry2;
UINT64 u8Offset, u8Offset1, u8Offset2, u8Item1, u8Item2, u8Lower, u8Upper;
UINT64* pu8Product = _SnlMemory::Unsigned (u8Product);
if (u8Product > 0)
{
for (u8Offset1 = 0; u8Offset1 < u8Length1; u8Offset1++)
{
u8Offset = u8Offset1;
u8Item1 = pu8Factor1 [u8Offset1];
u8Item2 = 0;
c1Carry1 = 0;
c1Carry2 = 0;
for (u8Offset2 = 0; u8Offset2 < u8Length2; u8Offset2++)
{
u8Lower = _umul128 (u8Item1, pu8Factor2 [u8Offset2], &u8Upper);
c1Carry1 = _addcarry_u64 (c1Carry1, u8Item2, u8Lower,
&u8Item2);
c1Carry2 = _addcarry_u64 (c1Carry2, u8Item2,
pu8Product [u8Offset],
pu8Product + u8Offset);
u8Item2 = u8Upper;
u8Offset++;
}
if (c1Carry1 != 0)
{
c1Carry2 = _addcarry_u64 (c1Carry2, u8Item2 + 1,
pu8Product [u8Offset],
pu8Product + u8Offset);
}
else if (u8Item2 != 0)
{
c1Carry2 = _addcarry_u64 (c1Carry2, u8Item2,
pu8Product [u8Offset],
pu8Product + u8Offset);
}
}
if (pu8Product [u8Product - 1] == 0)
{
u8Product--;
}
}
return pu8Product;
}
Создает буфер результатов в куче, достаточно большой, чтобы вместить максимальный размер продукта, и выполняет элементарные 64-zo-128-битные _umul128 умножения в двух вложенных циклах в сочетании с двумя потоками переноса ряби.дополнения с использованием _addcarry_u64.Производительность этой версии - лучшая из всех, что я пробовал до сих пор!Это примерно в 10 раз быстрее, чем эквивалентное умножение .NET BigInteger, поэтому в итоге я добился ускорения в 20 раз.