Быстрое умножение / деление на 2 для чисел с плавающей запятой и двойных чисел (C/C++)
В программном обеспечении, которое я пишу, я делаю миллионы умножения или деления на 2 (или степени 2) моих значений. Мне бы очень хотелось, чтобы эти значения были int чтобы я мог получить доступ к операторам bithift
int a = 1;
int b = a<<24
Тем не менее, я не могу, и я должен придерживаться двойников.
У меня такой вопрос: как существует стандартное представление двойных чисел (знак, экспонента, мантисса), есть ли способ поиграть с показателем степени, чтобы получить быстрое умножение / деление на степень 2?
Я даже могу предположить, что число бит будет фиксированным (программное обеспечение будет работать на машинах, которые всегда будут иметь 64-битные двойные числа)
PS: И да, алгоритм в основном выполняет только эти операции. Это узкое место (оно уже многопоточное).
Изменить: Или я полностью ошибаюсь, и умные компиляторы уже оптимизируют вещи для меня?
Временные результаты (с помощью Qt для измерения времени, перебор, но мне все равно):
#include <QtCore/QCoreApplication>
#include <QtCore/QElapsedTimer>
#include <QtCore/QDebug>
#include <iostream>
#include <math.h>
using namespace std;
int main(int argc, char *argv[])
{
QCoreApplication a(argc, argv);
while(true)
{
QElapsedTimer timer;
timer.start();
int n=100000000;
volatile double d=12.4;
volatile double D;
for(unsigned int i=0; i<n; ++i)
{
//D = d*32; // 200 ms
//D = d*(1<<5); // 200 ms
D = ldexp (d,5); // 6000 ms
}
qDebug() << "The operation took" << timer.elapsed() << "milliseconds";
}
return a.exec();
}
Прогоны предполагают, что D = d*(1<<5); а также D = d*32; работать в одно и то же время (200 мс), тогда как D = ldexp (d,5); намного медленнее (6000 мс). Я знаю, что это микро-эталонный тест, и что внезапно моя оперативная память взорвалась, потому что Chrome внезапно просил вычислять Pi в моей спине каждый раз, когда я запускаю ldexp()так что этот тест ничего не стоит. Но я все равно буду держать это.
С другой стороны, у меня возникли проблемы с выполнением reinterpret_cast<uint64_t *> потому что есть const нарушение (кажется volatile ключевое слово мешает)
8 ответов
Вы можете довольно смело предположить форматирование IEEE 754, детали которого могут стать довольно грубыми (особенно, когда вы попадаете в субнормалы). В общих случаях, однако, это должно работать:
const int DOUBLE_EXP_SHIFT = 52;
const unsigned long long DOUBLE_MANT_MASK = (1ull << DOUBLE_EXP_SHIFT) - 1ull;
const unsigned long long DOUBLE_EXP_MASK = ((1ull << 63) - 1) & ~DOUBLE_MANT_MASK;
void unsafe_shl(double* d, int shift) {
unsigned long long* i = (unsigned long long*)d;
if ((*i & DOUBLE_EXP_MASK) && ((*i & DOUBLE_EXP_MASK) != DOUBLE_EXP_MASK)) {
*i += (unsigned long long)shift << DOUBLE_EXP_SHIFT;
} else if (*i) {
*d *= (1 << shift);
}
}
РЕДАКТИРОВАТЬ: После некоторого времени, этот метод странно медленнее, чем метод double на моем компиляторе и машине, даже обрезается до минимально исполняемого кода:
double ds[0x1000];
for (int i = 0; i != 0x1000; i++)
ds[i] = 1.2;
clock_t t = clock();
for (int j = 0; j != 1000000; j++)
for (int i = 0; i != 0x1000; i++)
#if DOUBLE_SHIFT
ds[i] *= 1 << 4;
#else
((unsigned int*)&ds[i])[1] += 4 << 20;
#endif
clock_t e = clock();
printf("%g\n", (float)(e - t) / CLOCKS_PER_SEC);
В DOUBLE_SHIFT завершается через 1,6 секунды, с внутренней петлей
movupd xmm0,xmmword ptr [ecx]
lea ecx,[ecx+10h]
mulpd xmm0,xmm1
movupd xmmword ptr [ecx-10h],xmm0
По сравнению с 2,4 секундами в противном случае с внутренней петлей:
add dword ptr [ecx],400000h
lea ecx, [ecx+8]
Действительно неожиданно!
РЕДАКТИРОВАТЬ 2: Тайна раскрыта! Одним из изменений для VC11 является то, что теперь он всегда векторизует циклы с плавающей запятой, эффективно форсируя /arch:SSE2, хотя VC10, даже с /arch:SSE2, все еще хуже с 3,0 секундами с внутренним циклом:
movsd xmm1,mmword ptr [esp+eax*8+38h]
mulsd xmm1,xmm0
movsd mmword ptr [esp+eax*8+38h],xmm1
inc eax
VC10 без /arch:SSE2 (даже с /arch:SSE) составляет 5,3 секунды... с 1/100 итерациями!!, внутренний цикл:
fld qword ptr [esp+eax*8+38h]
inc eax
fmul st,st(1)
fstp qword ptr [esp+eax*8+30h]
Я знал, что стек x87 FP ужасен, но в 500 раз хуже, это просто смешно. Вы, вероятно, не увидите такого рода ускорения, например, матричных операций в SSE или int-хаков, поскольку это наихудший случай загрузки в стек FP, выполнения одной операции и сохранения из нее, но это хороший пример того, почему x87 это не способ пойти на что-либо перф. связанные с.
Это одна из тех вещей, которые сильно зависят от приложения. Это может помочь в некоторых случаях, а не в других. (В подавляющем большинстве случаев прямое умножение все еще лучше.)
"Интуитивно понятный" способ сделать это - просто извлечь биты в 64-битное целое и добавить значение смещения непосредственно в экспоненту. (это будет работать до тех пор, пока вы не нажмете NAN или INF)
Так что-то вроде этого:
union{
uint64 i;
double f;
};
f = 123.;
i += 0x0010000000000000ull;
// Check for zero. And if it matters, denormals as well.
Обратите внимание, что этот код не является C-совместимым в любом случае, и показан только для иллюстрации идеи. Любая попытка реализовать это должна быть сделана непосредственно в ассемблере или в SSE.
Однако в большинстве случаев затраты на перемещение данных из блока FP в целочисленный блок (и обратно) будут стоить гораздо дороже, чем просто выполнить прямое умножение. Это особенно актуально для эпохи, предшествующей SSE, когда значение необходимо сохранить из FPU x87 в память и затем прочитать обратно в целочисленные регистры.
В эпоху SSE Integer SSE и FP SSE используют одни и те же регистры ISA (хотя они по-прежнему имеют отдельные файлы регистров). Согласно Agner Fog, для перемещения данных между исполнительными блоками Integer SSE и FP SSE предусмотрены штрафы от 1 до 2 циклов. Таким образом, стоимость намного лучше, чем эпоха x87, но она все еще там.
В целом, это будет зависеть от того, что еще у вас есть на вашем конвейере. Но в большинстве случаев умножение все равно будет быстрее. Раньше я сталкивался с той же самой проблемой, поэтому говорю из личного опыта.
Теперь с 256-битными инструкциями AVX, которые поддерживают только инструкции FP, стимулов для подобных трюков еще меньше.
Как насчет ldexp?
Любой полуприличный компилятор сгенерирует оптимальный код на вашей платформе.
Но, как указывает @Clinton, простое написание "очевидным" способом должно быть так же хорошо. Умножение и деление на степени двух - это детская игра для современного компилятора.
Непосредственное манипулирование представлением с плавающей запятой, помимо непереносимости, почти наверняка будет не быстрее (и вполне может быть медленнее).
И, конечно, вы не должны тратить время даже на размышления над этим вопросом, если ваш инструмент профилирования не скажет вам об этом. Но людям, которые слушают этот совет, он никогда не понадобится, а тем, кто в нем нуждается, никогда не прислушаются.
[Обновить]
Итак, я только что попробовал ldexp с g++ 4.5.2. cmath заголовок встраивает это как вызов __builtin_ldexp, которые в свою очередь...
... испускает вызов в libm ldexp функция. Я бы подумал, что это встроенное средство будет тривиальным для оптимизации, но я думаю, что разработчики GCC никогда не удосужились об этом.
Итак, умножение на 1 << p это, вероятно, ваш лучший выбор, как вы обнаружили.
Самый быстрый способ сделать это, вероятно:
x *= (1 << p);
Такого рода вещи можно просто сделать, вызвав машинную инструкцию для добавления p экспоненте. Указание компилятору извлекать некоторые биты с помощью маски и делать что-то вручную, вероятно, замедлит, а не ускорит процесс.
Помните, C/C++ не является языком ассемблера. Использование оператора битового сдвига не обязательно компилируется в операцию сборки битового сдвига, и при этом использование умножения не обязательно компилируется в умножение. Происходит всякое странное и удивительное, например, какие регистры используются и какие инструкции можно запускать одновременно, что я недостаточно умен, чтобы понять. Но ваш компилятор, имеющий много человеко-лет знаний и опыта и много вычислительных мощностей, намного лучше делает эти суждения.
ps Имейте в виду, что если ваши двойники находятся в массиве или какой-либо другой плоской структуре данных, ваш компилятор может быть очень умным и использовать SSE для нескольких 2 или даже 4 двойных одновременно. Тем не менее, выполнение большого сдвига битов может запутать ваш компилятор и помешать этой оптимизации.
Начиная с C++17 вы также можете использовать шестнадцатеричные плавающие литералы. Таким образом, вы можете умножить на более высокую степень 2. Например:
d *= 0x1p64;
умножится dна 2^64. Я использую его для реализации моей быстрой целочисленной арифметики при преобразовании в double.
Хотя практическая польза для обработки степеней двух специально для чисел с плавающей запятой невелика / не имеет практического смысла, есть основания для двойных-двойных типов. Двойное двойное умножение и деление сложны в общем, но тривиально для умножения и деления на степень двух.
Например, для
typedef struct {double hi; double lo;} doubledouble;
doubledouble x;
x.hi*=2, x.lo*=2; //multiply x by 2
x.hi/=2, x.lo/=2; //divide x by 2
На самом деле я перегружен << а также >> за doubledouble так что это аналог целых чисел.
//x is a doubledouble type
x << 2 // multiply x by four;
x >> 3 // divide x by eight.
Умножение на 2 можно заменить на дополнение: x *= 2 эквивалентно x += x,
Деление на 2 можно заменить умножением на 0,5. Умножение обычно значительно быстрее деления.
Какие еще операции требует этот алгоритм? Вы можете разбить ваши поплавки на пары int (знак / мантисса и величина), выполнить обработку и воссоздать их в конце.
В зависимости от того, что вы умножаете, если у вас есть достаточно повторяющиеся данные, справочная таблица может обеспечить более высокую производительность за счет памяти.