Какой самый быстрый алгоритм разделения на языке ассемблера x86-64 для огромных чисел?

Я пишу библиотеку кода на языке ассемблера x86-64, чтобы обеспечить все обычные побитовые, сдвиговые, логические, сравнивающие, арифметические и математические функции для s0128, s0256, s0512, s1024 целочисленных типов со знаком и f0128, f0256, f0512, f1024 с плавающей точкой. -точечные типы. Пока я работаю над целочисленными типами со знаком, потому что функции с плавающей запятой, скорее всего, будут вызывать некоторые внутренние процедуры, написанные для целочисленных типов.

До сих пор я написал и протестировал функции для выполнения различных унарных операторов, операторов сравнения, а также операторов сложения, вычитания и умножения.

Теперь пытаюсь решить, как реализовать функции для операторов деления.

Моей первой мыслью было: «Ньютон-Рафсон, должно быть, лучший подход». Почему? Потому что он сходится очень быстро при хорошем начальном значении (начальное предположение), и я полагаю, что смогу выяснить, как выполнить встроенную 64-битную инструкцию деления для операндов, чтобы получить отличное начальное значение. Фактически, если начальное значение является точным до 64 бит, для получения правильных ответов нужно всего лишь взять:

`s0128` : 1~2 iterations : (or 1 iteration  plus 1~2 "test subtracts")
`s0256` : 2~3 iterations : (or 2 iterations plus 1~2 "test subtracts")
`s0512` : 3~4 iterations : (or 3 iterations plus 1~2 "test subtracts")
`s1024` : 4~5 iterations : (or 4 iterations plus 1~2 "test subtracts")

Однако немного больше обдумывая этот вопрос, заставляет меня задуматься. Например, я вспоминаю написанную мной основную подпрограмму, которая выполняет операцию умножения для всех больших целых типов:

s0128 :   4 iterations ==   4 (128-bit = 64-bit * 64-bit) multiplies +  12 adds
s0256 :  16 iterations ==  16 (128-bit = 64-bit * 64-bit) multiplies +  48 adds
s0512 :  64 iterations ==  64 (128-bit = 64-bit * 64-bit) multiplies + 192 adds
s1024 : 256 iterations == 256 (128-bit = 64-bit * 64-bit) multiplies + 768 adds

Рост числа операций для более широких типов данных весьма существенен, даже несмотря на то, что цикл довольно короткий и эффективный (в том числе с точки зрения кеширования). Этот цикл записывает каждую 64-битную часть результата только один раз и никогда не считывает какую-либо часть результата для дальнейшей обработки.

Это заставило меня задуматься о том, могут ли более традиционные алгоритмы разделения типа сдвиг и вычитание быть быстрее, особенно для больших типов.

Моя первая мысль была такой:

result = dividend / divisor                  // if I remember my terminology
remainder = dividend - (result * divisor)    // or something along these lines

№1: Чтобы вычислить каждый бит, обычно делитель вычитается из делимого ЕСЛИ делитель меньше или равен делимому. Что ж, обычно мы можем определить, что делитель определенно меньше или определенно больше дивиденда, проверяя только их наиболее значимые 64-битные части. Только тогда, когда эти ms64-битные части равны, процедура должна проверять следующие младшие 64-битные части, и только когда они равны, мы должны проверять даже меньшие, и так далее. Следовательно, почти на каждой итерации (вычислении каждого бита результата) мы можем значительно сократить количество инструкций, выполняемых для вычисления этого теста.

№2: Однако ... в среднем, примерно в 50% случаев мы обнаруживаем, что нам нужно вычесть делитель из делимого, поэтому нам все равно придется вычитать их ширину полностью. В этом случае мы фактически выполнили больше инструкций, чем при обычном подходе (где мы сначала вычитаем их, а затем проверяем флаги, чтобы определить, равен ли делитель ‹делимому). Таким образом, половину времени мы реализуем экономию, а половину - убыток. На больших типах, таких как s1024 (который содержит -16- 64-битные компоненты), экономия значительна, а потери небольшие, поэтому такой подход должен обеспечить большую чистую экономию. На крошечных типах, таких как s0128 (который содержит -2- 64-битные компоненты), экономия мала, а потери значительны, но не огромны.

Итак, мой вопрос: «Каковы наиболее эффективные алгоритмы разделения», учитывая:

#1: modern x86-64 CPUs like FX-8350
#2: executing in 64-bit mode only (no 32-bit)
#3: implementation entirely in assembly-language
#4: 128-bit to 1024-bit integer operands (nominally signed, but...)

ПРИМЕЧАНИЕ. Я предполагаю, что реальная реализация будет работать только с целыми числами без знака. В случае умножения оказалось проще и эффективнее (возможно) преобразовать отрицательные операнды в положительные, затем выполнить беззнаковое умножение, а затем отменить результат, если ровно один исходный операнд был отрицательным. Однако я рассмотрю целочисленный алгоритм со знаком, если он эффективен.

ПРИМЕЧАНИЕ. Если лучшие ответы для моих типов с плавающей запятой (f0128, f0256, f0512, f1024) отличаются, объясните, почему.

ПРИМЕЧАНИЕ. Моя внутренняя подпрограмма беззнакового умножения, которая выполняет операцию умножения для всех этих целочисленных типов данных, дает результат с двойной шириной. Другими словами:

u0256 = u0128 * u0128     // cannot overflow
u0512 = u0256 * u0256     // cannot overflow
u1024 = u0512 * u0512     // cannot overflow
u2048 = u1024 * u1024     // cannot overflow

Моя библиотека кода предлагает две версии умножения для каждого целочисленного типа данных со знаком:

s0128 = s0128 * s0128     // can overflow (result not fit in s0128)
s0256 = s0256 * s0256     // can overflow (result not fit in s0256)
s0512 = s0512 * s0512     // can overflow (result not fit in s0512)
s1024 = s1024 * s1024     // can overflow (result not fit in s1024)

s0256 = s0128 * s0128     // cannot overflow
s0512 = s0256 * s0256     // cannot overflow
s1024 = s0512 * s0512     // cannot overflow
s2048 = s1024 * s1024     // cannot overflow

Это согласуется с политикой моей библиотеки кода «никогда не терять точности» и «никогда не переполняться» (ошибки возвращаются, когда ответ недействителен из-за потери точности или из-за переполнения / потери значимости). Однако при вызове функций, возвращающих значение двойной ширины, такие ошибки возникать не могут.