boost/multiprecision/complex_adaptor.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //  Copyright 2018 John Maddock. Distributed under the Boost
0003 //  Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MP_COMPLEX_ADAPTOR_HPP
0007 #define BOOST_MP_COMPLEX_ADAPTOR_HPP
0008
0009 #include <boost/multiprecision/number.hpp>
0010 #include <cstdint>
0011 #include <boost/multiprecision/detail/digits.hpp>
0012 #include <boost/multiprecision/detail/hash.hpp>
0013 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
0014 #include <cmath>
0015 #include <algorithm>
0016 #include <complex>
0017
0018 namespace boost {
0019 namespace multiprecision {
0020 namespace backends {
0021
0022 template <class Backend>
0023 struct complex_adaptor
0024 {
0025  protected:
0026    Backend m_real, m_imag;
0027
0028  public:
0029    Backend& real_data()
0030    {
0031       return m_real;
0032    }
0033    const Backend& real_data() const
0034    {
0035       return m_real;
0036    }
0037    Backend& imag_data()
0038    {
0039       return m_imag;
0040    }
0041    const Backend& imag_data() const
0042    {
0043       return m_imag;
0044    }
0045
0046    using signed_types = typename Backend::signed_types  ;
0047    using unsigned_types = typename Backend::unsigned_types;
0048    using float_types = typename Backend::float_types   ;
0049    using exponent_type = typename Backend::exponent_type ;
0050
0051    complex_adaptor() {}
0052    complex_adaptor(const complex_adaptor& o) : m_real(o.real_data()), m_imag(o.imag_data()) {}
0053    // Rvalue construct:
0054    complex_adaptor(complex_adaptor&& o) : m_real(std::move(o.real_data())), m_imag(std::move(o.imag_data()))
0055    {}
0056    complex_adaptor(const Backend& val)
0057        : m_real(val)
0058    {}
0059
0060    template <class T>
0061    complex_adaptor(const T& val, const typename std::enable_if<std::is_convertible<T, Backend>::value>::type* = nullptr)
0062        : m_real(val)
0063    {}
0064
0065    complex_adaptor(const std::complex<float>& val)
0066    {
0067       m_real = (long double)val.real();
0068       m_imag = (long double)val.imag();
0069    }
0070    complex_adaptor(const std::complex<double>& val)
0071    {
0072       m_real = (long double)val.real();
0073       m_imag = (long double)val.imag();
0074    }
0075    complex_adaptor(const std::complex<long double>& val)
0076    {
0077       m_real = val.real();
0078       m_imag = val.imag();
0079    }
0080    template <class T, class U>
0081    complex_adaptor(const T& a, const U& b, typename std::enable_if<std::is_constructible<Backend, T const&>::value&& std::is_constructible<Backend, U const&>::value>::type const* = nullptr)
0082       : m_real(a), m_imag(b) {}
0083    template <class T, class U>
0084    complex_adaptor(T&& a, const U& b, typename std::enable_if<std::is_constructible<Backend, T>::value&& std::is_constructible<Backend, U>::value>::type const* = nullptr)
0085       : m_real(static_cast<T&&>(a)), m_imag(b) {}
0086    template <class T, class U>
0087    complex_adaptor(T&& a, U&& b, typename std::enable_if<std::is_constructible<Backend, T>::value&& std::is_constructible<Backend, U>::value>::type const* = nullptr)
0088       : m_real(static_cast<T&&>(a)), m_imag(static_cast<U&&>(b)) {}
0089    template <class T, class U>
0090    complex_adaptor(const T& a, U&& b, typename std::enable_if<std::is_constructible<Backend, T>::value&& std::is_constructible<Backend, U>::value>::type const* = nullptr)
0091       : m_real(a), m_imag(static_cast<U&&>(b)) {}
0092
0093    complex_adaptor& operator=(const complex_adaptor& o)
0094    {
0095       m_real = o.real_data();
0096       m_imag = o.imag_data();
0097       return *this;
0098    }
0099    // rvalue assign:
0100    complex_adaptor& operator=(complex_adaptor&& o) noexcept
0101    {
0102       m_real = std::move(o.real_data());
0103       m_imag = std::move(o.imag_data());
0104       return *this;
0105    }
0106    template <class V>
0107    typename std::enable_if<std::is_assignable<Backend, V>::value, complex_adaptor&>::type operator=(const V& v)
0108    {
0109       using ui_type = typename std::tuple_element<0, unsigned_types>::type;
0110       m_real = v;
0111       m_imag = ui_type(0u);
0112       return *this;
0113    }
0114    template <class T>
0115    complex_adaptor& operator=(const std::complex<T>& val)
0116    {
0117       m_real = (long double)val.real();
0118       m_imag = (long double)val.imag();
0119       return *this;
0120    }
0121    complex_adaptor& operator=(const char* s)
0122    {
0123       using ui_type = typename std::tuple_element<0, unsigned_types>::type;
0124       ui_type                                           zero = 0u;
0125
0126       using default_ops::eval_fpclassify;
0127
0128       if (s && (*s == '('))
0129       {
0130          std::string part;
0131          const char* p = ++s;
0132          while (*p && (*p != ',') && (*p != ')'))
0133             ++p;
0134          part.assign(s, p);
0135          if (part.size())
0136             real_data() = part.c_str();
0137          else
0138             real_data() = zero;
0139          s = p;
0140          if (*p && (*p != ')'))
0141          {
0142             ++p;
0143             while (*p && (*p != ')'))
0144                ++p;
0145             part.assign(s + 1, p);
0146          }
0147          else
0148             part.erase();
0149          if (part.size())
0150             imag_data() = part.c_str();
0151          else
0152             imag_data() = zero;
0153
0154          if (eval_fpclassify(imag_data()) == static_cast<int>(FP_NAN))
0155          {
0156             real_data() = imag_data();
0157          }
0158       }
0159       else
0160       {
0161          real_data() = s;
0162          imag_data() = zero;
0163       }
0164       return *this;
0165    }
0166
0167    int compare(const complex_adaptor& o) const
0168    {
0169       // They are either equal or not:
0170       return (m_real.compare(o.real_data()) == 0) && (m_imag.compare(o.imag_data()) == 0) ? 0 : 1;
0171    }
0172    template <class T>
0173    int compare(const T& val) const
0174    {
0175       using default_ops::eval_is_zero;
0176       return (m_real.compare(val) == 0) && eval_is_zero(m_imag) ? 0 : 1;
0177    }
0178    void swap(complex_adaptor& o)
0179    {
0180       real_data().swap(o.real_data());
0181       imag_data().swap(o.imag_data());
0182    }
0183    std::string str(std::streamsize dig, std::ios_base::fmtflags f) const
0184    {
0185       using default_ops::eval_is_zero;
0186       if (eval_is_zero(imag_data()))
0187          return m_real.str(dig, f);
0188       return "(" + m_real.str(dig, f) + "," + m_imag.str(dig, f) + ")";
0189    }
0190    void negate()
0191    {
0192       m_real.negate();
0193       m_imag.negate();
0194    }
0195
0196    //
0197    // Default precision:
0198    //
0199    static BOOST_MP_CXX14_CONSTEXPR unsigned default_precision() noexcept
0200    {
0201       return Backend::default_precision();
0202    }
0203    static BOOST_MP_CXX14_CONSTEXPR void default_precision(unsigned digits10)
0204    {
0205       Backend::default_precision(digits10);
0206       Backend::thread_default_precision(digits10);
0207    }
0208    static BOOST_MP_CXX14_CONSTEXPR unsigned thread_default_precision() noexcept
0209    {
0210       return Backend::thread_default_precision();
0211    }
0212    static BOOST_MP_CXX14_CONSTEXPR void thread_default_precision(unsigned digits10)
0213    {
0214       Backend::thread_default_precision(digits10);
0215    }
0216    BOOST_MP_CXX14_CONSTEXPR unsigned precision() const noexcept
0217    {
0218       return m_real.precision();
0219    }
0220    BOOST_MP_CXX14_CONSTEXPR void precision(unsigned digits10)
0221    {
0222       m_real.precision(digits10);
0223       m_imag.precision(digits10);
0224    }
0225    //
0226    // Variable precision options:
0227    //
0228    static constexpr variable_precision_options default_variable_precision_options()noexcept
0229    {
0230       return Backend::default_variable_precision_options();
0231    }
0232    static constexpr variable_precision_options thread_default_variable_precision_options()noexcept
0233    {
0234       return Backend::thread_default_variable_precision_options();
0235    }
0236    static BOOST_MP_CXX14_CONSTEXPR void default_variable_precision_options(variable_precision_options opts)
0237    {
0238       Backend::default_variable_precision_options(opts);
0239       Backend::thread_default_variable_precision_options(opts);
0240    }
0241    static BOOST_MP_CXX14_CONSTEXPR void thread_default_variable_precision_options(variable_precision_options opts)
0242    {
0243       Backend::thread_default_variable_precision_options(opts);
0244    }
0245 };
0246
0247 template <class Backend, class T>
0248 inline typename std::enable_if<boost::multiprecision::detail::is_arithmetic<T>::value, bool>::type eval_eq(const complex_adaptor<Backend>& a, const T& b) noexcept
0249 {
0250    return a.compare(b) == 0;
0251 }
0252
0253 template <class Backend>
0254 inline void eval_add(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
0255 {
0256    eval_add(result.real_data(), o.real_data());
0257    eval_add(result.imag_data(), o.imag_data());
0258 }
0259 template <class Backend>
0260 inline void eval_subtract(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
0261 {
0262    eval_subtract(result.real_data(), o.real_data());
0263    eval_subtract(result.imag_data(), o.imag_data());
0264 }
0265 template <class Backend>
0266 inline void eval_multiply(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
0267 {
0268    Backend t1, t2, t3;
0269    eval_multiply(t1, result.real_data(), o.real_data());
0270    eval_multiply(t2, result.imag_data(), o.imag_data());
0271    eval_subtract(t3, t1, t2);
0272    eval_multiply(t1, result.real_data(), o.imag_data());
0273    eval_multiply(t2, result.imag_data(), o.real_data());
0274    eval_add(t1, t2);
0275    result.real_data() = std::move(t3);
0276    result.imag_data() = std::move(t1);
0277 }
0278 template <class Backend>
0279 inline void eval_divide(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& z)
0280 {
0281    // (a+bi) / (c + di)
0282    using default_ops::eval_add;
0283    using default_ops::eval_divide;
0284    using default_ops::eval_fabs;
0285    using default_ops::eval_is_zero;
0286    using default_ops::eval_multiply;
0287    using default_ops::eval_subtract;
0288    Backend t1, t2;
0289
0290    //
0291    // Backup sign bits for later, so we can fix up
0292    // signed zeros at the end:
0293    //
0294    int a_sign = eval_signbit(result.real_data());
0295    int b_sign = eval_signbit(result.imag_data());
0296    int c_sign = eval_signbit(z.real_data());
0297    int d_sign = eval_signbit(z.imag_data());
0298
0299    if (eval_is_zero(z.imag_data()))
0300    {
0301       eval_divide(result.real_data(), z.real_data());
0302       eval_divide(result.imag_data(), z.real_data());
0303    }
0304    else
0305    {
0306       eval_fabs(t1, z.real_data());
0307       eval_fabs(t2, z.imag_data());
0308       if (t1.compare(t2) < 0)
0309       {
0310          eval_divide(t1, z.real_data(), z.imag_data()); // t1 = c/d
0311          eval_multiply(t2, z.real_data(), t1);
0312          eval_add(t2, z.imag_data()); // denom = c * (c/d) + d
0313          Backend t_real(result.real_data());
0314          // real = (a * (c/d) + b) / (denom)
0315          eval_multiply(result.real_data(), t1);
0316          eval_add(result.real_data(), result.imag_data());
0317          eval_divide(result.real_data(), t2);
0318          // imag = (b * c/d - a) / denom
0319          eval_multiply(result.imag_data(), t1);
0320          eval_subtract(result.imag_data(), t_real);
0321          eval_divide(result.imag_data(), t2);
0322       }
0323       else
0324       {
0325          eval_divide(t1, z.imag_data(), z.real_data()); // t1 = d/c
0326          eval_multiply(t2, z.imag_data(), t1);
0327          eval_add(t2, z.real_data()); // denom = d * d/c + c
0328
0329          Backend r_t(result.real_data());
0330          Backend i_t(result.imag_data());
0331
0332          // real = (b * d/c + a) / denom
0333          eval_multiply(result.real_data(), result.imag_data(), t1);
0334          eval_add(result.real_data(), r_t);
0335          eval_divide(result.real_data(), t2);
0336          // imag = (-a * d/c + b) / denom
0337          eval_multiply(result.imag_data(), r_t, t1);
0338          result.imag_data().negate();
0339          eval_add(result.imag_data(), i_t);
0340          eval_divide(result.imag_data(), t2);
0341       }
0342    }
0343    //
0344    // Finish off by fixing up signed zeros.
0345    //
0346    // This sets the signs "as if" we had evaluated the result using:
0347    //
0348    // real = (ac + bd) / (c^2 + d^2)
0349    // imag = (bc - ad) / (c^2 + d^2)
0350    //
0351    // ie a zero is negative only if the two parts of the numerator
0352    // are both negative and zero.
0353    //
0354    if (eval_is_zero(result.real_data()))
0355    {
0356       int r_sign = eval_signbit(result.real_data());
0357       int r_required = (a_sign != c_sign) && (b_sign != d_sign);
0358       if (r_required != r_sign)
0359          result.real_data().negate();
0360    }
0361    if (eval_is_zero(result.imag_data()))
0362    {
0363       int i_sign = eval_signbit(result.imag_data());
0364       int i_required = (b_sign != c_sign) && (a_sign == d_sign);
0365       if (i_required != i_sign)
0366          result.imag_data().negate();
0367    }
0368 }
0369 template <class Backend, class T>
0370 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_add(complex_adaptor<Backend>& result, const T& scalar)
0371 {
0372    using default_ops::eval_add;
0373    eval_add(result.real_data(), scalar);
0374 }
0375 template <class Backend, class T>
0376 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_subtract(complex_adaptor<Backend>& result, const T& scalar)
0377 {
0378    using default_ops::eval_subtract;
0379    eval_subtract(result.real_data(), scalar);
0380 }
0381 template <class Backend, class T>
0382 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_multiply(complex_adaptor<Backend>& result, const T& scalar)
0383 {
0384    using default_ops::eval_multiply;
0385    eval_multiply(result.real_data(), scalar);
0386    eval_multiply(result.imag_data(), scalar);
0387 }
0388 template <class Backend, class T>
0389 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_divide(complex_adaptor<Backend>& result, const T& scalar)
0390 {
0391    using default_ops::eval_divide;
0392    eval_divide(result.real_data(), scalar);
0393    eval_divide(result.imag_data(), scalar);
0394 }
0395 // Optimised 3 arg versions:
0396 template <class Backend, class T>
0397 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_add(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& a, const T& scalar)
0398 {
0399    using default_ops::eval_add;
0400    eval_add(result.real_data(), a.real_data(), scalar);
0401    result.imag_data() = a.imag_data();
0402 }
0403 template <class Backend, class T>
0404 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_subtract(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& a, const T& scalar)
0405 {
0406    using default_ops::eval_subtract;
0407    eval_subtract(result.real_data(), a.real_data(), scalar);
0408    result.imag_data() = a.imag_data();
0409 }
0410 template <class Backend, class T>
0411 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_multiply(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& a, const T& scalar)
0412 {
0413    using default_ops::eval_multiply;
0414    eval_multiply(result.real_data(), a.real_data(), scalar);
0415    eval_multiply(result.imag_data(), a.imag_data(), scalar);
0416 }
0417 template <class Backend, class T>
0418 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_divide(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& a, const T& scalar)
0419 {
0420    using default_ops::eval_divide;
0421    eval_divide(result.real_data(), a.real_data(), scalar);
0422    eval_divide(result.imag_data(), a.imag_data(), scalar);
0423 }
0424
0425 template <class Backend>
0426 inline bool eval_is_zero(const complex_adaptor<Backend>& val) noexcept
0427 {
0428    using default_ops::eval_is_zero;
0429    return eval_is_zero(val.real_data()) && eval_is_zero(val.imag_data());
0430 }
0431 template <class Backend>
0432 inline int eval_get_sign(const complex_adaptor<Backend>&)
0433 {
0434    static_assert(sizeof(Backend) == UINT_MAX, "Complex numbers have no sign bit."); // designed to always fail
0435    return 0;
0436 }
0437
0438 template <class Result, class Backend>
0439 inline typename std::enable_if< !boost::multiprecision::detail::is_complex<Result>::value>::type eval_convert_to(Result* result, const complex_adaptor<Backend>& val)
0440 {
0441    using default_ops::eval_convert_to;
0442    using default_ops::eval_is_zero;
0443    if (!eval_is_zero(val.imag_data()))
0444    {
0445       BOOST_MP_THROW_EXCEPTION(std::runtime_error("Could not convert imaginary number to scalar."));
0446    }
0447    eval_convert_to(result, val.real_data());
0448 }
0449
0450 template <class Backend, class T>
0451 inline void assign_components(complex_adaptor<Backend>& result, const T& a, const T& b)
0452 {
0453    result.real_data() = a;
0454    result.imag_data() = b;
0455 }
0456
0457 //
0458 // Native non-member operations:
0459 //
0460 template <class Backend>
0461 inline void eval_sqrt(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& val)
0462 {
0463    // Use the following:
0464    // sqrt(z) = (s, zi / 2s)       for zr >= 0
0465    //           (|zi| / 2s, +-s)   for zr <  0
0466    // where s = sqrt{ [ |zr| + sqrt(zr^2 + zi^2) ] / 2 },
0467    // and the +- sign is the same as the sign of zi.
0468    using default_ops::eval_abs;
0469    using default_ops::eval_add;
0470    using default_ops::eval_divide;
0471    using default_ops::eval_get_sign;
0472    using default_ops::eval_is_zero;
0473
0474    if (eval_is_zero(val.imag_data()) && (eval_get_sign(val.real_data()) >= 0))
0475    {
0476       constexpr typename std::tuple_element<0, typename Backend::unsigned_types>::type zero = 0u;
0477       eval_sqrt(result.real_data(), val.real_data());
0478       result.imag_data() = zero;
0479       return;
0480    }
0481
0482    const bool __my_real_part_is_neg(eval_get_sign(val.real_data()) < 0);
0483
0484    Backend __my_real_part_fabs(val.real_data());
0485    if (__my_real_part_is_neg)
0486       __my_real_part_fabs.negate();
0487
0488    Backend t, __my_sqrt_part;
0489    eval_abs(__my_sqrt_part, val);
0490    eval_add(__my_sqrt_part, __my_real_part_fabs);
0491    eval_ldexp(t, __my_sqrt_part, -1);
0492    eval_sqrt(__my_sqrt_part, t);
0493
0494    if (__my_real_part_is_neg == false)
0495    {
0496       eval_ldexp(t, __my_sqrt_part, 1);
0497       eval_divide(result.imag_data(), val.imag_data(), t);
0498       result.real_data() = __my_sqrt_part;
0499    }
0500    else
0501    {
0502       const bool __my_imag_part_is_neg(eval_get_sign(val.imag_data()) < 0);
0503
0504       Backend __my_imag_part_fabs(val.imag_data());
0505       if (__my_imag_part_is_neg)
0506          __my_imag_part_fabs.negate();
0507
0508       eval_ldexp(t, __my_sqrt_part, 1);
0509       eval_divide(result.real_data(), __my_imag_part_fabs, t);
0510       if (__my_imag_part_is_neg)
0511          __my_sqrt_part.negate();
0512       result.imag_data() = __my_sqrt_part;
0513    }
0514 }
0515
0516 template <class Backend>
0517 inline void eval_abs(Backend& result, const complex_adaptor<Backend>& val)
0518 {
0519    Backend t1, t2;
0520    eval_multiply(t1, val.real_data(), val.real_data());
0521    eval_multiply(t2, val.imag_data(), val.imag_data());
0522    eval_add(t1, t2);
0523    eval_sqrt(result, t1);
0524 }
0525
0526 template <class Backend>
0527 inline void eval_pow(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& b, const complex_adaptor<Backend>& e)
0528 {
0529    using default_ops::eval_acos;
0530    using default_ops::eval_cos;
0531    using default_ops::eval_exp;
0532    using default_ops::eval_get_sign;
0533    using default_ops::eval_is_zero;
0534    using default_ops::eval_multiply;
0535    using default_ops::eval_sin;
0536
0537    if (eval_is_zero(e))
0538    {
0539       typename std::tuple_element<0, typename Backend::unsigned_types>::type one(1);
0540       result = one;
0541       return;
0542    }
0543    else if (eval_is_zero(b))
0544    {
0545       if (eval_is_zero(e.real_data()))
0546       {
0547          Backend n          = std::numeric_limits<number<Backend> >::quiet_NaN().backend();
0548          result.real_data() = n;
0549          result.imag_data() = n;
0550       }
0551       else if (eval_get_sign(e.real_data()) < 0)
0552       {
0553          Backend n          = std::numeric_limits<number<Backend> >::infinity().backend();
0554          result.real_data() = n;
0555          typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
0556          if (eval_is_zero(e.imag_data()))
0557             result.imag_data() = zero;
0558          else
0559             result.imag_data() = n;
0560       }
0561       else
0562       {
0563          typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
0564          result = zero;
0565       }
0566       return;
0567    }
0568    complex_adaptor<Backend> t;
0569    eval_log(t, b);
0570    eval_multiply(t, e);
0571    eval_exp(result, t);
0572 }
0573
0574 template <class Backend>
0575 inline void eval_exp(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0576 {
0577    using default_ops::eval_cos;
0578    using default_ops::eval_exp;
0579    using default_ops::eval_is_zero;
0580    using default_ops::eval_multiply;
0581    using default_ops::eval_sin;
0582
0583    if (eval_is_zero(arg.imag_data()))
0584    {
0585       eval_exp(result.real_data(), arg.real_data());
0586       typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
0587       result.imag_data() = zero;
0588       return;
0589    }
0590    eval_cos(result.real_data(), arg.imag_data());
0591    eval_sin(result.imag_data(), arg.imag_data());
0592    Backend e;
0593    eval_exp(e, arg.real_data());
0594    if (eval_is_zero(result.real_data()))
0595       eval_multiply(result.imag_data(), e);
0596    else if (eval_is_zero(result.imag_data()))
0597       eval_multiply(result.real_data(), e);
0598    else
0599       eval_multiply(result, e);
0600 }
0601
0602 template <class Backend>
0603 inline void eval_log(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0604 {
0605    using default_ops::eval_add;
0606    using default_ops::eval_atan2;
0607    using default_ops::eval_get_sign;
0608    using default_ops::eval_is_zero;
0609    using default_ops::eval_log;
0610    using default_ops::eval_multiply;
0611
0612    if (eval_is_zero(arg.imag_data()) && (eval_get_sign(arg.real_data()) >= 0))
0613    {
0614       eval_log(result.real_data(), arg.real_data());
0615       typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
0616       result.imag_data() = zero;
0617       return;
0618    }
0619
0620    Backend t1, t2;
0621    eval_multiply(t1, arg.real_data(), arg.real_data());
0622    eval_multiply(t2, arg.imag_data(), arg.imag_data());
0623    eval_add(t1, t2);
0624    eval_log(t2, t1);
0625    eval_ldexp(result.real_data(), t2, -1);
0626    eval_atan2(result.imag_data(), arg.imag_data(), arg.real_data());
0627 }
0628
0629 template <class Backend>
0630 inline void eval_log10(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0631 {
0632    using default_ops::eval_divide;
0633    using default_ops::eval_log;
0634
0635    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0636
0637    Backend ten;
0638    ten = ui_type(10);
0639    Backend l_ten;
0640    eval_log(l_ten, ten);
0641    eval_log(result, arg);
0642    eval_divide(result, l_ten);
0643 }
0644
0645 template <class Backend>
0646 inline void eval_sin(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0647 {
0648    using default_ops::eval_cos;
0649    using default_ops::eval_cosh;
0650    using default_ops::eval_sin;
0651    using default_ops::eval_sinh;
0652
0653    Backend t1, t2, t3;
0654    eval_sin(t1, arg.real_data());
0655    eval_cosh(t2, arg.imag_data());
0656    eval_multiply(t3, t1, t2);
0657
0658    eval_cos(t1, arg.real_data());
0659    eval_sinh(t2, arg.imag_data());
0660    eval_multiply(result.imag_data(), t1, t2);
0661    result.real_data() = t3;
0662 }
0663
0664 template <class Backend>
0665 inline void eval_cos(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0666 {
0667    using default_ops::eval_cos;
0668    using default_ops::eval_cosh;
0669    using default_ops::eval_sin;
0670    using default_ops::eval_sinh;
0671
0672    Backend t1, t2, t3;
0673    eval_cos(t1, arg.real_data());
0674    eval_cosh(t2, arg.imag_data());
0675    eval_multiply(t3, t1, t2);
0676
0677    eval_sin(t1, arg.real_data());
0678    eval_sinh(t2, arg.imag_data());
0679    eval_multiply(result.imag_data(), t1, t2);
0680    result.imag_data().negate();
0681    result.real_data() = t3;
0682 }
0683
0684 template <class T>
0685 void tanh_imp(const T& r, const T& i, T& r_result, T& i_result)
0686 {
0687    using default_ops::eval_tan;
0688    using default_ops::eval_sinh;
0689    using default_ops::eval_add;
0690    using default_ops::eval_fpclassify;
0691    using default_ops::eval_get_sign;
0692
0693    using ui_type = typename std::tuple_element<0, typename T::unsigned_types>::type;
0694    ui_type one(1);
0695    //
0696    // Set:
0697    // t = tan(i);
0698    // s = sinh(r);
0699    // b = s * (1 + t^2);
0700    // d = 1 + b * s;
0701    //
0702    T t, s, b, d;
0703    eval_tan(t, i);
0704    eval_sinh(s, r);
0705    eval_multiply(d, t, t);
0706    eval_add(d, one);
0707    eval_multiply(b, d, s);
0708    eval_multiply(d, b, s);
0709    eval_add(d, one);
0710
0711    if (eval_fpclassify(d) == FP_INFINITE)
0712    {
0713       r_result = one;
0714       if (eval_get_sign(s) < 0)
0715          r_result.negate();
0716       //
0717       // Imaginary part is a signed zero:
0718       //
0719       ui_type zero(0);
0720       i_result = zero;
0721       if (eval_get_sign(t) < 0)
0722          i_result.negate();
0723    }
0724    //
0725    // Real part is sqrt(1 + s^2) * b / d;
0726    // Imaginary part is t / d;
0727    //
0728    eval_divide(i_result, t, d);
0729    //
0730    // variable t is now spare, as is r_result.
0731    //
0732    eval_multiply(t, s, s);
0733    eval_add(t, one);
0734    eval_sqrt(r_result, t);
0735    eval_multiply(t, r_result, b);
0736    eval_divide(r_result, t, d);
0737 }
0738
0739 template <class Backend>
0740 inline void eval_tanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0741 {
0742    tanh_imp(arg.real_data(), arg.imag_data(), result.real_data(), result.imag_data());
0743 }
0744 template <class Backend>
0745 inline void eval_tan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0746 {
0747    Backend t(arg.imag_data());
0748    t.negate();
0749    tanh_imp(t, arg.real_data(), result.imag_data(), result.real_data());
0750    result.imag_data().negate();
0751 }
0752
0753 template <class Backend>
0754 inline void eval_asin(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0755 {
0756    using default_ops::eval_add;
0757    using default_ops::eval_multiply;
0758
0759    if (eval_is_zero(arg))
0760    {
0761       result = arg;
0762       return;
0763    }
0764
0765    complex_adaptor<Backend> t1, t2;
0766    assign_components(t1, arg.imag_data(), arg.real_data());
0767    t1.real_data().negate();
0768    eval_asinh(t2, t1);
0769
0770    assign_components(result, t2.imag_data(), t2.real_data());
0771    result.imag_data().negate();
0772 }
0773
0774 template <class Backend>
0775 inline void eval_acos(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0776 {
0777    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0778
0779    using default_ops::eval_asin;
0780
0781    Backend half_pi, t1;
0782    t1 = static_cast<ui_type>(1u);
0783    eval_asin(half_pi, t1);
0784    eval_asin(result, arg);
0785    result.negate();
0786    eval_add(result.real_data(), half_pi);
0787 }
0788
0789 template <class Backend>
0790 inline void eval_atan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0791 {
0792    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0793    ui_type                                                             one = (ui_type)1u;
0794
0795    using default_ops::eval_add;
0796    using default_ops::eval_is_zero;
0797    using default_ops::eval_log;
0798    using default_ops::eval_subtract;
0799
0800    complex_adaptor<Backend> __my_z_times_i, t1, t2, t3;
0801    assign_components(__my_z_times_i, arg.imag_data(), arg.real_data());
0802    __my_z_times_i.real_data().negate();
0803
0804    eval_add(t1, __my_z_times_i, one);
0805    eval_log(t2, t1);
0806    eval_subtract(t1, one, __my_z_times_i);
0807    eval_log(t3, t1);
0808    eval_subtract(t1, t3, t2);
0809
0810    eval_ldexp(result.real_data(), t1.imag_data(), -1);
0811    eval_ldexp(result.imag_data(), t1.real_data(), -1);
0812    if (!eval_is_zero(result.real_data()))
0813       result.real_data().negate();
0814 }
0815
0816 template <class Backend>
0817 inline void eval_sinh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0818 {
0819    using default_ops::eval_cos;
0820    using default_ops::eval_cosh;
0821    using default_ops::eval_sin;
0822    using default_ops::eval_sinh;
0823
0824    Backend t1, t2, t3;
0825    eval_cos(t1, arg.imag_data());
0826    eval_sinh(t2, arg.real_data());
0827    eval_multiply(t3, t1, t2);
0828
0829    eval_cosh(t1, arg.real_data());
0830    eval_sin(t2, arg.imag_data());
0831    eval_multiply(result.imag_data(), t1, t2);
0832    result.real_data() = t3;
0833 }
0834
0835 template <class Backend>
0836 inline void eval_cosh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0837 {
0838    using default_ops::eval_cos;
0839    using default_ops::eval_cosh;
0840    using default_ops::eval_sin;
0841    using default_ops::eval_sinh;
0842
0843    Backend t1, t2, t3;
0844    eval_cos(t1, arg.imag_data());
0845    eval_cosh(t2, arg.real_data());
0846    eval_multiply(t3, t1, t2);
0847
0848    eval_sin(t1, arg.imag_data());
0849    eval_sinh(t2, arg.real_data());
0850    eval_multiply(result.imag_data(), t1, t2);
0851    result.real_data() = t3;
0852 }
0853
0854 template <class Backend>
0855 inline void eval_asinh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0856 {
0857    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0858    ui_type                                                             one = (ui_type)1u;
0859
0860    using default_ops::eval_add;
0861    using default_ops::eval_log;
0862    using default_ops::eval_multiply;
0863
0864    complex_adaptor<Backend> t1, t2;
0865    eval_multiply(t1, arg, arg);
0866    eval_add(t1, one);
0867    eval_sqrt(t2, t1);
0868    eval_add(t2, arg);
0869    eval_log(result, t2);
0870 }
0871
0872 template <class Backend>
0873 inline void eval_acosh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0874 {
0875    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0876    ui_type                                                             one = (ui_type)1u;
0877
0878    using default_ops::eval_add;
0879    using default_ops::eval_divide;
0880    using default_ops::eval_log;
0881    using default_ops::eval_multiply;
0882    using default_ops::eval_subtract;
0883
0884    complex_adaptor<Backend> __my_zp(arg);
0885    eval_add(__my_zp.real_data(), one);
0886    complex_adaptor<Backend> __my_zm(arg);
0887    eval_subtract(__my_zm.real_data(), one);
0888
0889    complex_adaptor<Backend> t1, t2;
0890    eval_divide(t1, __my_zm, __my_zp);
0891    eval_sqrt(t2, t1);
0892    eval_multiply(t2, __my_zp);
0893    eval_add(t2, arg);
0894    eval_log(result, t2);
0895 }
0896
0897 template <class Backend>
0898 inline void eval_atanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0899 {
0900    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0901    ui_type                                                             one = (ui_type)1u;
0902
0903    using default_ops::eval_add;
0904    using default_ops::eval_divide;
0905    using default_ops::eval_log;
0906    using default_ops::eval_multiply;
0907    using default_ops::eval_subtract;
0908
0909    complex_adaptor<Backend> t1, t2, t3;
0910    eval_add(t1, arg, one);
0911    eval_log(t2, t1);
0912    eval_subtract(t1, one, arg);
0913    eval_log(t3, t1);
0914    eval_subtract(t2, t3);
0915
0916    eval_ldexp(result.real_data(), t2.real_data(), -1);
0917    eval_ldexp(result.imag_data(), t2.imag_data(), -1);
0918 }
0919
0920 template <class Backend>
0921 inline void eval_conj(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0922 {
0923    result = arg;
0924    result.imag_data().negate();
0925 }
0926
0927 template <class Backend>
0928 inline void eval_proj(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0929 {
0930    using default_ops::eval_get_sign;
0931
0932    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0933    ui_type                                                             zero = (ui_type)0u;
0934
0935    int c1 = eval_fpclassify(arg.real_data());
0936    int c2 = eval_fpclassify(arg.imag_data());
0937    if (c1 == FP_INFINITE)
0938    {
0939       result.real_data() = arg.real_data();
0940       if (eval_get_sign(result.real_data()) < 0)
0941          result.real_data().negate();
0942       result.imag_data() = zero;
0943       if (eval_get_sign(arg.imag_data()) < 0)
0944          result.imag_data().negate();
0945    }
0946    else if (c2 == FP_INFINITE)
0947    {
0948       result.real_data() = arg.imag_data();
0949       if (eval_get_sign(result.real_data()) < 0)
0950          result.real_data().negate();
0951       result.imag_data() = zero;
0952       if (eval_get_sign(arg.imag_data()) < 0)
0953          result.imag_data().negate();
0954    }
0955    else
0956       result = arg;
0957 }
0958
0959 template <class Backend>
0960 inline void eval_real(Backend& result, const complex_adaptor<Backend>& arg)
0961 {
0962    result = arg.real_data();
0963 }
0964 template <class Backend>
0965 inline void eval_imag(Backend& result, const complex_adaptor<Backend>& arg)
0966 {
0967    result = arg.imag_data();
0968 }
0969
0970 template <class Backend, class T>
0971 inline void eval_set_imag(complex_adaptor<Backend>& result, const T& arg)
0972 {
0973    result.imag_data() = arg;
0974 }
0975
0976 template <class Backend, class T>
0977 inline void eval_set_real(complex_adaptor<Backend>& result, const T& arg)
0978 {
0979    result.real_data() = arg;
0980 }
0981
0982 template <class Backend>
0983 inline std::size_t hash_value(const complex_adaptor<Backend>& val)
0984 {
0985    std::size_t result  = hash_value(val.real_data());
0986    std::size_t result2 = hash_value(val.imag_data());
0987    boost::multiprecision::detail::hash_combine(result, result2);
0988    return result;
0989 }
0990
0991 } // namespace backends
0992
0993 template <class Backend>
0994 struct number_category<complex_adaptor<Backend> > : public std::integral_constant<int, boost::multiprecision::number_kind_complex>
0995 {};
0996
0997 template <class Backend, expression_template_option ExpressionTemplates>
0998 struct component_type<number<complex_adaptor<Backend>, ExpressionTemplates> >
0999 {
1000    using type = number<Backend, ExpressionTemplates>;
1001 };
1002
1003 template <class Backend, expression_template_option ExpressionTemplates>
1004 struct complex_result_from_scalar<number<Backend, ExpressionTemplates> >
1005 {
1006    using type = number<complex_adaptor<Backend>, ExpressionTemplates>;
1007 };
1008
1009 namespace detail {
1010    template <class Backend>
1011    struct is_variable_precision<complex_adaptor<Backend> > : public is_variable_precision<Backend>
1012    {};
1013 #ifdef BOOST_HAS_INT128
1014    template <class Backend>
1015    struct is_convertible_arithmetic<int128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<int128_type, Backend>
1016    {};
1017    template <class Backend>
1018    struct is_convertible_arithmetic<uint128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<uint128_type, Backend>
1019    {};
1020 #endif
1021 #ifdef BOOST_HAS_FLOAT128
1022    template <class Backend>
1023    struct is_convertible_arithmetic<float128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<float128_type, Backend>
1024    {};
1025 #endif
1026    } // namespace detail
1027
1028
1029
1030 template <class Backend, expression_template_option ExpressionTemplates>
1031 struct complex_result_from_scalar<number<backends::debug_adaptor<Backend>, ExpressionTemplates> >
1032 {
1033    using type = number<backends::debug_adaptor<complex_adaptor<Backend> >, ExpressionTemplates>;
1034 };
1035
1036 template <class Backend, expression_template_option ExpressionTemplates>
1037 struct complex_result_from_scalar<number<backends::logged_adaptor<Backend>, ExpressionTemplates> >
1038 {
1039    using type = number<backends::logged_adaptor<complex_adaptor<Backend> >, ExpressionTemplates>;
1040 };
1041
1042 }
1043
1044 } // namespace boost::multiprecision
1045
1046 #endif