【问题标题】:Fast multiplication/division by 2 for floats and doubles (C/C++)浮点数和双精度数的快速乘法/除法 (C/C++)
【发布时间】:2011-12-04 23:03:35
【问题描述】:

在我正在编写的软件中,我正在对我的值进行数百万次乘法或除以 2(或 2 的幂)。我真的希望这些值是 int 以便我可以访问位移运算符

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

但是,我不能,我必须坚持双打。

我的问题是:由于有标准的双精度表示(符号、指数、尾数),有没有办法使用指数来获得 2 的幂的快速乘法/除法 ?

我什至可以假设位数将是固定的(该软件将在始终具有 64 位长双精度的机器上运行)

P.S : 是的,算法大多只做这些操作。这是瓶颈(它已经是多线程的了)。

编辑:还是我完全弄错了,聪明的编译器已经为我优化了?


临时结果(用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&lt;&lt;5);D = d*32; 同时运行(200 毫秒),而D = ldexp (d,5); 则慢得多(6000 毫秒)。我知道这是一个微基准测试,突然间,我的 RAM 爆炸了,因为 Chrome 突然要求在我每次运行 ldexp() 时计算 Pi,所以这个基准测试毫无价值.不过我还是会留着的。

另一方面,我在执行 reinterpret_cast&lt;uint64_t *&gt; 时遇到问题,因为存在 const 违规(似乎 volatile 关键字干扰)

【问题讨论】:

  • 不要仅仅因为它是多线程的就认为它是瓶颈。我们有一个多线程应用程序,我们发现它在许多地方都出现了与我们预期不同的瓶颈?您的分析准确度如何?
  • 与往常一样,对应用程序的分析永远不够。我的意思是,我使用了 CacheGrind,而且似乎我大部分时间都花在了一个主要用于乘法运算的函数上。它似乎。但我写道这是瓶颈,因为我对乘以 2 背后的理论思想更感兴趣,而不是“琐碎的考虑”(当然,我可以优化我的 SQL 请求,但老实说,我很确定它会与乘法相比毫无意义,而且大多数情况下,我不在乎 ^^)
  • 是的,我“知道” ^^ 但是有人说你可以用奇怪的东西混淆编译器(基本上,一些操作归结为*32,但编译器没有“看到”它)。这只是一行=)
  • 除了I used CacheGrind, and it seems I spend most of my time in a function that does mostly multiplications,我还能说什么呢?是的,当然,这个函数也做加法,并在栈上分配数据,但我想我对情况有一个很好的估计。
  • @Fezvez: 1) 你的循环需要展开。 2) CacheGrind 说你主要是在一些数学例程中?那是警报!在汇编语言级别单步执行代码,并确保它没有超出您的预期。它不应该调用anything。 3)多线程不会使代码更快。它充其量只是将其分散在更多的处理器上。 4) 如果性能是你真正关心的,learn this technique.

标签: c++ c optimization division multiplication


【解决方案1】:

从 c++17 开始,您还可以使用十六进制浮动文字。这样你就可以乘以 2 的高次幂。例如:

d *= 0x1p64;

d 乘以 2^64。我用它来实现我的快速整数算术转换为双精度。

【讨论】:

    【解决方案2】:

    根据您要相乘的内容,如果您有足够重复的数据,查找表可能会提供更好的性能,但会消耗内存。

    【讨论】:

    • 我不确定查找实际上是否比乘法更快。例如,了解查找表对于三角表的使用,但对于乘法?
    【解决方案3】:

    虽然对于 double 类型的 float 专门处理 2 的幂几乎没有/没有实际好处,但 double-double 类型有这种情况。 double-double 乘法和除法通常很复杂,但对于乘以和除以 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 重载了&lt;&lt;&gt;&gt;,所以它类似于整数。

    //x is a doubledouble type
    x << 2 // multiply x by four;
    x >> 3 // divide x by eight.
    

    【讨论】:

      【解决方案4】:

      您可以非常安全地假设 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 hack,因为这是加载到 FP 堆栈、执行一个运算并从中存储的最坏情况,但它是为什么 x87 的一个很好的例子不是追求任何性能的方法。相关。

      【讨论】:

      • 我试试看有没有效率!
      • 不知何故,我认为条件分支不会比 FP 乘法更快。
      • 我倾向于同意,但你总是会感到惊讶(好吧,我想感到惊讶!)
      • 请注意,使用普通的旧乘法允许编译器向量化循环 - 它一次处理两个双精度数。这可能就是它跑得更快的原因 - 让所有通过这种方式的人吸取教训! ;)
      • @SimonBuchan:两个 addls 比一个mulpd 慢 50%,不过 - 请记住,您正在为循环变量使用整数执行单元也会增加,因此使用单独的执行单元会带来一些好处。
      【解决方案5】:

      ldexp怎么样?

      任何半体面的编译器都会在您的平台上生成最佳代码。

      但正如@Clinton 指出的那样,简单地以“明显”的方式编写它应该也可以。对于现代编译器来说,乘以和除以 2 的幂是小菜一碟。

      直接修改浮点表示,除了不可移植之外,几乎肯定不会更快(而且很可能会更慢)。

      当然,您甚至不应该浪费时间去思考这个问题,除非您的分析工具告诉您这样做。但是听这个建议的人永远不会需要它,而那些需要它的人永远不会听。

      [更新]

      好的,所以我刚刚用 g++ 4.5.2 尝试了 ldexp。 cmath 标头将其内联为对 __builtin_ldexp 的调用,而后者又...

      ...发出对 libm ldexp 函数的调用。我原以为这个内置函数优化起来很简单,但我猜 GCC 开发人员从来没有考虑过它。

      因此,正如您所发现的,乘以 1 &lt;&lt; p 可能是您最好的选择。

      【讨论】:

      • VC 对大多数浮点运算都做同样的事情——我相信它可以尊重精度控制(_control87()_controlfp() 等...)。尝试摆弄浮动精度编译器开关...
      • ldexp 比常规 x*pow(2,exp) 源慢 6 倍:在 Intel Xeon 上进行基准测试
      【解决方案6】:

      这是高度应用特定的东西之一。它在某些情况下可能会有所帮助,而在其他情况下则无济于事。 (在绝大多数情况下,直接乘法仍然是最好的。)

      执行此操作的“直观”方法是将位提取为 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 时代,Integer SSE 和 FP SSE 使用相同的 ISA 寄存器(尽管它们仍然有单独的寄存器文件)。根据Agner Fog,在整数 SSE 和 FP SSE 执行单元之间移动数据会有 1 到 2 个周期的惩罚。所以成本比x87时代要好很多,但还是有的。

      总而言之,这将取决于您的管道上还有什么。但在大多数情况下,乘法仍然会更快。我以前也遇到过同样的问题,所以我是从第一手经验说的。

      现在有了只支持 FP 指令的 256 位 AVX 指令,玩这种把戏的动机就更小了。

      【讨论】:

      • 关于“只要……这将起作用”:根本不能保证它会起作用。它可能适用于给定的实现,但标准明确指出这(设置一种类型的联合并使用不同的类型来读回它)不是强制性的。由于我们已经进入了优化/非标准安全的行为,因此没有对此投反对票。
      • @paxdiablo:正确。我们已经远远超出了标准。我抛出这个例子只是为了展示这个想法。使用 SSE 寄存器更实际地完成。
      • @paxdiablo:当您知道机器的浮点表示时,它就可以工作。只要您知道这不会在 VAX 上运行(可能还有一些较旧的 IBM 大型机),您就知道那将是 IEEE 754。
      • SSE = 流式 SIMD 扩展,ISA = 指令集架构,FP = 浮点,AVX = 高级向量扩展
      • @paxdiablo:当然,别名规则尤其令人痛苦——我只是反对使用不可移植代码在某种程度上是不道德的(显然在有意义的地方)这种有点普遍的想法。跨度>
      【解决方案7】:

      乘以 2 可以用加法代替:x *= 2 等价于 x += x

      除以 2 可以替换为乘以 0.5。乘法通常比除法快得多。

      【讨论】:

      • 完全正确,但是一旦我想做x *= 33554432之类的事情,它就会变得难以处理
      • @Fezvez,浮点单元完成乘法的速度可能比您想出的任何优化都快。
      • 嗯,它只是在指数上加上 25,所以我猜我的审讯背后有一些意义 =)
      • @Fezvez,不要低估现代乘法指令的速度。如果您怀疑我,请测量并查看。
      • 我的意思是,这篇文章的重点是检查是否有更快的方法来进行一些非常具体的乘法运算。我并不是说存在更快的方法,但我觉得这样问是合理的。
      【解决方案8】:

      最快的方法可能是:

      x *= (1 << p);
      

      这种事情可以简单地通过调用机器指令将p 添加到指数来完成。告诉编译器改为使用掩码提取某些位并手动对其执行某些操作可能会使事情变得更慢,而不是更快。

      请记住,C/C++ 不是汇编语言。使用位移运算符不一定编译为位移汇编操作,使用乘法不一定编译为乘法。发生了各种奇怪而奇妙的事情,例如正在使用哪些寄存器以及可以同时运行哪些指令,这些我不够聪明,无法理解。但是您的编译器拥有多年的知识和经验以及强大的计算能力,在做出这些判断方面要好得多。

      ps 请记住,如果您的双精度数在数组或其他一些平面数据结构中,您的编译器可能非常聪明,并且使用 SSE 来同时多个 2 甚至 4 个双精度数.但是,进行大量位移可能会使您的编译器感到困惑并阻止这种优化。

      【讨论】:

      • 我不知道任何具有“将 p 添加到指数的机器指令”的架构。但也许应该有一个。
      • @masterxilo: x87 fscale 正是这样做的,但仅适用于 double 输入,而不是整数。它确实是x * (1&lt;&lt;trunc(y))x_exponent += trunc(y)。但它并不快:比fmul 慢得多,例如经典 P5 Pentium 上的 20 到 32 个时钟周期与 fmul 上的 3 个时钟周期,并且在现代 x86 上也好不了多少。 (agner.org/optimize)。所以fmul0.5 的常量要好得多
      • 但是正如这个和 Mysticial 的回答所指出的那样,现代 SIMD ISA 通常对 FP 和向量整数使用相同的寄存器。可以使用paddd xmm0, xmm1 之类的指令在mulps xmm0, xmm2 之类的指令之间对FP 位模式进行整数加法。当然,这并不能处理指数超出范围的情况,FP 乘法会给你无穷大,或者正确的下溢到次正规。 (或将偏置指数从 0(次正规)包装为全一(NaN 或无穷大,取决于有效数字。)
      【解决方案9】:

      该算法还需要哪些其他操作?您也许可以将浮点数分解为 int 对(符号/尾数和幅度),进行处理,并在最后重构它们。

      【讨论】:

      • 呃,好吧,我在这里和那里做了一些事情(矩阵乘法等......)我想这可能是一个好主意,但我认为这将是一个工作量(重新定义+-*,...)
      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2010-09-14
      • 1970-01-01
      • 1970-01-01
      • 2011-10-24
      • 1970-01-01
      • 2015-04-30
      • 2010-12-19
      相关资源
      最近更新 更多