I'm trying to write std::complex as an HLSL library. For this, I set out to implement the most basic functionality, starting with arithmetic operators. For finite numbers, all is well as expected. But things get weird at infinity. I won't be touching any HLSL here, it's all C++. All the following code was ran on Godbolt.
If I run the following code
std::complex<float> a(-3.0, 4.0);
std::complex<float> inf(1.0f / 0.0f, 0.0f);
a /= inf;
std::cout << a << std::endl;
Then the output will be (-0,0), which is expected. To implement this in HLSL I copied the definition for operator/= found here (same gcc version as Godbolt to make sure).
Now running this code:
std::complex<float> a(-3.0, 4.0);
std::complex<float> inf(1.0f / 0.0f, 0.0f);
float r = a.real() * inf.real() + a.imag() * inf.imag();
float n = std::norm(inf);
float imag = (a.imag() * inf.real() - a.real() * inf.imag()) / n;
float real = r / n;
std::cout << std::complex<float>(real, imag) << std::endl;
The output is now (-nan, -nan) which is different than what I got from just doing /=, but is what I expect from reading the code
Is the standard library doing something I'm not seeing in the source code, like some sort of sanity check? What am I missing?
You are looking at the wrong definition of operator/=. You are looking at the definition for the primary class template template<typename T> complex<T>. There is a specialization template<> complex<float> with its own operator/=. This function uses the built in __complex__ float type of GCC (__complex__ here is a keyword modifying float, like unsigned in unsigned int). Essentially, GCC implements C++'s std::complex in terms of its existing support for C's _Complex/complex types (which don't exist in standard C++).
std::complex<float> a(-3.0, 4.0);
std::complex<float> inf(1.0f / 0.0f, 0.0f);
// copy a and inf into the compiler-supported complex types
__complex__ float a_prim, inf_prim;
__real__ a_prim = a.real(); __imag__ a_prim = a.imag();
__real__ inf_prim = inf.real(); __imag__ inf_prim = inf.imag();
// use primitive complex division
a_prim /= inf_prim;
std::cout << std::complex<float>(__real__ a_prim, __imag__ a_prim) << "\n";
A comment in the GCC source suggests that __complex__ division eventually does compile down to a call to library function. (Also borne out by looking at the generated code for the above snippet.) That library is libgcc, the GCC runtime library. Indeed, in libgcc, there's documentation for a function that computes a complex quotient. You can stare at the (very long) definition here. (Note that there's some macro magic going on to define all the variants of the function for all the floating point types at once).