Эффективное вычисление 2 ** 64/делителя через быстродействующую обратную точку с плавающей запятой
В настоящее время я изучаю способы использования быстрой одноточечной функции с плавающей запятой различных современных процессоров для вычисления начального приближения для 64-разрядного беззнакового целочисленного деления на основе итераций Newton-Raphson с фиксированной точкой. Это требует вычисления 2 64/divisor, насколько это возможно, где начальное приближение должно быть меньше или равно математическому результату на основе требований следующих итераций с фиксированной точкой. Это означает, что это вычисление должно обеспечить недооценку. В настоящее время у меня есть следующий код, который хорошо работает, основываясь на обширном тестировании:
#include <stdint.h> // import uint64_t
#include <math.h> // import nextafterf()
uint64_t divisor, recip;
float r, s, t;
t = uint64_to_float_ru (divisor); // ensure t >= divisor
r = 1.0f / t;
s = 0x1.0p64f * nextafterf (r, 0.0f);
recip = (uint64_t)s; // underestimate of 2**64 / divisor
Хотя этот код является функциональным, на большинстве платформ он не совсем быстр. Одно очевидное улучшение, требующее немного машинного кода, заключается в замене деления r = 1.0f / t на код, который использует быструю с плавающей запятой, предоставляемую аппаратным обеспечением. Это может быть дополнено итерацией для получения результата, который находится в пределах 1 ulp математического результата, поэтому недооценка создается в контексте существующего кода. Пример реализации для x86_64:
#include <xmmintrin.h>
/* Compute 1.0f/a almost correctly rounded. Halley iteration with cubic convergence */
inline float fast_recip_f32 (float a)
{
__m128 t;
float e, r;
t = _mm_set_ss (a);
t = _mm_rcp_ss (t);
_mm_store_ss (&r, t);
e = fmaf (r, -a, 1.0f);
e = fmaf (e, e, e);
r = fmaf (e, r, r);
return r;
}
Реализации nextafterf() обычно не оптимизированы по производительности. На платформах, где есть средства для быстрого повторного преобразования IEEE 754 binary32 в int32 и наоборот, через intrinsics float_as_int() и int_as_float(), мы можем комбинировать использование nextafterf() и масштабирования следующим образом:
s = int_as_float (float_as_int (r) + 0x1fffffff);
Предполагая, что эти подходы возможны на данной платформе, это оставляет нас с преобразованиями между float и uint64_t в качестве основных препятствий. Большинство платформ не предоставляют команду, которая выполняет преобразование с uint64_t в float со статическим режимом округления (здесь: в направлении положительной бесконечности = вверх), а некоторые не предлагают никаких команд для преобразования между uint64_t и плавающей точечные типы, что делает это узким местом производительности.
t = uint64_to_float_ru (divisor);
r = fast_recip_f32 (t);
s = int_as_float (float_as_int (r) + 0x1fffffff);
recip = (uint64_t)s; /* underestimate of 2**64 / divisor */
Портативная, но медленная реализация uint64_to_float_ru использует динамические изменения в режиме округления FPU:
#include <fenv.h>
#pragma STDC FENV_ACCESS ON
float uint64_to_float_ru (uint64_t a)
{
float res;
int curr_mode = fegetround ();
fesetround (FE_UPWARD);
res = (float)a;
fesetround (curr_mode);
return res;
}
Я рассмотрел различные методы разделения и бит-скругления для обработки конверсий (например, округление на целочисленной стороне, затем используйте обычное преобразование в float, которое использует круглый режим IEEE 754 с округлением до ближайшего -или-четное), но накладные расходы, которые это создает, делают это вычисление быстрым обратным возвратом с плавающей запятой с точки зрения производительности. Как бы то ни было, похоже, что мне лучше было бы начать начальное приближение, используя классический LUT с интерполяцией или полиномиальную аппроксимацию с фиксированной точкой, и последуйте за ними с помощью 32-битного шага Newton-Raphson с фиксированной точкой.
Есть ли способы повысить эффективность моего текущего подхода? Интересны портативные и полупортативные способы, связанные с встроенными функциями для конкретных платформ (в частности, для x86 и ARM в качестве доминирующих в настоящее время архитектур ЦП). Компиляция для x86_64 с использованием компилятора Intel при очень высокой оптимизации (/O3 /QxCORE-AVX2 /Qprec-div-) вычисление начального приближения занимает больше инструкций, чем итерация, которая занимает около 20 инструкций. Ниже приведен полный код деления для ссылки, показывающий приближение в контексте.
uint64_t udiv64 (uint64_t dividend, uint64_t divisor)
{
uint64_t temp, quot, rem, recip, neg_divisor = 0ULL - divisor;
float r, s, t;
/* compute initial approximation for reciprocal; must be underestimate! */
t = uint64_to_float_ru (divisor);
r = 1.0f / t;
s = 0x1.0p64f * nextafterf (r, 0.0f);
recip = (uint64_t)s; /* underestimate of 2**64 / divisor */
/* perform Halley iteration with cubic convergence to refine reciprocal */
temp = neg_divisor * recip;
temp = umul64hi (temp, temp) + temp;
recip = umul64hi (recip, temp) + recip;
/* compute preliminary quotient and remainder */
quot = umul64hi (dividend, recip);
rem = dividend - divisor * quot;
/* adjust quotient if too small; quotient off by 2 at most */
if (rem >= divisor) quot += ((rem - divisor) >= divisor) ? 2 : 1;
/* handle division by zero */
if (divisor == 0ULL) quot = ~0ULL;
return quot;
}
umul64hi(), как правило, относится к встроенной платформе или к биту встроенного ассемблерного кода. В x86_64 я в настоящее время использую эту реализацию:
inline uint64_t umul64hi (uint64_t a, uint64_t b)
{
uint64_t res;
__asm__ (
"movq %1, %%rax;\n\t" // rax = a
"mulq %2;\n\t" // rdx:rax = a * b
"movq %%rdx, %0;\n\t" // res = (a * b)<63:32>
: "=rm" (res)
: "rm"(a), "rm"(b)
: "%rax", "%rdx");
return res;
}
Ответы
Ответ 1
Это решение объединяет две идеи:
- Вы можете преобразовать в плавающую точку, просто переинтерпретируя биты как с плавающей запятой и вычитая константу, если число находится в определенном диапазоне. Поэтому добавьте константу, переинтерпретируйте, а затем вычтите эту константу. Это даст усеченный результат (который поэтому всегда меньше или равен желаемому значению).
- Вы можете аппроксимировать взаимность, отрицая как экспонента, так и мантиссу. Это может быть достигнуто путем интерпретации битов как int.
Вариант 1 здесь работает только в определенном диапазоне, поэтому мы проверяем диапазон и корректируем используемые константы. Это работает в 64 битах, потому что желаемый поплавок имеет только 23 бита точности.
Результат этого кода будет двойным, но преобразование в float тривиально и может быть выполнено на битах или напрямую, в зависимости от аппаратного обеспечения.
После этого вы захотите выполнить итерацию Newton-Raphson.
Большая часть этого кода просто преобразуется в магические числа.
double
u64tod_inv( uint64_t u64 ) {
__asm__( "#annot0" );
union {
double f;
struct {
unsigned long m:52; // careful here with endianess
unsigned long x:11;
unsigned long s:1;
} u64;
uint64_t u64i;
} z,
magic0 = { .u64 = { 0, (1<<10)-1 + 52, 0 } },
magic1 = { .u64 = { 0, (1<<10)-1 + (52+12), 0 } },
magic2 = { .u64 = { 0, 2046, 0 } };
__asm__( "#annot1" );
if( u64 < (1UL << 52UL ) ) {
z.u64i = u64 + magic0.u64i;
z.f -= magic0.f;
} else {
z.u64i = ( u64 >> 12 ) + magic1.u64i;
z.f -= magic1.f;
}
__asm__( "#annot2" );
z.u64i = magic2.u64i - z.u64i;
return z.f;
}
Компиляция этого на ядре Intel 7 дает несколько инструкций (и ветку), но, разумеется, не размножается и не делит вообще. Если броски между int и double быстрые, это должно выполняться довольно быстро.
Я подозреваю, что float (всего 23 бита точности) потребует более 1 или 2 итераций Newton-Raphson, чтобы получить нужную вам точность, но я не выполнил математику...
