mpcomplex.hpp

Go to the documentation of this file.

#pragma once
 
#include "mpreal.hpp"
#if !UNI20_HAS_MPREAL
#error "uni20/core/mpcomplex.hpp requires UNI20_ENABLE_MPREAL=ON (GMP, MPFR and MPC)"
#endif
#include <mpc.h>
 
namespace uni20
{
namespace detail
{
struct mpcomplex_access;
struct exact_complex
{
    exact_rational re;
    exact_rational im;
};
struct mpc_value
{
    mpc_t value{};
    bool initialized = false;
    explicit mpc_value(Precision p)
    {
      require_mpfr_tls();
      mpc_init2(value, mpfr_precision(p));
      initialized = true;
      mpc_set_ui(value, 0, MPC_RNDNN);
    }
    mpc_value(mpc_value const& other) : mpc_value(Precision::bits(mpc_get_prec(other.value)))
    {
      mpc_set(value, other.value, MPC_RNDNN);
    }
    mpc_value(mpc_value&& other) noexcept
    {
      std::swap(value[0], other.value[0]);
      std::swap(initialized, other.initialized);
    }
    mpc_value& operator=(mpc_value other) noexcept
    {
      std::swap(value[0], other.value[0]);
      std::swap(initialized, other.initialized);
      return *this;
    }
    ~mpc_value()
    {
      if (initialized) mpc_clear(value);
    }
};
} // namespace detail
 
class mpcomplex {
  public:
    using value_type = mpreal;
    mpcomplex() : value_(std::in_place_type<detail::exact_complex>) {}
    explicit mpcomplex(uninitialized_t) noexcept : value_(std::monostate{}) {}
    explicit mpcomplex(Precision p) : mpcomplex(uninitialized)
    {
      if (p.is_exact())
        value_.emplace<detail::exact_complex>();
      else
        value_.emplace<detail::mpc_value>(p);
    }
    mpcomplex(mpreal const& re, mpreal const& im, Precision p) : mpcomplex(p)
    {
      if (p.is_exact())
      {
        if (!re.precision().is_exact() || !im.precision().is_exact())
          throw std::invalid_argument("complex<mpreal>: an approximation cannot be reclassified as exact");
        value_ = detail::exact_complex{re.exact_value(), im.exact_value()};
      }
      else
      {
        re.copy_to(mpc_realref(this->approximate()));
        im.copy_to(mpc_imagref(this->approximate()));
      }
    }
    mpcomplex(mpreal const& re, mpreal const& im) : mpcomplex(re, im, common_precision(re.precision(), im.precision()))
    {}
    mpcomplex(mpreal const& re) : mpcomplex(re, mpreal{}, re.precision()) {}
    mpcomplex(mpreal const& re, Precision p) : mpcomplex(re, mpreal{}, p) {}
    template <std::integral I>
      requires(!std::same_as<I, bool>)
    explicit mpcomplex(I value) : mpcomplex(mpreal(value))
    {}
    explicit mpcomplex(exact_rational const& value) : mpcomplex(mpreal(value)) {}
    explicit mpcomplex(decimal_literal value) : mpcomplex(mpreal(value)) {}
    template <std::integral I>
      requires(!std::same_as<I, bool>)
    mpcomplex(I value, Precision p) : mpcomplex(mpreal(value, p))
    {}
    mpcomplex(exact_rational const& value, Precision p) : mpcomplex(mpreal(value, p)) {}
    mpcomplex(decimal_literal value, Precision p) : mpcomplex(mpreal(value, p)) {}
    mpcomplex(std::string_view re, std::string_view im, Precision p) : mpcomplex(mpreal(re, p), mpreal(im, p)) {}
    mpcomplex(mpcomplex const& other, Precision p) : mpcomplex(p)
    {
      other.require_value();
      if (p.is_exact())
      {
        if (!other.is_exact())
          throw std::invalid_argument("complex<mpreal>: an approximation cannot be reclassified as exact");
        value_ = std::get<detail::exact_complex>(other.value_);
      }
      else
        other.copy_to(this->approximate());
    }
    explicit mpcomplex(mpc_srcptr value, Precision p) : mpcomplex(p) { mpc_set(this->approximate(), value, MPC_RNDNN); }
    mpcomplex(mpcomplex const&) = default;
    mpcomplex(mpcomplex&& other) noexcept : value_(std::move(other.value_)) { other.value_.emplace<std::monostate>(); }
    mpcomplex& operator=(mpcomplex other) noexcept
    {
      this->swap(other);
      return *this;
    }
    void swap(mpcomplex& other) noexcept { value_.swap(other.value_); }
    bool initialized() const noexcept { return !std::holds_alternative<std::monostate>(value_); }
    bool is_exact() const noexcept { return std::holds_alternative<detail::exact_complex>(value_); }
    Precision precision() const
    {
      this->require_value();
      return this->is_exact() ? Precision::exact() : Precision::bits(mpc_get_prec(this->native_handle()));
    }
    mpc_srcptr native_handle() const
    {
      this->require_value();
      if (this->is_exact())
        throw std::logic_error("complex<mpreal>: exact value needs finite precision before MPC access");
      return std::get<detail::mpc_value>(value_).value;
    }
    void copy_to(mpc_ptr destination) const
    {
      this->require_value();
      if (this->is_exact())
      {
        auto const& q = std::get<detail::exact_complex>(value_);
        mpfr_set_q(mpc_realref(destination), q.re.native_handle(), MPFR_RNDN);
        mpfr_set_q(mpc_imagref(destination), q.im.native_handle(), MPFR_RNDN);
      }
      else
        mpc_set(destination, this->native_handle(), MPC_RNDNN);
    }
    mpcomplex at(Precision p) const { return mpcomplex(*this, p); }
    mpreal real() const
    {
      if (this->is_exact()) return mpreal(std::get<detail::exact_complex>(value_).re);
      return mpreal(mpc_realref(this->native_handle()), this->precision());
    }
    mpreal imag() const
    {
      if (this->is_exact()) return mpreal(std::get<detail::exact_complex>(value_).im);
      return mpreal(mpc_imagref(this->native_handle()), this->precision());
    }
    void real(mpreal const& re)
    {
      auto p = this->precision();
      if (p.is_exact())
        *this = mpcomplex(re, this->imag());
      else
        re.copy_to(mpc_realref(this->approximate()));
    }
    void imag(mpreal const& im)
    {
      auto p = this->precision();
      if (p.is_exact())
        *this = mpcomplex(this->real(), im);
      else
        im.copy_to(mpc_imagref(this->approximate()));
    }
    mpcomplex& operator+=(mpcomplex const& rhs);
    mpcomplex& operator-=(mpcomplex const& rhs);
    mpcomplex& operator*=(mpcomplex const& rhs);
    mpcomplex& operator/=(mpcomplex const& rhs);
    mpcomplex& operator+=(exact_rational const& rhs);
    mpcomplex& operator-=(exact_rational const& rhs);
    mpcomplex& operator*=(exact_rational const& rhs);
    mpcomplex& operator/=(exact_rational const& rhs);
 
  private:
    friend struct detail::mpcomplex_access;
    explicit mpcomplex(detail::exact_complex q) : value_(std::move(q)) {}
    void require_value() const
    {
      if (!this->initialized()) throw std::logic_error("complex<mpreal>: numerical use of an unset value");
    }
    mpc_ptr approximate()
    {
      (void)this->native_handle();
      return std::get<detail::mpc_value>(value_).value;
    }
    std::variant<std::monostate, detail::exact_complex, detail::mpc_value> value_;
};
 
namespace detail
{
struct mpcomplex_access
{
    template <auto Operation> static mpcomplex unary(mpcomplex const& x, Precision p)
    {
      (void)p.bit_count();
      return unary<Operation>(x.at(p));
    }
    template <auto Operation> static mpcomplex binary(mpcomplex const& x, mpcomplex const& y, Precision p)
    {
      (void)p.bit_count();
      return binary<Operation>(x.at(p), y.at(p));
    }
    template <auto Operation> static mpreal real_result(mpcomplex const& x, Precision p)
    {
      (void)p.bit_count();
      return real_result<Operation>(x.at(p));
    }
 
    static mpcomplex exact_sqrt(mpcomplex const& x)
    {
      auto const& q = std::get<exact_complex>(x.value_);
      auto magnitude = uni20::sqrt(q.re * q.re + q.im * q.im);
      auto re = uni20::sqrt((magnitude + q.re) / 2);
      auto im = uni20::sqrt((magnitude - q.re) / 2);
      // Exact zero has no sign: the principal root uses the upper cut side.
      if (q.im < 0) im = -im;
      return mpcomplex(exact_complex{std::move(re), std::move(im)});
    }
    template <auto Operation> static mpcomplex unary(mpcomplex const& x)
    {
      auto p = x.precision();
      if (p.is_exact())
      {
        auto const& q = std::get<exact_complex>(x.value_);
        if constexpr (same_function_v<Operation, mpc_neg>)
          return mpcomplex(exact_complex{-q.re, -q.im});
        else if constexpr (same_function_v<Operation, mpc_conj>)
          return mpcomplex(exact_complex{q.re, -q.im});
        else if constexpr (same_function_v<Operation, mpc_sqrt>)
          return exact_sqrt(x);
        else
        {
          if constexpr (same_function_v<Operation, mpc_exp> || same_function_v<Operation, mpc_cos> ||
                        same_function_v<Operation, mpc_cosh>)
          {
            if (q.re == 0 && q.im == 0) return mpcomplex{1};
          }
          else if constexpr (same_function_v<Operation, mpc_log> || same_function_v<Operation, mpc_acosh>)
          {
            if (q.re == 1 && q.im == 0) return mpcomplex{};
          }
          else if constexpr (same_function_v<Operation, mpc_sin> || same_function_v<Operation, mpc_sinh>)
          {
            if (q.re == 0 && q.im == 0) return mpcomplex{};
          }
          throw std::logic_error("complex<mpreal>: elementary function requires finite working precision");
        }
      }
      mpcomplex out(p);
      Operation(out.approximate(), x.native_handle(), MPC_RNDNN);
      return out;
    }
    template <auto Operation> static mpcomplex mixed_finite(mpcomplex const& a, mpcomplex const& b, Precision p)
    {
      auto ar = a.real(), ai = a.imag(), br = b.real(), bi = b.imag();
      if constexpr (same_function_v<Operation, mpc_add>)
        return mpcomplex(ar + br, ai + bi, p);
      else if constexpr (same_function_v<Operation, mpc_sub>)
        return mpcomplex(ar - br, ai - bi, p);
      else
      {
        if (!mpreal_access::can_expand_rational(ar) || !mpreal_access::can_expand_rational(ai) ||
            !mpreal_access::can_expand_rational(br) || !mpreal_access::can_expand_rational(bi))
        {
          auto xr = mpreal_access::stored_scaled(ar), xi = mpreal_access::stored_scaled(ai);
          auto yr = mpreal_access::stored_scaled(br), yi = mpreal_access::stored_scaled(bi);
          auto rr = xr * yr, ii = xi * yi, ri = xr * yi, ir = xi * yr;
          mpcomplex result(p);
          if constexpr (same_function_v<Operation, mpc_mul>)
          {
            scaled_rational_ratio({rr, -ii, {exact_rational{1}, 0}, {}}).round_to(mpc_realref(result.approximate()));
            scaled_rational_ratio({ri, ir, {exact_rational{1}, 0}, {}}).round_to(mpc_imagref(result.approximate()));
          }
          else
          {
            auto dr = yr * yr, di = yi * yi;
            scaled_rational_ratio({rr, ii, dr, di}).round_to(mpc_realref(result.approximate()));
            scaled_rational_ratio({ir, -ri, dr, di}).round_to(mpc_imagref(result.approximate()));
          }
          // Cancellation produces +0. Two zero products preserve the same sign
          // rules as the ordinary rational path; underflow signs stay intact.
          bool const rr_negative = signbit(ar) != signbit(br), ii_negative = signbit(ai) != signbit(bi);
          bool const ri_negative = signbit(ar) != signbit(bi), ir_negative = signbit(ai) != signbit(br);
          if (rr.coefficient == 0 && ii.coefficient == 0)
            mpfr_set_zero(mpc_realref(result.approximate()),
                          (rr_negative && (same_function_v<Operation, mpc_mul> ? !ii_negative : ii_negative)) ? -1 : 1);
          if (ri.coefficient == 0 && ir.coefficient == 0)
            mpfr_set_zero(mpc_imagref(result.approximate()),
                          (ir_negative && (same_function_v<Operation, mpc_mul> ? ri_negative : !ri_negative)) ? -1 : 1);
          return result;
        }
        auto xr = mpreal_access::stored_rational(ar), xi = mpreal_access::stored_rational(ai);
        auto yr = mpreal_access::stored_rational(br), yi = mpreal_access::stored_rational(bi);
        // Form each component exactly before its single final rounding. Keep
        // signed-zero rules from the component products: rational zero alone
        // does not distinguish the two sides of a complex branch cut.
        auto rr = xr * yr, ii = xi * yi, ri = xr * yi, ir = xi * yr;
        bool const rr_negative = signbit(ar) != signbit(br);
        bool const ii_negative = signbit(ai) != signbit(bi);
        bool const ri_negative = signbit(ar) != signbit(bi);
        bool const ir_negative = signbit(ai) != signbit(br);
        exact_complex value;
        bool re_negative_zero, im_negative_zero;
        if constexpr (same_function_v<Operation, mpc_mul>)
        {
          value = {rr - ii, ri + ir};
          re_negative_zero = rr == 0 && ii == 0 && rr_negative && !ii_negative;
          im_negative_zero = ri == 0 && ir == 0 && ri_negative && ir_negative;
        }
        else
        {
          auto d = yr * yr + yi * yi;
          value = {(rr + ii) / d, (ir - ri) / d};
          re_negative_zero = rr == 0 && ii == 0 && rr_negative && ii_negative;
          im_negative_zero = ir == 0 && ri == 0 && ir_negative && !ri_negative;
        }
        mpcomplex result(value);
        result = result.at(p);
        if (value.re == 0) mpfr_set_zero(mpc_realref(result.approximate()), re_negative_zero ? -1 : 1);
        if (value.im == 0) mpfr_set_zero(mpc_imagref(result.approximate()), im_negative_zero ? -1 : 1);
        return result;
      }
    }
 
    template <auto Operation> static mpcomplex binary(mpcomplex const& a, mpcomplex const& b)
    {
      auto p = common_precision(a.precision(), b.precision());
      if (p.is_exact())
      {
        auto const& x = std::get<exact_complex>(a.value_);
        auto const& y = std::get<exact_complex>(b.value_);
        if constexpr (same_function_v<Operation, mpc_add>)
          return mpcomplex(exact_complex{x.re + y.re, x.im + y.im});
        else if constexpr (same_function_v<Operation, mpc_sub>)
          return mpcomplex(exact_complex{x.re - y.re, x.im - y.im});
        else if constexpr (same_function_v<Operation, mpc_mul>)
          return mpcomplex(exact_complex{x.re * y.re - x.im * y.im, x.re * y.im + x.im * y.re});
        else if constexpr (same_function_v<Operation, mpc_div>)
        {
          auto d = y.re * y.re + y.im * y.im;
          return mpcomplex(exact_complex{(x.re * y.re + x.im * y.im) / d, (x.im * y.re - x.re * y.im) / d});
        }
        else
        {
          if constexpr (same_function_v<Operation, mpc_pow>)
          {
            if (y.im == 0)
            {
              auto exponent = y.re.native_handle();
              if (mpz_cmp_ui(mpq_denref(exponent), 1) == 0) return integer_power(a, mpq_numref(exponent));
              if (mpz_cmp_ui(mpq_denref(exponent), 2) == 0) return integer_power(exact_sqrt(a), mpq_numref(exponent));
              if (x.im == 0 && x.re >= 0) return mpcomplex(uni20::pow(x.re, y.re));
            }
            if (x.re == 1 && x.im == 0) return mpcomplex{1};
          }
          throw std::logic_error("complex<mpreal>: elementary function requires finite working precision");
        }
      }
      if constexpr (same_function_v<Operation, mpc_add> || same_function_v<Operation, mpc_sub> ||
                    same_function_v<Operation, mpc_mul> || same_function_v<Operation, mpc_div>)
      {
        if (a.is_exact() != b.is_exact())
        {
          auto const& approximate = a.is_exact() ? b : a;
          auto native = approximate.native_handle();
          bool const finite = mpfr_number_p(mpc_realref(native)) && mpfr_number_p(mpc_imagref(native));
          if (finite && (!same_function_v<Operation, mpc_div> || b.real() != 0 || b.imag() != 0))
            return mixed_finite<Operation>(a, b, p);
        }
      }
      mpcomplex av(uninitialized), bv(uninitialized);
      if (a.is_exact()) av = a.at(p);
      if (b.is_exact()) bv = b.at(p);
      mpcomplex out(p);
      Operation(out.approximate(), a.is_exact() ? av.native_handle() : a.native_handle(),
                b.is_exact() ? bv.native_handle() : b.native_handle(), MPC_RNDNN);
      return out;
    }
    template <auto Operation> static mpcomplex& assign(mpcomplex& a, mpcomplex const& b)
    {
      if (!a.is_exact() && !b.is_exact())
      {
        (void)common_precision(a.precision(), b.precision());
        Operation(a.approximate(), a.native_handle(), b.native_handle(), MPC_RNDNN);
      }
      else
        a = binary<Operation>(a, b);
      return a;
    }
    template <auto Operation> static mpreal real_result(mpcomplex const& x)
    {
      auto p = x.precision();
      if (p.is_exact())
      {
        if constexpr (same_function_v<Operation, mpc_norm>)
        {
          auto const& q = std::get<exact_complex>(x.value_);
          return mpreal(q.re * q.re + q.im * q.im);
        }
        else
        {
          auto const& q = std::get<exact_complex>(x.value_);
          if constexpr (same_function_v<Operation, mpc_abs>)
            return mpreal(uni20::sqrt(q.re * q.re + q.im * q.im));
          else if constexpr (same_function_v<Operation, mpc_arg>)
            if (q.im == 0 && q.re > 0) return mpreal{};
          throw std::logic_error("complex<mpreal>: elementary function requires finite working precision");
        }
      }
      mpreal out(p);
      Operation(out.approximate(), x.native_handle(), MPFR_RNDN);
      return out;
    }
};
} // namespace detail
 
inline mpcomplex operator+(mpcomplex const& a, mpcomplex const& b)
{
  return detail::mpcomplex_access::binary<mpc_add>(a, b);
}
inline mpcomplex operator-(mpcomplex const& a, mpcomplex const& b)
{
  return detail::mpcomplex_access::binary<mpc_sub>(a, b);
}
inline mpcomplex operator*(mpcomplex const& a, mpcomplex const& b)
{
  return detail::mpcomplex_access::binary<mpc_mul>(a, b);
}
inline mpcomplex operator/(mpcomplex const& a, mpcomplex const& b)
{
  return detail::mpcomplex_access::binary<mpc_div>(a, b);
}
inline mpcomplex operator-(mpcomplex const& x) { return detail::mpcomplex_access::unary<mpc_neg>(x); }
inline mpcomplex operator+(mpcomplex const& x)
{
  (void)x.precision();
  return x;
}
inline mpcomplex& mpcomplex::operator+=(mpcomplex const& rhs)
{
  return detail::mpcomplex_access::assign<mpc_add>(*this, rhs);
}
inline mpcomplex& mpcomplex::operator-=(mpcomplex const& rhs)
{
  return detail::mpcomplex_access::assign<mpc_sub>(*this, rhs);
}
inline mpcomplex& mpcomplex::operator*=(mpcomplex const& rhs)
{
  return detail::mpcomplex_access::assign<mpc_mul>(*this, rhs);
}
inline mpcomplex& mpcomplex::operator/=(mpcomplex const& rhs)
{
  return detail::mpcomplex_access::assign<mpc_div>(*this, rhs);
}
 
inline mpcomplex operator+(mpcomplex const& a, exact_rational const& b) { return a + mpcomplex(b); }
inline mpcomplex operator-(mpcomplex const& a, exact_rational const& b) { return a - mpcomplex(b); }
inline mpcomplex operator*(mpcomplex const& a, exact_rational const& b) { return a * mpcomplex(b); }
inline mpcomplex operator/(mpcomplex const& a, exact_rational const& b) { return a / mpcomplex(b); }
inline mpcomplex operator+(exact_rational const& a, mpcomplex const& b) { return mpcomplex(a) + b; }
inline mpcomplex operator-(exact_rational const& a, mpcomplex const& b) { return mpcomplex(a) - b; }
inline mpcomplex operator*(exact_rational const& a, mpcomplex const& b) { return mpcomplex(a) * b; }
inline mpcomplex operator/(exact_rational const& a, mpcomplex const& b) { return mpcomplex(a) / b; }
inline mpcomplex& mpcomplex::operator+=(exact_rational const& rhs) { return *this += mpcomplex(rhs); }
inline mpcomplex& mpcomplex::operator-=(exact_rational const& rhs) { return *this -= mpcomplex(rhs); }
inline mpcomplex& mpcomplex::operator*=(exact_rational const& rhs) { return *this *= mpcomplex(rhs); }
inline mpcomplex& mpcomplex::operator/=(exact_rational const& rhs) { return *this /= mpcomplex(rhs); }
 
inline bool operator==(mpcomplex const& a, mpcomplex const& b) { return a.real() == b.real() && a.imag() == b.imag(); }
inline bool operator==(mpcomplex const& a, exact_rational const& b) { return a.real() == b && a.imag() == 0; }
inline mpcomplex conj(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_conj>(z); }
inline mpreal abs(mpcomplex const& z) { return detail::mpcomplex_access::real_result<mpc_abs>(z); }
inline mpreal norm(mpcomplex const& z) { return detail::mpcomplex_access::real_result<mpc_norm>(z); }
inline mpreal arg(mpcomplex const& z) { return detail::mpcomplex_access::real_result<mpc_arg>(z); }
inline mpreal real(mpcomplex const& z) { return z.real(); }
inline mpreal imag(mpcomplex const& z) { return z.imag(); }
inline bool isfinite(mpcomplex const& z) { return isfinite(z.real()) && isfinite(z.imag()); }
inline bool isnan(mpcomplex const& z) { return isnan(z.real()) || isnan(z.imag()); }
inline bool isinf(mpcomplex const& z) { return isinf(z.real()) || isinf(z.imag()); }
inline mpcomplex sqrt(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_sqrt>(z); }
inline mpcomplex exp(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_exp>(z); }
inline mpcomplex log(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_log>(z); }
inline mpcomplex sin(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_sin>(z); }
inline mpcomplex cos(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_cos>(z); }
inline mpcomplex sinh(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_sinh>(z); }
inline mpcomplex cosh(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_cosh>(z); }
inline mpcomplex acosh(mpcomplex const& z) { return detail::mpcomplex_access::unary<mpc_acosh>(z); }
inline mpcomplex pow(mpcomplex const& a, mpcomplex const& b) { return detail::mpcomplex_access::binary<mpc_pow>(a, b); }
inline mpcomplex sqrt(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::unary<mpc_sqrt>(z, p); }
inline mpcomplex exp(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::unary<mpc_exp>(z, p); }
inline mpcomplex log(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::unary<mpc_log>(z, p); }
inline mpcomplex sin(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::unary<mpc_sin>(z, p); }
inline mpcomplex cos(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::unary<mpc_cos>(z, p); }
inline mpcomplex sinh(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::unary<mpc_sinh>(z, p); }
inline mpcomplex cosh(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::unary<mpc_cosh>(z, p); }
inline mpcomplex acosh(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::unary<mpc_acosh>(z, p); }
inline mpreal abs(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::real_result<mpc_abs>(z, p); }
inline mpreal norm(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::real_result<mpc_norm>(z, p); }
inline mpreal arg(mpcomplex const& z, Precision p) { return detail::mpcomplex_access::real_result<mpc_arg>(z, p); }
inline mpcomplex pow(mpcomplex const& a, mpcomplex const& b, Precision p)
{
  return detail::mpcomplex_access::binary<mpc_pow>(a, b, p);
}
 
inline std::ostream& operator<<(std::ostream& out, mpcomplex const& z)
{
  return out << '(' << z.real() << ',' << z.imag() << ')';
}
} // namespace uni20