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 template <class Backend, class T>
0198 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
0199 {
0200    return a.compare(b) == 0;
0201 }
0202
0203 template <class Backend>
0204 inline void eval_add(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
0205 {
0206    eval_add(result.real_data(), o.real_data());
0207    eval_add(result.imag_data(), o.imag_data());
0208 }
0209 template <class Backend>
0210 inline void eval_subtract(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
0211 {
0212    eval_subtract(result.real_data(), o.real_data());
0213    eval_subtract(result.imag_data(), o.imag_data());
0214 }
0215 template <class Backend>
0216 inline void eval_multiply(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
0217 {
0218    Backend t1, t2, t3;
0219    eval_multiply(t1, result.real_data(), o.real_data());
0220    eval_multiply(t2, result.imag_data(), o.imag_data());
0221    eval_subtract(t3, t1, t2);
0222    eval_multiply(t1, result.real_data(), o.imag_data());
0223    eval_multiply(t2, result.imag_data(), o.real_data());
0224    eval_add(t1, t2);
0225    result.real_data() = std::move(t3);
0226    result.imag_data() = std::move(t1);
0227 }
0228 template <class Backend>
0229 inline void eval_divide(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& z)
0230 {
0231    // (a+bi) / (c + di)
0232    using default_ops::eval_add;
0233    using default_ops::eval_divide;
0234    using default_ops::eval_fabs;
0235    using default_ops::eval_is_zero;
0236    using default_ops::eval_multiply;
0237    using default_ops::eval_subtract;
0238    Backend t1, t2;
0239
0240    //
0241    // Backup sign bits for later, so we can fix up
0242    // signed zeros at the end:
0243    //
0244    int a_sign = eval_signbit(result.real_data());
0245    int b_sign = eval_signbit(result.imag_data());
0246    int c_sign = eval_signbit(z.real_data());
0247    int d_sign = eval_signbit(z.imag_data());
0248
0249    if (eval_is_zero(z.imag_data()))
0250    {
0251       eval_divide(result.real_data(), z.real_data());
0252       eval_divide(result.imag_data(), z.real_data());
0253    }
0254    else
0255    {
0256       eval_fabs(t1, z.real_data());
0257       eval_fabs(t2, z.imag_data());
0258       if (t1.compare(t2) < 0)
0259       {
0260          eval_divide(t1, z.real_data(), z.imag_data()); // t1 = c/d
0261          eval_multiply(t2, z.real_data(), t1);
0262          eval_add(t2, z.imag_data()); // denom = c * (c/d) + d
0263          Backend t_real(result.real_data());
0264          // real = (a * (c/d) + b) / (denom)
0265          eval_multiply(result.real_data(), t1);
0266          eval_add(result.real_data(), result.imag_data());
0267          eval_divide(result.real_data(), t2);
0268          // imag = (b * c/d - a) / denom
0269          eval_multiply(result.imag_data(), t1);
0270          eval_subtract(result.imag_data(), t_real);
0271          eval_divide(result.imag_data(), t2);
0272       }
0273       else
0274       {
0275          eval_divide(t1, z.imag_data(), z.real_data()); // t1 = d/c
0276          eval_multiply(t2, z.imag_data(), t1);
0277          eval_add(t2, z.real_data()); // denom = d * d/c + c
0278
0279          Backend r_t(result.real_data());
0280          Backend i_t(result.imag_data());
0281
0282          // real = (b * d/c + a) / denom
0283          eval_multiply(result.real_data(), result.imag_data(), t1);
0284          eval_add(result.real_data(), r_t);
0285          eval_divide(result.real_data(), t2);
0286          // imag = (-a * d/c + b) / denom
0287          eval_multiply(result.imag_data(), r_t, t1);
0288          result.imag_data().negate();
0289          eval_add(result.imag_data(), i_t);
0290          eval_divide(result.imag_data(), t2);
0291       }
0292    }
0293    //
0294    // Finish off by fixing up signed zeros.
0295    //
0296    // This sets the signs "as if" we had evaluated the result using:
0297    //
0298    // real = (ac + bd) / (c^2 + d^2)
0299    // imag = (bc - ad) / (c^2 + d^2)
0300    //
0301    // ie a zero is negative only if the two parts of the numerator
0302    // are both negative and zero.
0303    //
0304    if (eval_is_zero(result.real_data()))
0305    {
0306       int r_sign = eval_signbit(result.real_data());
0307       int r_required = (a_sign != c_sign) && (b_sign != d_sign);
0308       if (r_required != r_sign)
0309          result.real_data().negate();
0310    }
0311    if (eval_is_zero(result.imag_data()))
0312    {
0313       int i_sign = eval_signbit(result.imag_data());
0314       int i_required = (b_sign != c_sign) && (a_sign == d_sign);
0315       if (i_required != i_sign)
0316          result.imag_data().negate();
0317    }
0318 }
0319 template <class Backend, class T>
0320 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_add(complex_adaptor<Backend>& result, const T& scalar)
0321 {
0322    using default_ops::eval_add;
0323    eval_add(result.real_data(), scalar);
0324 }
0325 template <class Backend, class T>
0326 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_subtract(complex_adaptor<Backend>& result, const T& scalar)
0327 {
0328    using default_ops::eval_subtract;
0329    eval_subtract(result.real_data(), scalar);
0330 }
0331 template <class Backend, class T>
0332 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_multiply(complex_adaptor<Backend>& result, const T& scalar)
0333 {
0334    using default_ops::eval_multiply;
0335    eval_multiply(result.real_data(), scalar);
0336    eval_multiply(result.imag_data(), scalar);
0337 }
0338 template <class Backend, class T>
0339 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_divide(complex_adaptor<Backend>& result, const T& scalar)
0340 {
0341    using default_ops::eval_divide;
0342    eval_divide(result.real_data(), scalar);
0343    eval_divide(result.imag_data(), scalar);
0344 }
0345 // Optimised 3 arg versions:
0346 template <class Backend, class T>
0347 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)
0348 {
0349    using default_ops::eval_add;
0350    eval_add(result.real_data(), a.real_data(), scalar);
0351    result.imag_data() = a.imag_data();
0352 }
0353 template <class Backend, class T>
0354 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)
0355 {
0356    using default_ops::eval_subtract;
0357    eval_subtract(result.real_data(), a.real_data(), scalar);
0358    result.imag_data() = a.imag_data();
0359 }
0360 template <class Backend, class T>
0361 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)
0362 {
0363    using default_ops::eval_multiply;
0364    eval_multiply(result.real_data(), a.real_data(), scalar);
0365    eval_multiply(result.imag_data(), a.imag_data(), scalar);
0366 }
0367 template <class Backend, class T>
0368 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)
0369 {
0370    using default_ops::eval_divide;
0371    eval_divide(result.real_data(), a.real_data(), scalar);
0372    eval_divide(result.imag_data(), a.imag_data(), scalar);
0373 }
0374
0375 template <class Backend>
0376 inline bool eval_is_zero(const complex_adaptor<Backend>& val) noexcept
0377 {
0378    using default_ops::eval_is_zero;
0379    return eval_is_zero(val.real_data()) && eval_is_zero(val.imag_data());
0380 }
0381 template <class Backend>
0382 inline int eval_get_sign(const complex_adaptor<Backend>&)
0383 {
0384    static_assert(sizeof(Backend) == UINT_MAX, "Complex numbers have no sign bit."); // designed to always fail
0385    return 0;
0386 }
0387
0388 template <class Result, class Backend>
0389 inline typename std::enable_if< !boost::multiprecision::detail::is_complex<Result>::value>::type eval_convert_to(Result* result, const complex_adaptor<Backend>& val)
0390 {
0391    using default_ops::eval_convert_to;
0392    using default_ops::eval_is_zero;
0393    if (!eval_is_zero(val.imag_data()))
0394    {
0395       BOOST_MP_THROW_EXCEPTION(std::runtime_error("Could not convert imaginary number to scalar."));
0396    }
0397    eval_convert_to(result, val.real_data());
0398 }
0399
0400 template <class Backend, class T>
0401 inline void assign_components(complex_adaptor<Backend>& result, const T& a, const T& b)
0402 {
0403    result.real_data() = a;
0404    result.imag_data() = b;
0405 }
0406
0407 //
0408 // Native non-member operations:
0409 //
0410 template <class Backend>
0411 inline void eval_sqrt(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& val)
0412 {
0413    // Use the following:
0414    // sqrt(z) = (s, zi / 2s)       for zr >= 0
0415    //           (|zi| / 2s, +-s)   for zr <  0
0416    // where s = sqrt{ [ |zr| + sqrt(zr^2 + zi^2) ] / 2 },
0417    // and the +- sign is the same as the sign of zi.
0418    using default_ops::eval_abs;
0419    using default_ops::eval_add;
0420    using default_ops::eval_divide;
0421    using default_ops::eval_get_sign;
0422    using default_ops::eval_is_zero;
0423
0424    if (eval_is_zero(val.imag_data()) && (eval_get_sign(val.real_data()) >= 0))
0425    {
0426       constexpr typename std::tuple_element<0, typename Backend::unsigned_types>::type zero = 0u;
0427       eval_sqrt(result.real_data(), val.real_data());
0428       result.imag_data() = zero;
0429       return;
0430    }
0431
0432    const bool __my_real_part_is_neg(eval_get_sign(val.real_data()) < 0);
0433
0434    Backend __my_real_part_fabs(val.real_data());
0435    if (__my_real_part_is_neg)
0436       __my_real_part_fabs.negate();
0437
0438    Backend t, __my_sqrt_part;
0439    eval_abs(__my_sqrt_part, val);
0440    eval_add(__my_sqrt_part, __my_real_part_fabs);
0441    eval_ldexp(t, __my_sqrt_part, -1);
0442    eval_sqrt(__my_sqrt_part, t);
0443
0444    if (__my_real_part_is_neg == false)
0445    {
0446       eval_ldexp(t, __my_sqrt_part, 1);
0447       eval_divide(result.imag_data(), val.imag_data(), t);
0448       result.real_data() = __my_sqrt_part;
0449    }
0450    else
0451    {
0452       const bool __my_imag_part_is_neg(eval_get_sign(val.imag_data()) < 0);
0453
0454       Backend __my_imag_part_fabs(val.imag_data());
0455       if (__my_imag_part_is_neg)
0456          __my_imag_part_fabs.negate();
0457
0458       eval_ldexp(t, __my_sqrt_part, 1);
0459       eval_divide(result.real_data(), __my_imag_part_fabs, t);
0460       if (__my_imag_part_is_neg)
0461          __my_sqrt_part.negate();
0462       result.imag_data() = __my_sqrt_part;
0463    }
0464 }
0465
0466 template <class Backend>
0467 inline void eval_abs(Backend& result, const complex_adaptor<Backend>& val)
0468 {
0469    Backend t1, t2;
0470    eval_multiply(t1, val.real_data(), val.real_data());
0471    eval_multiply(t2, val.imag_data(), val.imag_data());
0472    eval_add(t1, t2);
0473    eval_sqrt(result, t1);
0474 }
0475
0476 template <class Backend>
0477 inline void eval_pow(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& b, const complex_adaptor<Backend>& e)
0478 {
0479    using default_ops::eval_acos;
0480    using default_ops::eval_cos;
0481    using default_ops::eval_exp;
0482    using default_ops::eval_get_sign;
0483    using default_ops::eval_is_zero;
0484    using default_ops::eval_multiply;
0485    using default_ops::eval_sin;
0486
0487    if (eval_is_zero(e))
0488    {
0489       typename std::tuple_element<0, typename Backend::unsigned_types>::type one(1);
0490       result = one;
0491       return;
0492    }
0493    else if (eval_is_zero(b))
0494    {
0495       if (eval_is_zero(e.real_data()))
0496       {
0497          Backend n          = std::numeric_limits<number<Backend> >::quiet_NaN().backend();
0498          result.real_data() = n;
0499          result.imag_data() = n;
0500       }
0501       else if (eval_get_sign(e.real_data()) < 0)
0502       {
0503          Backend n          = std::numeric_limits<number<Backend> >::infinity().backend();
0504          result.real_data() = n;
0505          typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
0506          if (eval_is_zero(e.imag_data()))
0507             result.imag_data() = zero;
0508          else
0509             result.imag_data() = n;
0510       }
0511       else
0512       {
0513          typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
0514          result = zero;
0515       }
0516       return;
0517    }
0518    complex_adaptor<Backend> t;
0519    eval_log(t, b);
0520    eval_multiply(t, e);
0521    eval_exp(result, t);
0522 }
0523
0524 template <class Backend>
0525 inline void eval_exp(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0526 {
0527    using default_ops::eval_cos;
0528    using default_ops::eval_exp;
0529    using default_ops::eval_is_zero;
0530    using default_ops::eval_multiply;
0531    using default_ops::eval_sin;
0532
0533    if (eval_is_zero(arg.imag_data()))
0534    {
0535       eval_exp(result.real_data(), arg.real_data());
0536       typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
0537       result.imag_data() = zero;
0538       return;
0539    }
0540    eval_cos(result.real_data(), arg.imag_data());
0541    eval_sin(result.imag_data(), arg.imag_data());
0542    Backend e;
0543    eval_exp(e, arg.real_data());
0544    if (eval_is_zero(result.real_data()))
0545       eval_multiply(result.imag_data(), e);
0546    else if (eval_is_zero(result.imag_data()))
0547       eval_multiply(result.real_data(), e);
0548    else
0549       eval_multiply(result, e);
0550 }
0551
0552 template <class Backend>
0553 inline void eval_log(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0554 {
0555    using default_ops::eval_add;
0556    using default_ops::eval_atan2;
0557    using default_ops::eval_get_sign;
0558    using default_ops::eval_is_zero;
0559    using default_ops::eval_log;
0560    using default_ops::eval_multiply;
0561
0562    if (eval_is_zero(arg.imag_data()) && (eval_get_sign(arg.real_data()) >= 0))
0563    {
0564       eval_log(result.real_data(), arg.real_data());
0565       typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
0566       result.imag_data() = zero;
0567       return;
0568    }
0569
0570    Backend t1, t2;
0571    eval_multiply(t1, arg.real_data(), arg.real_data());
0572    eval_multiply(t2, arg.imag_data(), arg.imag_data());
0573    eval_add(t1, t2);
0574    eval_log(t2, t1);
0575    eval_ldexp(result.real_data(), t2, -1);
0576    eval_atan2(result.imag_data(), arg.imag_data(), arg.real_data());
0577 }
0578
0579 template <class Backend>
0580 inline void eval_log10(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0581 {
0582    using default_ops::eval_divide;
0583    using default_ops::eval_log;
0584
0585    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0586
0587    Backend ten;
0588    ten = ui_type(10);
0589    Backend l_ten;
0590    eval_log(l_ten, ten);
0591    eval_log(result, arg);
0592    eval_divide(result, l_ten);
0593 }
0594
0595 template <class Backend>
0596 inline void eval_sin(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0597 {
0598    using default_ops::eval_cos;
0599    using default_ops::eval_cosh;
0600    using default_ops::eval_sin;
0601    using default_ops::eval_sinh;
0602
0603    Backend t1, t2, t3;
0604    eval_sin(t1, arg.real_data());
0605    eval_cosh(t2, arg.imag_data());
0606    eval_multiply(t3, t1, t2);
0607
0608    eval_cos(t1, arg.real_data());
0609    eval_sinh(t2, arg.imag_data());
0610    eval_multiply(result.imag_data(), t1, t2);
0611    result.real_data() = t3;
0612 }
0613
0614 template <class Backend>
0615 inline void eval_cos(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0616 {
0617    using default_ops::eval_cos;
0618    using default_ops::eval_cosh;
0619    using default_ops::eval_sin;
0620    using default_ops::eval_sinh;
0621
0622    Backend t1, t2, t3;
0623    eval_cos(t1, arg.real_data());
0624    eval_cosh(t2, arg.imag_data());
0625    eval_multiply(t3, t1, t2);
0626
0627    eval_sin(t1, arg.real_data());
0628    eval_sinh(t2, arg.imag_data());
0629    eval_multiply(result.imag_data(), t1, t2);
0630    result.imag_data().negate();
0631    result.real_data() = t3;
0632 }
0633
0634 template <class T>
0635 void tanh_imp(const T& r, const T& i, T& r_result, T& i_result)
0636 {
0637    using default_ops::eval_tan;
0638    using default_ops::eval_sinh;
0639    using default_ops::eval_add;
0640    using default_ops::eval_fpclassify;
0641    using default_ops::eval_get_sign;
0642
0643    using ui_type = typename std::tuple_element<0, typename T::unsigned_types>::type;
0644    ui_type one(1);
0645    //
0646    // Set:
0647    // t = tan(i);
0648    // s = sinh(r);
0649    // b = s * (1 + t^2);
0650    // d = 1 + b * s;
0651    //
0652    T t, s, b, d;
0653    eval_tan(t, i);
0654    eval_sinh(s, r);
0655    eval_multiply(d, t, t);
0656    eval_add(d, one);
0657    eval_multiply(b, d, s);
0658    eval_multiply(d, b, s);
0659    eval_add(d, one);
0660
0661    if (eval_fpclassify(d) == FP_INFINITE)
0662    {
0663       r_result = one;
0664       if (eval_get_sign(s) < 0)
0665          r_result.negate();
0666       //
0667       // Imaginary part is a signed zero:
0668       //
0669       ui_type zero(0);
0670       i_result = zero;
0671       if (eval_get_sign(t) < 0)
0672          i_result.negate();
0673    }
0674    //
0675    // Real part is sqrt(1 + s^2) * b / d;
0676    // Imaginary part is t / d;
0677    //
0678    eval_divide(i_result, t, d);
0679    //
0680    // variable t is now spare, as is r_result.
0681    //
0682    eval_multiply(t, s, s);
0683    eval_add(t, one);
0684    eval_sqrt(r_result, t);
0685    eval_multiply(t, r_result, b);
0686    eval_divide(r_result, t, d);
0687 }
0688
0689 template <class Backend>
0690 inline void eval_tanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0691 {
0692    tanh_imp(arg.real_data(), arg.imag_data(), result.real_data(), result.imag_data());
0693 }
0694 template <class Backend>
0695 inline void eval_tan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0696 {
0697    Backend t(arg.imag_data());
0698    t.negate();
0699    tanh_imp(t, arg.real_data(), result.imag_data(), result.real_data());
0700    result.imag_data().negate();
0701 }
0702
0703 template <class Backend>
0704 inline void eval_asin(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0705 {
0706    using default_ops::eval_add;
0707    using default_ops::eval_multiply;
0708
0709    if (eval_is_zero(arg))
0710    {
0711       result = arg;
0712       return;
0713    }
0714
0715    complex_adaptor<Backend> t1, t2;
0716    assign_components(t1, arg.imag_data(), arg.real_data());
0717    t1.real_data().negate();
0718    eval_asinh(t2, t1);
0719
0720    assign_components(result, t2.imag_data(), t2.real_data());
0721    result.imag_data().negate();
0722 }
0723
0724 template <class Backend>
0725 inline void eval_acos(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0726 {
0727    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0728
0729    using default_ops::eval_asin;
0730
0731    Backend half_pi, t1;
0732    t1 = static_cast<ui_type>(1u);
0733    eval_asin(half_pi, t1);
0734    eval_asin(result, arg);
0735    result.negate();
0736    eval_add(result.real_data(), half_pi);
0737 }
0738
0739 template <class Backend>
0740 inline void eval_atan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0741 {
0742    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0743    ui_type                                                             one = (ui_type)1u;
0744
0745    using default_ops::eval_add;
0746    using default_ops::eval_is_zero;
0747    using default_ops::eval_log;
0748    using default_ops::eval_subtract;
0749
0750    complex_adaptor<Backend> __my_z_times_i, t1, t2, t3;
0751    assign_components(__my_z_times_i, arg.imag_data(), arg.real_data());
0752    __my_z_times_i.real_data().negate();
0753
0754    eval_add(t1, __my_z_times_i, one);
0755    eval_log(t2, t1);
0756    eval_subtract(t1, one, __my_z_times_i);
0757    eval_log(t3, t1);
0758    eval_subtract(t1, t3, t2);
0759
0760    eval_ldexp(result.real_data(), t1.imag_data(), -1);
0761    eval_ldexp(result.imag_data(), t1.real_data(), -1);
0762    if (!eval_is_zero(result.real_data()))
0763       result.real_data().negate();
0764 }
0765
0766 template <class Backend>
0767 inline void eval_sinh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0768 {
0769    using default_ops::eval_cos;
0770    using default_ops::eval_cosh;
0771    using default_ops::eval_sin;
0772    using default_ops::eval_sinh;
0773
0774    Backend t1, t2, t3;
0775    eval_cos(t1, arg.imag_data());
0776    eval_sinh(t2, arg.real_data());
0777    eval_multiply(t3, t1, t2);
0778
0779    eval_cosh(t1, arg.real_data());
0780    eval_sin(t2, arg.imag_data());
0781    eval_multiply(result.imag_data(), t1, t2);
0782    result.real_data() = t3;
0783 }
0784
0785 template <class Backend>
0786 inline void eval_cosh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0787 {
0788    using default_ops::eval_cos;
0789    using default_ops::eval_cosh;
0790    using default_ops::eval_sin;
0791    using default_ops::eval_sinh;
0792
0793    Backend t1, t2, t3;
0794    eval_cos(t1, arg.imag_data());
0795    eval_cosh(t2, arg.real_data());
0796    eval_multiply(t3, t1, t2);
0797
0798    eval_sin(t1, arg.imag_data());
0799    eval_sinh(t2, arg.real_data());
0800    eval_multiply(result.imag_data(), t1, t2);
0801    result.real_data() = t3;
0802 }
0803
0804 template <class Backend>
0805 inline void eval_asinh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0806 {
0807    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0808    ui_type                                                             one = (ui_type)1u;
0809
0810    using default_ops::eval_add;
0811    using default_ops::eval_log;
0812    using default_ops::eval_multiply;
0813
0814    complex_adaptor<Backend> t1, t2;
0815    eval_multiply(t1, arg, arg);
0816    eval_add(t1, one);
0817    eval_sqrt(t2, t1);
0818    eval_add(t2, arg);
0819    eval_log(result, t2);
0820 }
0821
0822 template <class Backend>
0823 inline void eval_acosh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0824 {
0825    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0826    ui_type                                                             one = (ui_type)1u;
0827
0828    using default_ops::eval_add;
0829    using default_ops::eval_divide;
0830    using default_ops::eval_log;
0831    using default_ops::eval_multiply;
0832    using default_ops::eval_subtract;
0833
0834    complex_adaptor<Backend> __my_zp(arg);
0835    eval_add(__my_zp.real_data(), one);
0836    complex_adaptor<Backend> __my_zm(arg);
0837    eval_subtract(__my_zm.real_data(), one);
0838
0839    complex_adaptor<Backend> t1, t2;
0840    eval_divide(t1, __my_zm, __my_zp);
0841    eval_sqrt(t2, t1);
0842    eval_multiply(t2, __my_zp);
0843    eval_add(t2, arg);
0844    eval_log(result, t2);
0845 }
0846
0847 template <class Backend>
0848 inline void eval_atanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0849 {
0850    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0851    ui_type                                                             one = (ui_type)1u;
0852
0853    using default_ops::eval_add;
0854    using default_ops::eval_divide;
0855    using default_ops::eval_log;
0856    using default_ops::eval_multiply;
0857    using default_ops::eval_subtract;
0858
0859    complex_adaptor<Backend> t1, t2, t3;
0860    eval_add(t1, arg, one);
0861    eval_log(t2, t1);
0862    eval_subtract(t1, one, arg);
0863    eval_log(t3, t1);
0864    eval_subtract(t2, t3);
0865
0866    eval_ldexp(result.real_data(), t2.real_data(), -1);
0867    eval_ldexp(result.imag_data(), t2.imag_data(), -1);
0868 }
0869
0870 template <class Backend>
0871 inline void eval_conj(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0872 {
0873    result = arg;
0874    result.imag_data().negate();
0875 }
0876
0877 template <class Backend>
0878 inline void eval_proj(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
0879 {
0880    using default_ops::eval_get_sign;
0881
0882    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
0883    ui_type                                                             zero = (ui_type)0u;
0884
0885    int c1 = eval_fpclassify(arg.real_data());
0886    int c2 = eval_fpclassify(arg.imag_data());
0887    if (c1 == FP_INFINITE)
0888    {
0889       result.real_data() = arg.real_data();
0890       if (eval_get_sign(result.real_data()) < 0)
0891          result.real_data().negate();
0892       result.imag_data() = zero;
0893       if (eval_get_sign(arg.imag_data()) < 0)
0894          result.imag_data().negate();
0895    }
0896    else if (c2 == FP_INFINITE)
0897    {
0898       result.real_data() = arg.imag_data();
0899       if (eval_get_sign(result.real_data()) < 0)
0900          result.real_data().negate();
0901       result.imag_data() = zero;
0902       if (eval_get_sign(arg.imag_data()) < 0)
0903          result.imag_data().negate();
0904    }
0905    else
0906       result = arg;
0907 }
0908
0909 template <class Backend>
0910 inline void eval_real(Backend& result, const complex_adaptor<Backend>& arg)
0911 {
0912    result = arg.real_data();
0913 }
0914 template <class Backend>
0915 inline void eval_imag(Backend& result, const complex_adaptor<Backend>& arg)
0916 {
0917    result = arg.imag_data();
0918 }
0919
0920 template <class Backend, class T>
0921 inline void eval_set_imag(complex_adaptor<Backend>& result, const T& arg)
0922 {
0923    result.imag_data() = arg;
0924 }
0925
0926 template <class Backend, class T>
0927 inline void eval_set_real(complex_adaptor<Backend>& result, const T& arg)
0928 {
0929    result.real_data() = arg;
0930 }
0931
0932 template <class Backend>
0933 inline std::size_t hash_value(const complex_adaptor<Backend>& val)
0934 {
0935    std::size_t result  = hash_value(val.real_data());
0936    std::size_t result2 = hash_value(val.imag_data());
0937    boost::multiprecision::detail::hash_combine(result, result2);
0938    return result;
0939 }
0940
0941 } // namespace backends
0942
0943 template <class Backend>
0944 struct number_category<complex_adaptor<Backend> > : public std::integral_constant<int, boost::multiprecision::number_kind_complex>
0945 {};
0946
0947 template <class Backend, expression_template_option ExpressionTemplates>
0948 struct component_type<number<complex_adaptor<Backend>, ExpressionTemplates> >
0949 {
0950    using type = number<Backend, ExpressionTemplates>;
0951 };
0952
0953 template <class Backend, expression_template_option ExpressionTemplates>
0954 struct complex_result_from_scalar<number<Backend, ExpressionTemplates> >
0955 {
0956    using type = number<complex_adaptor<Backend>, ExpressionTemplates>;
0957 };
0958
0959 namespace detail {
0960    template <class Backend>
0961    struct is_variable_precision<complex_adaptor<Backend> > : public is_variable_precision<Backend>
0962    {};
0963 #ifdef BOOST_HAS_INT128
0964    template <class Backend>
0965    struct is_convertible_arithmetic<int128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<int128_type, Backend>
0966    {};
0967    template <class Backend>
0968    struct is_convertible_arithmetic<uint128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<uint128_type, Backend>
0969    {};
0970 #endif
0971 #ifdef BOOST_HAS_FLOAT128
0972    template <class Backend>
0973    struct is_convertible_arithmetic<float128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<float128_type, Backend>
0974    {};
0975 #endif
0976    } // namespace detail
0977
0978
0979
0980 template <class Backend, expression_template_option ExpressionTemplates>
0981 struct complex_result_from_scalar<number<backends::debug_adaptor<Backend>, ExpressionTemplates> >
0982 {
0983    using type = number<backends::debug_adaptor<complex_adaptor<Backend> >, ExpressionTemplates>;
0984 };
0985
0986 template <class Backend, expression_template_option ExpressionTemplates>
0987 struct complex_result_from_scalar<number<backends::logged_adaptor<Backend>, ExpressionTemplates> >
0988 {
0989    using type = number<backends::logged_adaptor<complex_adaptor<Backend> >, ExpressionTemplates>;
0990 };
0991
0992 }
0993
0994 } // namespace boost::multiprecision
0995
0996 #endif