Мне нужен быстрый 96-битный на 64-битный алгоритм деления для математической библиотеки с фиксированной точкой

Question

В настоящее время я пишу быструю математическую библиотеку с фиксированной точкой 32.32. Мне удалось правильно добавить сложение, вычитание и умножение, но я довольно застрял в делении.

Небольшое напоминание для тех, кто не может вспомнить: число с фиксированной точкой 32.32 - это число, содержащее 32 бита целой части и 32 бита дробной части.

Лучший алгоритм, который я придумал, требует 96-битного целочисленного деления, для чего компиляторы обычно не имеют встроенных модулей.

Так или иначе, вот оно:

G = 2^32

notation: x is the 64-bit fixed-point number, x1 is its low nibble and x2 is its high

G*(a/b) = ((a1 + a2*G) / (b1 + b2*G))*G      // Decompose this

G*(a/b) = (a1*G) / (b1*G + b2) + (a2*G*G) / (b1*G + b2)

Как видите, (a2*G*G) гарантированно будет больше, чем обычное 64-битное целое число. Если бы uint128_t действительно поддерживался моим компилятором, я бы просто сделал следующее:

((uint128_t)x << 32) / y)

Ну, это не так, и мне нужно решение. Спасибо за вашу помощь.

Answer 1

Вы можете разбить большее деление на несколько кусков, которые делят с меньшим количеством битов. Как уже упоминалось в другом постере, алгоритм можно найти в TAOCP от Кнута.

Однако, не нужно покупать книгу!

На веб-сайте восхищения хакерами есть код, который реализует алгоритм на языке C. Он написан для выполнения 64-битных беззнаковых делений, используя только 32-битную арифметику, поэтому вы не можете напрямую вырезать и вставить код. Чтобы получить от 64 до 128 бит, вы должны расширить все типы, маски и константы на два, например. короткое становится int, 0xffff становится 0xffffffffll ect.

После этого простого и легкого изменения вы сможете делать 128-битные деления.

Код здесь: http://www.hackersdelight.org/HDcode/divlu.c (может плохо переноситься в веб-браузере из-за окончания строки. Если это так, просто сохраните код и откройте его в блокноте или около того).

Поскольку вашим самым большим значениям требуется только 96 бит, одно из 64-битных делений будет всегда возвращать ноль, так что вы даже можете немного упростить код.

Да, и прежде, чем я это забуду: код работает только со значениями без знака. Чтобы преобразовать деление со знаком в беззнаковое, вы можете сделать что-то вроде этого (стиль псевдокода):

fixpoint Divide (fixpoint a, fixpoint b)
{
  // check if the integers are of different sign:
  fixpoint sign_difference = a ^ b; 

  // do unsigned division:
  fixpoint x = unsigned_divide (abs(a), abs(b));

  // if the signs have been different: negate the result.
  if (sign_difference < 0)
  {
     x = -x;
  }

  return x;
}

Стоит также проверить сам сайт: http://www.hackersdelight.org/

Надеюсь, это поможет.

Кстати - отличное задание, над которым вы работаете ... Не могли бы вы рассказать нам, для чего вам нужна библиотека с фиксированной точкой?

---

Кстати - обычный алгоритм сдвига и вычитания для деления тоже подойдет.

Если вы нацелены на x86, вы можете реализовать его с помощью встроенных функций MMX или SSE. Алгоритм опирается только на примитивные операции, поэтому он также может выполнять довольно быстро.

Answer 2

Быстро-грязно.

Делайте A / B деление с плавающей запятой двойной точности. Это дает вам C ~ = A / B. Это только приблизительно из-за точности с плавающей точкой и 53 битов мантиссы.

Округлите C до представимого числа в вашей системе с фиксированной точкой.

Теперь вычислите (снова с вашей фиксированной точкой) D = A-C * B. Это должно иметь значительно меньшую величину, чем А.

Повторите, теперь вычисление D / B с плавающей запятой. Опять округлите ответ до целого числа. Сложите каждый результат деления вместе, как вы идете. Вы можете остановиться, когда ваш остаток настолько мал, что ваше деление с плавающей запятой возвращает 0 после округления.

Ты еще не закончил. Теперь вы очень близки к ответу, но разделения не были точными. Чтобы завершить, вам нужно выполнить бинарный поиск. Используя (очень хорошую) начальную оценку, посмотрите, не увеличивает ли она ошибку. Вы в основном хотите заключить в скобки правильный ответ и продолжать делить диапазон пополам с помощью новых тестов.

Да, вы можете выполнить итерацию Ньютона здесь, но бинарный поиск, вероятно, будет проще, поскольку вам нужны только простые умножения и добавления с использованием существующего инструментария точности 32.32.

Это не самый эффективный метод, но на данный момент его проще всего кодировать.

Answer 3

Лучший самонастраивающийся ответ :
Простите за ответ C #, но во всех случаях должно работать следующее. Вероятно, существует решение, которое находит правильные сдвиги для более быстрого использования, но сейчас мне придется думать гораздо глубже, чем я могу. Это должно быть разумно эффективно, хотя:

int upshift = 32;
ulong mask = 0xFFFFFFFF00000000;
ulong mod = x % y;
while ((mod & mask) != 0)
{
     // Current upshift of the remainder would overflow... so adjust
     y >>= 1;
     mask <<= 1;
     upshift--;

     mod = x % y;
}
ulong div = ((x / y) << upshift) + (mod << upshift) / y;

Простой, но небезопасный ответ :
Это вычисление может вызвать переполнение при повышении смещения остатка x % y, если для этого остатка установлены какие-либо биты в старших 32 битах, что приведет к неправильному ответу.

((x / y) << 32) + ((x % y) << 32) / y

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

Вторая часть вычисляет младшие биты из остатка от старшего разряда (бит, который не может быть разделен далее), сдвинут вверх и затем разделен.

Answer 4

Мне нравится ответ Нильса, который, вероятно, самый лучший. Это просто длинное деление, как мы все учили в начальной школе, за исключением того, что цифры - это база 2 ^ 32 вместо базы 10.

Однако вы можете также рассмотреть возможность использования метода аппроксимации Ньютона для деления:

  x := x (N + N - N * D * x)

где N - числитель, а D - демонинатор.

Это просто использует умножения и сложения, которые у вас уже есть, и очень быстро приближается к 1 ULP точности. С другой стороны, вы не сможете получить точный ответ 0,5-ULP во всех случаях.

В любом случае хитрый бит обнаруживает и обрабатывает переполнения.