Ответ 1
Я считаю, что то, что он искал, - это более эффективный способ приближения 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];
}
Второе возможное решение гораздо более приемлемо, потому что оно полностью избавляет от профсоюзов. Однако это решение будет намного медленнее, если оно не будет должным образом оптимизировано компилятором. Но, с другой стороны, решение ниже будет полностью независимым от порядка следования байтов:
- что байты имеют ширину 8 бит
- эти байты являются наименьшей атомарной единицей на целевой машине.
- двойники имеют ширину 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++ для всех моих долгие годы (если, конечно, я не попаду в автомобильную аварию завтра и умру).