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]

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).

❗ 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

@@            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	Coverage Δ
test/test_complex_class.cpp	100.0% <100.0%> (ø)
include/boost/multiprecision/complex.hpp	90.3% <90.3%> (ø)

... and 3 files with indirect coverage changes

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 7646f8e...903d23c. Read the comment docs.

[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: f61a9e1. 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: 612f4fe

[jzmaddock]

@mborland you're re-inventing the wheel here? All these complex number functions are already implemented, it's just that number<whatever> is already a complex number type with the same interface as std::complex for a suitably defined whatever. See https://www.boost.org/doc/libs/1_85_0/libs/multiprecision/doc/html/boost_multiprecision/tut/complex.html

[mborland]

In @NAThompson's original issue he was using cpp_bin_float_100 as the number type. The templates could be enable_if for the already complex types, but if we want to allow a wrapper of any floating-point multiprecision type I don't see any other way around this.

[jzmaddock]

I guess my point is that adding a boost::multiprecision::complex class duplicates what is already there and adds confusion (and maintenance), what if you had (caution untested code ahead!!):

template <typename T, boost::multiprecision::expression_template_option ET>
class complex<boost::multiprecision::number<T, ET>>
{
  using self_type = complex<boost::multiprecision::number<T, ET>>;
  using complex_type = decltype(polar(std::declval<boost::multiprecision::number<T, ET>>(), std::declval<boost::multiprecision::number<T, ET>>()));
  complex_type m_data;

  complex_type& data(){ return m_data; }
  const complex_type& data()const{ return m_data; }
public:
  // public members here..

  // and then:
  friend inline self_type operator+(const self_type& a, const self_type& b)
  { return a.data() + b.data(); }
  // etc
  friend inline self_type exp(const self_type& val) 
  { return exp(val.data()); }
  // etc
};

Now there's no duplication, so bugs fixed in one place are fixed everywhere. Plus if the user has included mpc.hpp for example then mpc will be used for mpfr complex number types rather than synthesising our own type.

It's still a bind writing out all the forwarding functions, but at least they're trivial and improvements in one place propagate throughout.

[mborland]

I can give that a go tomorrow. Forwarding functions will be monotonous which is fine. What is non-obvious to me is how to overload say self_type + complex type like is here: https://github.com/boostorg/multiprecision/pull/620/files#diff-db4d91a1fc843b60d38c832c49b4b5ec9b8c22f772d5c5d0a56eea563ae59a5dR50. I'll play with it

[jzmaddock]

I think you have the form of the function already there?

BTW those functions are definitely NOT noexcept in the multiprecision case: memory allocations usually are involved!

[mborland]

@jzmaddock is 9b3575b 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 9b3575b 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