Complex Overloads

[mborland]

@NAThompson I started poking around on this after you asked. I found the only real way to have viable overloads for multiprecision types is to inject them into namespace std. The STLs use unconstrained templates for their functions unlike <cmath> From libc++:

template<class _Tp>
inline _LIBCPP_INLINE_VISIBILITY
_Tp
abs(const complex<_Tp>& __c)
{
    return std::hypot(__c.real(), __c.imag());
}

I also tried using the derived class (the giant commented out class) like:

namespace std {

template <>
class complex<boost::multiprecision::cpp_bin_float_50> : public boost::multiprecision::complex<boost::multiprecision::cpp_bin_float_50> {};

I don't think this is explicit illegal since the standard only specifies T - the type of the real and imaginary parts. The behavior is unspecified (and may fail to compile) if T is not a cv-unqualified standard(until C++23) floating-point type and undefined if T is not NumericType

CC: @jzmaddock

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 98.18653% with 7 lines in your changes are missing coverage. Please review.

Project coverage is 94.2%. Comparing base (7646f8e) to head (ded45c0).

:exclamation: Current head ded45c0 differs from pull request most recent head 903d23c

Please upload reports for the commit 903d23c to get more accurate results.

Additional details and impacted files

[](https://app.codecov.io/gh/boostorg/multiprecision/pull/620?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg) ```diff @@ Coverage Diff @@ ## develop #620 +/- ## ========================================= + Coverage 94.1% 94.2% +0.1% ========================================= Files 274 276 +2 Lines 26766 27152 +386 ========================================= + Hits 25163 25553 +390 + Misses 1603 1599 -4 ``` | [Files](https://app.codecov.io/gh/boostorg/multiprecision/pull/620?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg) | Coverage Δ | | |---|---|---| | [test/test\_complex\_class.cpp](https://app.codecov.io/gh/boostorg/multiprecision/pull/620?src=pr&el=tree&filepath=test%2Ftest_complex_class.cpp&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg#diff-dGVzdC90ZXN0X2NvbXBsZXhfY2xhc3MuY3Bw) | `100.0% <100.0%> (ø)` | | | [include/boost/multiprecision/complex.hpp](https://app.codecov.io/gh/boostorg/multiprecision/pull/620?src=pr&el=tree&filepath=include%2Fboost%2Fmultiprecision%2Fcomplex.hpp&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg#diff-aW5jbHVkZS9ib29zdC9tdWx0aXByZWNpc2lvbi9jb21wbGV4LmhwcA==) | `90.3% <90.3%> (ø)` | | ... and [3 files with indirect coverage changes](https://app.codecov.io/gh/boostorg/multiprecision/pull/620/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg) ------ [Continue to review full report in Codecov by Sentry](https://app.codecov.io/gh/boostorg/multiprecision/pull/620?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg). > **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg) > `Δ = absolute (impact)`, `ø = not affected`, `? = missing data` > Powered by [Codecov](https://app.codecov.io/gh/boostorg/multiprecision/pull/620?dropdown=coverage&src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg). Last update [7646f8e...903d23c](https://app.codecov.io/gh/boostorg/multiprecision/pull/620?dropdown=coverage&src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=boostorg).

[jzmaddock]

I suspect we get into trouble whatever here: about the only thing we're legally allowed to do is [partially-] specialise std::complex for type number. The non member operators I think we can put in an "associated" namespace ie boost::multiprecision and ADL will still find them. What we're not allowed to do, is add overloads for say std::real in namespace std. I have a hunch potential bear traps lurk if we do that, but I don't see any good alternative.

Implementation wise, since it's all already written (albeit under another name), I would just make std::complex<boost::multiprecision::number<T, ET>> a thin wrapper around our existing types (multiprecision does have an internal traits class for that somewhere - in order to implement polar). That would then preserve things like the mpc wrapper where the backend can do complex arithmetic better than we can.

[NAThompson]

Is there sense in trying to upstream some of this to clang? It appears that they simply haven't imagined this usecase.

[mborland]

Is there sense in trying to upstream some of this to clang? It appears that they simply haven't imagined this usecase.

I don't think so. Up to C++23 the standard specified standard floating point types, and since then it says trivially copyable literal numeric types. I assume they changed the wording to fit the types of <stdfloat>. I don't believe most of the multiprecision types (except maybe float128) meet that requirement.

[jzmaddock]

It's a long standing C++ library issue: and tightened up in C++23, so for example MSVC will now static_assert if you try and use std::complex with anything other than the mandated types.

[mborland]

I suspect we get into trouble whatever here: about the only thing we're legally allowed to do is [partially-] specialise std::complex for type number. The non member operators I think we can put in an "associated" namespace ie boost::multiprecision and ADL will still find them. What we're not allowed to do, is add overloads for say std::real in namespace std. I have a hunch potential bear traps lurk if we do that, but I don't see any good alternative.

Implementation wise, since it's all already written (albeit under another name), I would just make std::complex<boost::multiprecision::number<T, ET>> a thin wrapper around our existing types (multiprecision does have an internal traits class for that somewhere - in order to implement polar). That would then preserve things like the mpc wrapper where the backend can do complex arithmetic better than we can.

This works locally. Added in: https://github.com/boostorg/multiprecision/pull/620/commits/f61a9e1a976d21da39e689144e3491a818b615a1. How would you use the existing backend implementation of polar instead of the one provided? Simple enable_if in the event we know the type has that functionallity?

edit: Figured it out here: https://github.com/boostorg/multiprecision/pull/620/commits/612f4fe2dfbfd314bf734d212320631b83786b48

[mborland]

@jzmaddock is https://github.com/boostorg/multiprecision/pull/620/commits/9b3575b34522281283d6df3499e458fb02ae210e more what you were expecting? The only thing I don't see a way around is overloading polar to return std::complex. The existing overload in default ops takes two number types so the only difference between the two would be the return type which is not allowed.

[jzmaddock]

Apologies Matt, have been travelling, will try and look at this tomorrow...

[jzmaddock]

@jzmaddock is https://github.com/boostorg/multiprecision/commit/9b3575b34522281283d6df3499e458fb02ae210e more what you were expecting?

Yup, thanks!

The only thing I don't see a way around is overloading polar to return std::complex. The existing overload in default ops takes two number types so the only difference between the two would be the return type which is not allowed.

Oh :( Are they not OK as long as they are in different namespaces?

[mborland]

Oh :( Are they not OK as long as they are in different namespaces?

Interestingly you would think it would be fine. If I move all the new stuff into boost::multiprecision::complex and then add using namespace boost::multiprecision::complex to the test file I start getting ADL errors where multiprecision types fail during substituion into the STL functions. Maybe I am missing something?

[jzmaddock]

Interestingly you would think it would be fine. If I move all the new stuff into boost::multiprecision::complex and then add using namespace boost::multiprecision::complex to the test file I start getting ADL errors where multiprecision types fail during substituion into the STL functions. Maybe I am missing something?

The only solution I could find for now, was to place polar in namespace std, and make all calls to polar explicit: ADL will never work here, because the multiprecision overload will always be found via ADL, and that's the overload we don't want!

Tentatively pushed.

[mborland]

@jzmaddock since this directly overloads std::complex do we need to add a specialization of Eigen::NumTraits beyond what you have here: https://github.com/boostorg/multiprecision/blob/develop/include/boost/multiprecision/eigen.hpp

[jzmaddock]

since this directly overloads std::complex do we need to add a specialization of Eigen::NumTraits beyond what you have here

I don't think so but could be wrong, has this been tried in Nick's original problem?

[NAThompson]

I don't think so but could be wrong, has this been tried in Nick's original problem?

@mborland: The polynomial roots code is pretty trivial; want to just merge my branch into this one and close the other PR?

[mborland]

I don't think so but could be wrong, has this been tried in Nick's original problem?

@mborland: The polynomial roots code is pretty trivial; want to just merge my branch into this one and close the other PR?

Does it belong in multi-precision since right now the PR is against math. I added the test set from your comment in the most recent two commits.

[NAThompson]

@mborland : Yup, you're right . . . that was a silly comment earlier.

[ckormanyos]

Hi Matt (@mborland) I just ran this thing through some local implementations I have of the complex-valued Riemann-Zeta function.

Just a heads-up I think the complex-adaption is missing overloads for $x^n$, where $x$ is multiprecision and real-valued, yet $n$ is either int or possibly a signed or unsigned integral type (where this means integral in the sense of std::is_integral).

MSVC when doing those integral pawers ended up calling std::pow and kicking out the multiple-precision. Whereby GCC's implemenation of <complex> somewhow got it right. Nonetheless, I think actually present clear overloads are needed for complex-to-power-of-integral.

Cc: @jzmaddock

[ckormanyos]

The so-called ladder method comes to mind, which I have appended from my own, self-written extended_complex class years ago.

    template<typename T, typename EnableType> complex<T, EnableType> pow(const complex<T, EnableType>& my_z, int my_pn)
    {
      if     (my_pn <  static_cast<int>(INT8_C(0))) { return  T(static_cast<unsigned>(UINT8_C(1))) / pow(my_z, -my_pn); }
      else if(my_pn == static_cast<int>(INT8_C(0))) { return  complex<T, EnableType>(T(static_cast<unsigned>(UINT8_C(1)))); }
      else if(my_pn == static_cast<int>(INT8_C(1))) { return  my_z; }
      else if(my_pn == static_cast<int>(INT8_C(2))) { return  my_z * my_z; }
      else if(my_pn == static_cast<int>(INT8_C(3))) { return (my_z * my_z) * my_z; }
      else if(my_pn == static_cast<int>(INT8_C(4))) { const complex<T, EnableType> my_z2(my_z * my_z); return (my_z2 * my_z2); }
      else
      {
        complex<T, EnableType> my_result = T(static_cast<unsigned>(UINT8_C(1)));

        complex<T, EnableType> y(my_z);

        auto p_local = static_cast<std::uint64_t>(my_pn);

        // Use the so-called ladder method for the power calculation.
        for(;;)
        {
          const auto lo_bit = static_cast<std::uint_fast8_t>(p_local & static_cast<unsigned>(UINT8_C(1)));

          const auto do_power_multiply = (lo_bit != static_cast<std::uint_fast8_t>(UINT8_C(0)));

          if(do_power_multiply)
          {
            my_result *= y;
          }

          p_local >>= static_cast<unsigned>(UINT8_C(1));

          if(p_local == static_cast<std::uint64_t>(UINT8_C(0)))
          {
            break;
          }

          y *= y;
        }

        return my_result;
      }
    }

Cc: @mborland and @jzmaddock