Быстрое 1 / Х деление (ответное) - PullRequest
22 голосов
/ 30 марта 2012

Есть ли способ улучшить обратную величину (деление 1 на X) относительно скорости, если точность не является критической?

Итак, мне нужно вычислить 1 / X. Есть ли обходной путь, чтобы я потерял точность, но делаю это быстрее?

Ответы [ 6 ]

7 голосов
/ 27 сентября 2016

Я считаю, что то, что он искал, - это более эффективный способ приблизить 1.0 / x вместо некоторого технического определения аппроксимации, в котором говорится, что вы могли бы использовать 1 как очень неточный ответ. Я также считаю, что это удовлетворяет это.

#ifdef __cplusplus
    #include <cstdint>
#else
    #include <stdint.h>
#endif

__inline__ double __attribute__((const)) reciprocal( double x ) {
    union {
        double dbl;
        #ifdef __cplusplus
            std::uint_least64_t ull;
        #else
            uint_least64_t ull;
        #endif
    } u;
    u.dbl = x;
    u.ull = ( 0xbfcdd6a18f6a6f52ULL - u.ull ) >> 1;
                                // pow( x, -0.5 )
    u.dbl *= u.dbl;             // pow( pow(x,-0.5), 2 ) = pow( x, -1 ) = 1.0 / x
    return u.dbl;
}
__inline__ float __attribute__((const)) reciprocal( float x ) {
    union {
        float single;
        #ifdef __cplusplus
            std::uint_least32_t uint;
        #else
            uint_least32_t uint;
        #endif
    } u;
    u.single = x;
    u.uint = ( 0xbe6eb3beU - u.uint ) >> 1;
                                // pow( x, -0.5 )
    u.single *= u.single;       // pow( pow(x,-0.5), 2 ) = pow( x, -1 ) = 1.0 / x
    return u.single;
}


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

Что касается тестирования производительности, инструкции по аппаратному обеспечению x 2 в сочетании с инструкциями по аппаратному вычитанию выполняются так же быстро, как и инструкция по оборудованию 1.0 / x на современных компьютерах (мои тесты были на Intel i7, но Я бы предположил, аналогичные результаты для других процессоров). Однако если бы этот алгоритм был внедрен в аппаратное обеспечение как новая инструкция по сборке, то увеличение скорости, вероятно, было бы достаточно хорошим, чтобы эта инструкция была довольно практичной.

Для получения дополнительной информации об этом методе, эта реализация основана на замечательном «быстром» алгоритме обратного квадратного корня .

Как Pharap обратил мое внимание, чтение неактивного свойства из объединения является неопределенным поведением, поэтому есть два возможных решения, которые я разработал из его полезных комментариев, чтобы избежать неопределенного поведения. Первое решение больше похоже на неприятную уловку, позволяющую обойти семантическую семантику, которая практически не лучше оригинального решения.

#ifdef __cplusplus
    #include <cstdint>
#else
    #include <stdint.h>
#endif
__inline__ double __attribute__((const)) reciprocal( double x ) {
    union {
        double dbl[2];
        #ifdef __cplusplus
            std::uint_least64_t ull[2];
        #else
            uint_least64_t ull[2];
        #endif
    } u;
    u.dbl[0] = x; // dbl is now the active property, so only dbl can be read now
    u.ull[1] = 0;//trick to set ull to the active property so that ull can be read
    u.ull][0] = ( 0xbfcdd6a18f6a6f52ULL - u.ull[0] ) >> 1;
    u.dbl[1] = 0; // now set dbl to the active property so that it can be read
    u.dbl[0] *= u.dbl[0];
    return u.dbl[0];
}
__inline__ float __attribute__((const)) reciprocal( float x ) {
    union {
        float flt[2];
        #ifdef __cplusplus
            std::uint_least32_t ull[2];
        #else
            uint_least32_t ull[2];
        #endif
    } u;
    u.flt[0] = x; // now flt is active
    u.uint[1] = 0; // set uint to be active for reading and writing
    u.uint[0] = ( 0xbe6eb3beU - u.uint[0] ) >> 1;
    u.flt[1] = 0; // set flt to be active for reading and writing
    u.flt[0] *= u.flt[0];
    return u.flt[0];
}

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

  1. , что байты имеют ширину 8 бит
  2. что байты являются наименьшей атомарной единицей на целевой машине.
  3. с двойным размером 8 байтов и с плавающей точкой 4 байта.

#ifdef __cplusplus
    #include <cstdint>
    #include <cstring>
    #define stdIntWithEightBits std::uint8_t
    #define stdIntSizeOfFloat std::uint32_t
    #define stdIntSizeOfDouble std::uint64_t
#else
    #include <stdint.h>
    #include <string.h>
    #define stdIntWithEightBits uint8_t
    #define stdIntSizeOfFloat uint32_t
    #define stdIntSizeOfDouble uint64_t
#endif

__inline__ double __attribute__((const)) reciprocal( double x ) {
    double byteIndexFloat = 1.1212798184631136e-308;//00 08 10 18 20 28 30 38 bits
    stdIntWithEightBits* byteIndexs = reinterpret_cast<stdIntWithEightBits*>(&byteIndexFloat);

    stdIntWithEightBits* inputBytes = reinterpret_cast<stdIntWithEightBits*>(&x);

    stdIntSizeOfDouble inputAsUll = (
        (inputBytes[0] << byteIndexs[0]) |
        (inputBytes[1] << byteIndexs[1]) |
        (inputBytes[2] << byteIndexs[2]) |
        (inputBytes[3] << byteIndexs[3]) |
        (inputBytes[4] << byteIndexs[4]) |
        (inputBytes[5] << byteIndexs[5]) |
        (inputBytes[6] << byteIndexs[6]) |
        (inputBytes[7] << byteIndexs[7])
    );
    inputAsUll = ( 0xbfcdd6a18f6a6f52ULL - inputAsUll ) >> 1;

    double outputDouble;

    const stdIntWithEightBits outputBytes[] = {
        inputAsUll >> byteIndexs[0],
        inputAsUll >> byteIndexs[1],
        inputAsUll >> byteIndexs[2],
        inputAsUll >> byteIndexs[3],
        inputAsUll >> byteIndexs[4],
        inputAsUll >> byteIndexs[5],
        inputAsUll >> byteIndexs[6],
        inputAsUll >> byteIndexs[7]
    };
    memcpy(&outputDouble, &outputBytes, 8);

    return outputDouble * outputDouble;
}

__inline__ float __attribute__((const)) reciprocal( float x ) {
    float byteIndexFloat = 7.40457e-40; // 0x00 08 10 18 bits
    stdIntWithEightBits* byteIndexs = reinterpret_cast<stdIntWithEightBits*>(&byteIndexFloat);

    stdIntWithEightBits* inputBytes = reinterpret_cast<stdIntWithEightBits*>(&x);

    stdIntSizeOfFloat inputAsInt = (
        (inputBytes[0] << byteIndexs[0]) |
        (inputBytes[1] << byteIndexs[1]) |
        (inputBytes[2] << byteIndexs[2]) |
        (inputBytes[3] << byteIndexs[3])
    );
    inputAsInt = ( 0xbe6eb3beU - inputAsInt ) >> 1;

    float outputFloat;

    const stdIntWithEightBits outputBytes[] = {
        inputAsInt >> byteIndexs[0],
        inputAsInt >> byteIndexs[1],
        inputAsInt >> byteIndexs[2],
        inputAsInt >> byteIndexs[3]
    };
    memcpy(&outputFloat, &outputBytes, 4);

    return outputFloat * outputFloat;
}

Отказ от ответственности: Наконец, обратите внимание, что я больше новичок в C ++. В связи с этим, я с широко распростертыми объятиями приветствую любые передовые практики, правильное форматирование или внесение ясности в смысл до конца как улучшения качества этого ответа для всех, кто его читает, так и расширения моих знаний о C ++ за все мои годы, чтобы прийти (если, конечно, я не попаду в автомобильную аварию завтра и умру).

5 голосов
/ 30 марта 2012

Во-первых, убедитесь, что это не случай преждевременной оптимизации.Вы знаете, что это ваше узкое место?

Как говорит Mystical, 1 / x можно вычислить очень быстро.Убедитесь, что вы не используете double тип данных для 1 или делителя.Поплавки намного быстрее.

Тем не менее, эталон, эталон, эталон.Не тратьте свое время, тратя часы на теорию чисел, просто чтобы обнаружить источник низкой производительности - доступ к IO.

3 голосов
/ 30 марта 2012

Прежде всего, если вы включите оптимизацию компилятора, компилятор, вероятно, оптимизирует вычисления, если это возможно (например, чтобы вытащить их из цикла). Чтобы увидеть эту оптимизацию, вам нужно собрать и запустить ее в режиме Release.

Деление может быть тяжелее, чем умножение (но комментатор указал, что взаимные вычисления столь же быстры, как умножение на современных процессорах, в этом случае это не правильно для вашего случая), поэтому, если у вас действительно появляется 1/X где-то внутри цикла (и более одного раза) вы можете помочь, кэшируя результат внутри цикла (float Y = 1.0f/X;), а затем используя Y. (Оптимизация компилятора может сделать это в любом случае.)

Кроме того, некоторые формулы могут быть переработаны для удаления деления или других неэффективных вычислений. Для этого вы можете опубликовать большее вычисление. Даже в этом случае саму программу или алгоритм иногда можно реструктурировать, чтобы исключить необходимость частого попадания на трудоемкие циклы.

Сколько точности можно пожертвовать? Если по случайности вам нужен только порядок величины, вы можете легко получить это, используя оператор модуля или побитовые операции.

Однако, в общем, нет способа ускорить деление. Если бы они были, компиляторы уже делали бы это.

1 голос
/ 08 июня 2019

Я тестировал эти методы на Arduino NANO на скорость и «точность».
Основной расчет был для установки переменных, Y = 15 000 000 и Z = 65 535
(в моем реальном случае Y - это константа, а Z может варьироваться от 65353 до 3000, поэтому полезный тест)
Время вычислений на Arduino было установлено путем установки минимального значения пин-кода, затем высокого значения, полученного при вычислении, а затем снова низкого уровня и сравнения времени с логическим анализатором. ЗА 100 ЦИКЛОВ. С переменными в виде целых чисел без знака: -

Y * Z takes 0.231 msec
Y / Z takes  3.867 msec.  
With variables as floats:-  
Y * Z takes  1.066 msec
Y / Z takes  4.113 msec.  
Basic Bench Mark  and ( 15,000,000/65535 = 228.885 via calculator.) 

Использование обратного алгоритма с плавающей точкой {Jack Giffin's}:

Y * reciprocal(Z)  takes  1.937msec  which is a good improvement, but accuracy less so 213.68.  

Использование {nimig18's} float inv_fast:

Y* inv_fast(Z)  takes  5.501 msec  accuracy 228.116  with single iteration  
Y* inv_fast(Z)  takes  7.895 msec  accuracy 228.883  with second iteration 

Использование Q_rsqrt из Википедии (на что указывает {Джек Гиффин})

Y * Q*rsqrt(Z) takes  6.104 msec  accuracy   228.116  with single iteration  
All entertaining but ultimately disappointing!
1 голос
/ 06 апреля 2017

Это должно быть сделано с несколькими предварительно развернутыми итерациями Ньютона, оцененными как полином Хорнера, который использует операции слияния-умножения с накоплением, которые большинство современных ЦП выполняют за один цикл Clk (каждый раз):

float inv_fast(float x) {
    union { float f; int i; } v;
    float w, sx;
    int m;

    sx = (x < 0) ? -1:1;
    x = sx * x;

    v.i = (int)(0x7EF127EA - *(uint32_t *)&x);
    w = x * v.f;

    // Efficient Iterative Approximation Improvement in horner polynomial form.
    v.f = v.f * (2 - w);     // Single iteration, Err = -3.36e-3 * 2^(-flr(log2(x)))
    // v.f = v.f * ( 4 + w * (-6 + w * (4 - w)));  // Second iteration, Err = -1.13e-5 * 2^(-flr(log2(x)))
    // v.f = v.f * (8 + w * (-28 + w * (56 + w * (-70 + w *(56 + w * (-28 + w * (8 - w)))))));  // Third Iteration, Err = +-6.8e-8 *  2^(-flr(log2(x)))

    return v.f * sx;
}

Fine Print: При приближении к 0 аппроксимация работает не так хорошо, поэтому либо вам, программисту, нужно проверить производительность, либо ограничить ввод от низкого до того, как прибегнуть к аппаратному делению. т.е. быть ответственным!

0 голосов
/ 16 января 2014

Самый быстрый способ, который я знаю, - это использовать SIMD-операции. http://msdn.microsoft.com/en-us/library/796k1tty(v=vs.90).aspx

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...