对浮点数和双精度快速乘法/除2（C / C ++）
在我正在编写的软件中，我正在进行数百万次乘法或除以2（或2的幂）的值。 我真的希望这些值为int以便我可以访问bitshift运算符
int a = 1; int b = a<<24
但是，我不能，而且我必须坚持双打。
我的问题是： 由于存在双精度（符号，指数，尾数）的标准表示，有没有办法使用指数来获得2的幂的快速乘法/除法 ？
我甚至可以假设位数将被修复（该软件将在总是具有64位长的双倍的机器上工作）
PS：是的，该算法主要只执行这些操作。 这是瓶颈（它已经是multithreading的）。
编辑：或者我完全错了，聪明的编译器已经为我优化了一些东西？
临时结果（用Qt测量时间，矫枉过正，但我不在乎）：
#include #include #include #include #include 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 ms），而D = ldexp (d,5); 慢得多（6000毫秒）。 我知道这是一个微型基准测试，突然之间，我的RAM已经爆炸，因为我每次运行ldexp()时都会突然要求我在后面计算Pi，所以这个基准测试没什么价值。 但我会保留它。
另一方面，我在执行reinterpret_cast遇到问题，因为存在const违规（似乎volatile关键字干扰）
你可以非常安全地假设IEEE 754格式化，其细节可以得到相当的gnarley（特别是当你进入subnormals时）。 但是，在常见情况下，这应该有效：
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的一个变化是它现在总是向量化浮点循环，有效地强制/拱形：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
没有 / arch的VC10：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 hacks，因为这是加载到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之前的时代尤其如此，需要将值从x87 FPU存储到存储器中，然后再读回整数寄存器。
在SSE时代，整数SSE和FP SSE使用相同的ISA寄存器（尽管它们仍然具有单独的寄存器文件）。 根据Agner Fog的说法 ，在整数SSE和FP SSE执行单元之间移动数据会有1到2个周期的惩罚。 所以成本比x87时代要好得多，但它仍然存在。
总而言之，它取决于你在管道上还有什么。 但在大多数情况下，乘法仍然会更快。 我之前遇到过这个完全相同的问题，所以我从第一手经验开始说话。
现在使用仅支持FP指令的256位AVX指令，更不用说像这样的技巧了。
ldexp怎么样 ？
任何半合适的编译器都会在您的平台上生成最佳代码。
但正如@Clinton所指出的那样，只需以“明显”的方式编写它也应该做得很好。 乘以和除以2的幂是现代编译器的儿童游戏。
除了不可移植之外，直接修改浮点表示几乎肯定不会更快（并且可能更慢）。
当然，除非你的分析工具告诉你，否则你不应该浪费时间思考这个问题。 但听取这种建议的人永远不会需要它，而那些需要它的人永远不会听。
[更新]
好的，所以我只是尝试使用g ++ 4.5.2进行ldexp。 cmath标题将其内联为对__builtin_ldexp的调用，而__builtin_ldexp依次调用…
…发出对libm ldexp函数的调用。 我本以为这个内置版本很容易进行优化，但我想GCC开发人员从未接触过它。
因此，正如您所发现的那样，乘以1 << p可能是您最好的选择。
最快的方法可能是：
x *= (1 << p);
这种事情可以简单地通过调用机器指令来将p添加到指数来完成。 告诉编译器改为使用掩码提取一些位并手动执行某些操作可能会使事情变慢，而不是更快。
请记住，C / C ++不是汇编语言。 使用bitshift运算符不一定编译为bitshift汇编操作，使用乘法不一定编译为乘法。 有各种奇怪和奇妙的事情正在发生，例如正在使用的寄存器以及可以同时运行哪些指令，我不够聪明。 但是你的编译器具有许多人的知识和经验以及大量的计算能力，在做出这些判断方面要好得多。
ps请记住，如果你的双打是一个数组或其他平面数据结构，你的编译器可能非常聪明，并且同时使用SSE多次2，甚至4个双打。 但是，进行大量的位移可能会使编译器混淆并阻止这种优化。
该算法还需要哪些其他操作？ 您可以将浮点数分成int对（符号/尾数和幅度），进行处理，最后重构它们。
乘以2可以用加法代替： x *= 2相当于x += x 。
除以2可以乘以0.5来代替。 乘法通常比除法快得多。
虽然对于两种类型的浮子来说处理两种力量几乎没有实际好处，但对于双重类型来说存在这种情况。 双倍乘法和除法一般是复杂的，但乘法和除以2的幂是微不足道的。
例如
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.
根据您所乘的数据，如果您的数据足够重复，则查找表可能会以内存为代价提供更好的性能。