Эффективное вычисление 2**64 / делитель с помощью быстрой обратной величины с плавающей точкой
В настоящее время я изучаю способы использования быстрой взаимной возможности с плавающей запятой одинарной точности различных современных процессоров для вычисления начального приближения для 64-разрядного целочисленного деления без знака на основе итераций Ньютона-Рафсона с фиксированной запятой. Это требует вычисления 264 / делитель, с максимально возможной точностью, где начальное приближение должно быть меньше или равно математическому результату, основываясь на требованиях следующих итераций с фиксированной точкой. Это означает, что это вычисление должно давать занижение. В настоящее время у меня есть следующий код, который хорошо работает на основе обширного тестирования:
#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 и наоборот, через встроенные 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-битный шаг Ньютона-Рафсона с фиксированной точкой.
Есть ли способы повысить эффективность моего нынешнего подхода? Представляют интерес портативные и полупортативные способы с использованием встроенных функций для конкретных платформ (в частности, для 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 бита точности.
Результат в этом коде будет двойным, но преобразование в число с плавающей запятой тривиально и может быть выполнено по битам или напрямую, в зависимости от аппаратного обеспечения.
После этого вы захотите выполнить итерацию Ньютона-Рафсона.
Большая часть этого кода просто конвертируется в магические числа.
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 Core 7 дает ряд инструкций (и ветвь), но, конечно, не умножает или делит вообще. Если приведение между int и double происходит быстро, оно должно выполняться довольно быстро.
Я подозреваю, что float (только с 23 битами точности) потребует более 1 или 2 итераций Ньютона-Рафсона, чтобы получить желаемую точность, но я не сделал математику...