gcc 当两个std::complex相乘时,为什么调用__muldc3?

zf9nrax1  于 2023-06-23  发布在  其他
关注(0)|答案(4)|浏览(109)

我天真地假设,复数乘法将被编译器内联,例如,对于这个函数:

#include <complex>

void mult(std::complex<double> &a, std::complex<double> &b){
    a*=b;
}

然而,当由gcc(使用-O2)编译时,resulting assembler令人惊讶(至少对我来说):

mult(std::complex<double>&, std::complex<double>&):
        pushq   %rbx
        movsd   8(%rdi), %xmm3
        movsd   (%rdi), %xmm2
        movq    %rdi, %rbx
        movsd   8(%rsi), %xmm1
        movsd   (%rsi), %xmm0
        call    __muldc3
        movsd   %xmm0, (%rbx)
        movsd   %xmm1, 8(%rbx)
        popq    %rbx
        ret

有一个对这个函数__multdc3的调用,它以某种方式取代了对operator*=的调用(它的错误名称将是_ZNSt7complexIdEmLIdEERS0_RKS_IT_E,复数将被每个引用传递)。
然而,在operator*=的实现中似乎没有什么特别的东西可以解释这个魔术:

// 26.2.5/13
  // XXX: This is a grammar school implementation.
  template<typename _Tp>
    template<typename _Up>
    complex<_Tp>&
    complex<_Tp>::operator*=(const complex<_Up>& __z)
    {
      const _Tp __r = _M_real * __z.real() - _M_imag * __z.imag();
      _M_imag = _M_real * __z.imag() + _M_imag * __z.real();
      _M_real = __r;
      return *this;
}

我一定错过了什么,所以我的问题是:产生汇编程序的原因是什么?

7lrncoxx

7lrncoxx1#

你应该注意到,严格地说,用公式实现复杂的浮点乘法是“错误的”

(a+i*b)*(c + i*d) = a*c - b*d + i*(b*c + a*d)

我写错了引号,因为C++标准实际上并不需要正确的实现。C在某个附录中明确说明了这一点。
简单的实现在输入中使用Inf和/或NaN时不会给予正确的结果。
考虑(Inf + 0*i)*(Inf + 0*i):显然,对于一致的行为,结果应该与真实的浮点相同,即分别为Inf(Inf + 0*i)。然而,上面的公式给出Inf + i*NaN
因此,我可以想象__muldc3是一个更长的函数,可以正确处理这些情况。
当调用与-ffast-math一起消失时,这很可能是解释。

s1ag04yj

s1ag04yj2#

至少对于纯C,gcc遵循有点疯狂的C99附件f(?)复数乘/除的规则; C++也可以。尝试-fast-math或-fcx-fortran-rules。

voase2hg

voase2hg3#

Andreas H.给出了方向,在这个答案中,我只是试着把这些点联系起来。
首先,我对operator*=的定义是错误的:对于浮点类型有专门化,例如对于double,这个运算符定义为:

template<typename _Tp>
    complex&
    operator*=(const complex<_Tp>& __z)
{
  _ComplexT __t;
  __real__ __t = __z.real();
  __imag__ __t = __z.imag();
  _M_value *= __t;
  return *this;
}

_ComplexT定义为

typedef __complex__ double _ComplexT;

_M_value定义为

private:
   _ComplexT _M_value;

这意味着,school-formula用于std::complex<int>等类型,但对于浮点类型,它归结为C乘法(__complex__是gcc扩展),这是标准的一部分,因为C99,附录G(最重要的部分在答案的末尾)。
其中一个问题是:

(0.0+1.0*j)*(inf+inf*j) = (0.0*inf-1*inf)+(0.0*inf+1.0*inf)j
                        =  nan + nan*j

由于约定(G.3.1),第一个因子是非零的,有限的,而第二个因子是无限的,因此由于(G.5.1.4),结果必须是有限的,而nan+nan*j不是这种情况。
__multdc3中的计算策略很简单:使用school公式,如果这不起作用(两者都是nan),将使用更复杂的方法-这只是太复杂了,无法内联。(顺便说一下,clang会内联学校公式,如果这个公式不起作用,就会调用__multdc3)。
另一个问题可能是一个产品上溢/下溢,但两个产品的总和不会。这可能包含在“可能引发虚假浮点异常”的表述中,但我不是100%确定。
G.3.1节规定:
一个复数或虚数至少有一个无穷部分被认为是无穷大(即使它的另一部分是NaN)。
复数或虚数是一个有限数,如果它的每个部分都是有限数(既不是无限的也不是NaN)。
如果复值或虚数的每个部分都是零,那么它就是零。
第G.5节涉及二元运算符和状态:
对于大多数操作数类型,具有虚操作数或复操作数的二元运算符的结果的值完全由通常的数学公式参考真实的算术确定。对于某些操作数类型,通常的数学公式是有问题的,因为它对无穷大的处理和由于不适当的上溢或下溢;在这些情况下,结果满足某些性质(在G.5.1中规定),但不完全确定。
第G.5.1节涉及乘法算子和状态:
1.* 和/运算符满足所有真实的、虚数和复数操作数的下列无穷大属性:

  • 如果一个操作数是无穷大并且另一个操作数是非零有限数或无穷大,则运算符的结果是无穷大;
    ...
    1.如果 * 运算符的两个操作数都是复数,或者如果/运算符的第二个操作数是复数,则该运算符在适合计算结果部分时引发浮点异常,并且可能引发虚假浮点异常。
vc9ivgsu

vc9ivgsu4#

这可能是因为软件浮点运算在默认情况下是开启的。根据gcc手册页:

-msoft-float
This option ignored; it is provided for compatibility purposes only.    
software floating point code is emitted by default,
and this default can overridden by FPX options; mspfp, mspfp-compact, 
or mspfp-fast for single precision, and mdpfp, mdpfp-compact, or mdpfp-fast 
for double precision.

尝试使用-mdpfp或任何其他替代方法进行编译。

相关问题