Подтвердить что ты не робот

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

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

int a = 1;
int b = a<<24

Однако я не могу, и я должен придерживаться двойников.

Мой вопрос: , поскольку существует стандартное представление удвоений (знак, экспонента, мантисса), есть ли способ играть с показателем, чтобы получить быстрые умножения/деления на мощность 2

Я даже могу предположить, что количество бит будет исправлено (программное обеспечение будет работать на машинах, которые всегда будут иметь удвоение на 64 бита)

P.S: И да, алгоритм в основном выполняет только эти операции. Это узкое место (оно уже многопоточное).

Изменить: Или я полностью ошибаюсь, и умные компиляторы уже оптимизируют для меня вещи?


Временные результаты (с Qt для измерения времени, overkill, но мне все равно):

#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 вмешивается)

4b9b3361

Ответ 1

Вы можете с уверенностью предполагать форматирование IEEE 754, детали которого могут получить довольно gnarley (особенно, когда вы попадаете в субнормальные файлы). Однако в общих случаях это должно работать:

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);
    }
} 

EDIT: после некоторого времени этот метод странно медленнее, чем двойной метод моего компилятора и машины, даже лишенный минимального кода:

    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 раз хуже, это просто смешно. Вероятно, вы не увидите таких преобразований ускорения, то есть матрицы ops для SSE или int hacks, так как это худший случай загрузки в стек FP, выполнение одного op и сохранение из него, но это хороший пример для почему x87 это не способ пойти на что-то перф. связаны между собой.

Ответ 2

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

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

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

union{
    uint64 i;
    double f;
};

f = 123.;
i += 0x0010000000000000ull;

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

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

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

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

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

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

Ответ 3

Как насчет ldexp?

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

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

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

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

[обновление]

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

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

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

Ответ 4

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

x *= (1 << p);

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

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

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

Ответ 5

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

Ответ 6

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

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

Ответ 7

Несмотря на то, что практического преимущества для рассмотрения полномочий двух специально для float двойных типов мало, нет для этого случая для 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.

Ответ 8

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