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

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

Однако нет необходимости покупать книгу!

На веб-сайте hackers delight есть код, который реализует алгоритм на C. Он написан для выполнения 64-битное деление без знака с использованием только 32-битной арифметики, поэтому вы не можете напрямую вырезать и вставить код. Чтобы перейти от 64 к 128 битам, вы должны расширить все типы, маски и константы на два, например, short становится 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/

Hope это помогает.

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

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

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

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

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

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

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

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

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

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

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

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

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

Quick -n-dirty.

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

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

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

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

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

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

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