Быстрое умножение / деление на 2 для чисел с плавающей точкой и двойных чисел (C / C ++)

Question

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

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 мешает)

Answer 1

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

«Интуитивно понятный» способ сделать это - просто извлечь биты в 64-разрядное целое число и напрямую добавить значение сдвига.в экспоненте.(это будет работать до тех пор, пока вы не нажмете NAN или INF)

Так что-то вроде этого:

union{
    uint64 i;
    double f;
};

f = 123.;
i += 0x0010000000000000ull;

//  Check for zero. And if it matters, denormals as well.

Обратите внимание, что этот код не совместим с C в любомКстати, и показано только для иллюстрации идеи.Любая попытка реализовать это должна быть сделана непосредственно в ассемблере или в SSE.

Однако в большинстве происходит перегрузка данных из FPunit to integer unit (и обратно) будет стоить намного дороже, чем просто делать умножение сразу.Это особенно актуально для эпохи, предшествующей SSE, когда значение необходимо сохранить из FPU x87 в память и затем прочитать обратно в целочисленные регистры.

В эпоху SSE Integer SSE и FP SSE используютте же самые регистры ISA (хотя у них все еще есть отдельные файлы регистров).Согласно Agner Fog , для перемещения данных между исполнительными блоками Integer SSE и FP SSE существует штраф от 1 до 2 циклов.Таким образом, стоимость намного лучше, чем в эпоху x87, но она все еще там.

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

Теперь с 256-битными инструкциями AVX, которые поддерживают только инструкции FP, стимулов для таких трюков еще меньшеэто.

Answer 2

Как насчет ldexp ?

Любой полуприличный компилятор сгенерирует оптимальный код на вашей платформе.

Но, как указывает @Clinton, просто напишите его в«очевидный» способ должен делать то же самое.Умножение и деление на степени двух - это детская игра для современного компилятора.

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

И, конечно, вы не должны тратить время, даже размышляя над этим вопросом, если ваш инструмент профилирования не скажет вам об этом.Но людям, которые слушают этот совет, он никогда не понадобится, а тем, кто в нем нуждается, он никогда не прислушается.

[update]

ОК, поэтому я просто попытался выполнить ldexp с помощью g ++ 4.50,2.Заголовок cmath указывает на это как вызов __builtin_ldexp, который в свою очередь ...

... испускает вызов функции libm ldexp.Я бы подумал, что эта встроенная функция будет тривиальной для оптимизации, но я думаю, что разработчики GCC никогда не удосужились к ней.

Так что, как вы обнаружили, умножение на 1 << p, вероятно, является вашим лучшим выбором.

Answer 3

Вы можете довольно смело предположить форматирование 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 это не способ пойти на что-либо перф. связаны между собой.

Answer 4

Самый быстрый способ сделать это, вероятно:

x *= (1 << p);

Такого рода вещи можно просто сделать, вызвав машинную инструкцию для добавления p к показателю степени.Указание компилятору извлекать некоторые биты с маской и делать что-то вручную с ней, вероятно, сделает вещи медленнее, а не быстрее.

Помните, C / C ++ не является языком ассемблера.Использование оператора битового сдвига не обязательно компилируется в операцию сборки битового сдвига, и при этом использование умножения не обязательно компилируется в умножение.Происходит всякое странное и удивительное, например, какие регистры используются и какие инструкции можно выполнять одновременно, что я недостаточно умен, чтобы понять.Но ваш компилятор, обладающий многолетними знаниями и опытом и большим вычислительным потенциалом, намного лучше делает эти суждения.

ps Имейте в виду, если ваши двойники находятся в массивеили какая-то другая плоская структура данных, ваш компилятор может быть действительно умным и использовать SSE для кратных 2 или даже 4 двойных одновременно.Тем не менее, выполнение большого сдвига битов может запутать ваш компилятор и предотвратить эту оптимизацию.

Answer 5

Хотя практическая обработка двух степеней, особенно для чисел с плавающей запятой, имеет мало практических преимуществ или вообще не дает их, для этого есть случаи для double-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.

Answer 6

Умножение на 2 можно заменить на сложение: x *= 2 эквивалентно x += x.

Деление на 2 можно заменить умножением на 0,5.Умножение обычно значительно быстрее деления.

Answer 7

Какие еще операции требует этот алгоритм?Вы можете разбить ваши поплавки на пары int (знак / мантисса и величина), выполнить обработку и восстановить их в конце.

Answer 8

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