The answer of Andreas H. gives the direction, in this answer I just try to connect the dots.
First, I was wrong about the definition of the operator*=
: there is a specialization for floating types, for example for doubles this operator defined as:
template<typename _Tp>
complex&
operator*=(const complex<_Tp>& __z)
{
_ComplexT __t;
__real__ __t = __z.real();
__imag__ __t = __z.imag();
_M_value *= __t;
return *this;
}
With _ComplexT
defined as
typedef __complex__ double _ComplexT;
and _M_value
defined as
private:
_ComplexT _M_value;
That means, the school-formula is used for such types as std::complex<int>
but for the floating types it boils down to C-multiplication (__complex__
is a gcc-extension), which is a part of the standard since C99, appendix G (most important parts are at the end of the answer).
One of the problematic cases would be:
(0.0+1.0*j)*(inf+inf*j) = (0.0*inf-1*inf)+(0.0*inf+1.0*inf)j
= nan + nan*j
Due to the convention (G.3.1) the first factor is a nonzero, finite and the second is infinite, thus due to (G.5.1.4) the result must be finite, which is not the case for nan+nan*j
.
The strategy of the calculation in __multdc3
is simple: use the school formula, and if this doesn't work (both are nan
) a more complicated approach will be used - which is just too complex to be inlined. (By the way, clang inlines the school formula and calls __multdc3
if this formula didn't work).
Another problem could be that one product overflow/underflow but the sum of two products doesn't. This is probably covered by the formulation "may raise spurious floating-point exceptions", but I'm not 100% sure.
The section G.3.1 states:
A complex or imaginary value with at least one infinite part is regarded as an infinity (even if its other part is a NaN).
A complex or imaginary value is a finite number if each of its parts is a finite number (neither infinite nor NaN).
A complex or imaginary value is a zero if each of its parts is a zero.
Section G.5 is concerned with binary operators and states:
For most operand types, the value of the result of a binary operator
with an imaginary or complex operand is completely determined, with
reference to real arithmetic, by the usual mathematical formula. For
some operand types, the usual mathematical formula is
problematic because of its treatment of infinities and because
of undue overflow or underflow; in these cases the result satisfies
certain properties (specified in G.5.1), but is not completely
determined.
The section G.5.1 is concerned with multiplicative operators and states:
The * and / operators satisfy the following infinity properties for
all real, imaginary,and complex operands:
—if one operand is an infinity and the other operand is a
nonzero finite number or an infinity,then the result of the
operator is an infinity;
...
- If both operands of the * operator are complex or if the second operand of the / operator is complex, the operator raises floating-point exceptions if appropriate for the calculation of the parts of the result, and may raise spurious floating-point exceptions.