Description
Description
Currently, Stan's complex number support is entirely built on std::complex
, including autodiff, which
uses std::complex<stan::math::var>
.
This is, unfortunately, unspecified behavior in the C++ spec [26.4.2]:
The effect of instantiating the template
complex
for any type other thanfloat
,double
, orlong double
is unspecified. The specializationscomplex<float>
,complex<double>
, andcomplex<long double>
are literal types
For a reminder on what "unspecified behavior" means:
unspecified behavior - the behavior of the program varies between implementations, and the conforming implementation is not required to document the effects of each behavior. Each unspecified behavior results in one of a set of valid results.
Essentially, unspecified behavior is the same as "implementation-defined behavior" but without the requirement that implementations document what they are doing. This is also often taken to mean there are no backwards compatibility guarantees on any specific unspecified behavior.
This creates both a maintenance burden (each new libstdc++/libc++ release can create arbitrary amounts of work for our developers) and a stability hazard (the idea that "Stan X.Y will continue to work a year from now, without needing to update to Stan X.Z" is false as things stand today)
Problems
Recent versions of clang
/libstdc++
have made changes which they are fully within their rights to do by the spec, but have broken Stan builds.
- In libstdc++16, they changed the definition of
log(complex)
fromcomplex<T>(log(abs(x)), arg(x));
tocomplex<T>(std::log(std::abs(x)), std::arg(x));
. This broke argument dependent lookup for this function.
A similar change brokeoperator*
for our complex types.
This lead to Unable to completemake build
on Mac M1 cmdstan#1158, which was the reason we needed a 2.32.1 release.
@andrjohns provided the fix in Add overloads for complex multiply to fix clang16 #2892 - In libstdc++17, a similar change was made to
fabs
, which necessitated to Fix usages of fabs in check_symmetric with llvm17 #2991 - In libstdc++19, the internal structure of
pow
was rewritten such that several overloads lead to a static assert failing if the type passed was not arithmetic: Compilation failures under LLVM 19 #3106
What to do
This is less clear to me.
Option 1 - walk on egg shells
So far, all of the issues that have arisen from this have been due to argument dependent lookup breaking for these types. We can fix that by being much more explicit, as we did in #2892 and #2991. This requires auditing the existing usages, which probably requires a fair amount of C++ expertise to understand how the calls are being resolved.
Option 2 - our own type
We could rather trivially define our own stan::math::complex<T>
type. We could make it assignable from std::complex<double>
, and I think be off to the races? I believe the complex linear algebra we use in Eigen all support a template argument for the complex type, rather than assuming std::complex
.
This would require a fair amount of boilerplate to actually do any math on it, and in the case of double
we may lose out on some of the optimizations that having the type built in to the language grants, but we'd own it.