ceph/src/boost/boost/multiprecision/complex_adaptor.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //  Copyright 2018 John Maddock. Distributed under the Boost
   3 //  Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MULTIPRECISION_COMPLEX_ADAPTOR_HPP
   7 #define BOOST_MULTIPRECISION_COMPLEX_ADAPTOR_HPP
   8
   9 #include <boost/multiprecision/number.hpp>
  10 #include <cstdint>
  11 #include <boost/multiprecision/detail/digits.hpp>
  12 #include <boost/multiprecision/detail/hash.hpp>
  13 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
  14 #include <cmath>
  15 #include <algorithm>
  16 #include <complex>
  17
  18 namespace boost {
  19 namespace multiprecision {
  20 namespace backends {
  21
  22 template <class Backend>
  23 struct debug_adaptor;
  24
  25 template <class Backend>
  26 struct logged_adaptor;
  27
  28 template <class Backend>
  29 struct complex_adaptor
  30 {
  31  protected:
  32    Backend m_real, m_imag;
  33
  34  public:
  35    Backend& real_data()
  36    {
  37       return m_real;
  38    }
  39    const Backend& real_data() const
  40    {
  41       return m_real;
  42    }
  43    Backend& imag_data()
  44    {
  45       return m_imag;
  46    }
  47    const Backend& imag_data() const
  48    {
  49       return m_imag;
  50    }
  51
  52    using signed_types = typename Backend::signed_types  ;
  53    using unsigned_types = typename Backend::unsigned_types;
  54    using float_types = typename Backend::float_types   ;
  55    using exponent_type = typename Backend::exponent_type ;
  56
  57    complex_adaptor() {}
  58    complex_adaptor(const complex_adaptor& o) : m_real(o.real_data()), m_imag(o.imag_data()) {}
  59    // Rvalue construct:
  60    complex_adaptor(complex_adaptor&& o) : m_real(std::move(o.real_data())), m_imag(std::move(o.imag_data()))
  61    {}
  62    complex_adaptor(const Backend& val)
  63        : m_real(val)
  64    {}
  65
  66    template <class T>
  67    complex_adaptor(const T& val, const typename std::enable_if<std::is_convertible<T, Backend>::value>::type* = nullptr)
  68        : m_real(val)
  69    {}
  70
  71    complex_adaptor(const std::complex<float>& val)
  72    {
  73       m_real = (long double)val.real();
  74       m_imag = (long double)val.imag();
  75    }
  76    complex_adaptor(const std::complex<double>& val)
  77    {
  78       m_real = (long double)val.real();
  79       m_imag = (long double)val.imag();
  80    }
  81    complex_adaptor(const std::complex<long double>& val)
  82    {
  83       m_real = val.real();
  84       m_imag = val.imag();
  85    }
  86    template <class T, class U>
  87    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)
  88       : m_real(a), m_imag(b) {}
  89    template <class T, class U>
  90    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)
  91       : m_real(static_cast<T&&>(a)), m_imag(b) {}
  92    template <class T, class U>
  93    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)
  94       : m_real(static_cast<T&&>(a)), m_imag(static_cast<U&&>(b)) {}
  95    template <class T, class U>
  96    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)
  97       : m_real(a), m_imag(static_cast<U&&>(b)) {}
  98
  99    complex_adaptor& operator=(const complex_adaptor& o)
 100    {
 101       m_real = o.real_data();
 102       m_imag = o.imag_data();
 103       return *this;
 104    }
 105    // rvalue assign:
 106    complex_adaptor& operator=(complex_adaptor&& o) noexcept
 107    {
 108       m_real = std::move(o.real_data());
 109       m_imag = std::move(o.imag_data());
 110       return *this;
 111    }
 112    template <class V>
 113    typename std::enable_if<std::is_assignable<Backend, V>::value, complex_adaptor&>::type operator=(const V& v)
 114    {
 115       using ui_type = typename std::tuple_element<0, unsigned_types>::type;
 116       m_real = v;
 117       m_imag = ui_type(0u);
 118       return *this;
 119    }
 120    template <class T>
 121    complex_adaptor& operator=(const std::complex<T>& val)
 122    {
 123       m_real = (long double)val.real();
 124       m_imag = (long double)val.imag();
 125       return *this;
 126    }
 127    complex_adaptor& operator=(const char* s)
 128    {
 129       using ui_type = typename std::tuple_element<0, unsigned_types>::type;
 130       ui_type                                           zero = 0u;
 131
 132       using default_ops::eval_fpclassify;
 133
 134       if (s && (*s == '('))
 135       {
 136          std::string part;
 137          const char* p = ++s;
 138          while (*p && (*p != ',') && (*p != ')'))
 139             ++p;
 140          part.assign(s, p);
 141          if (part.size())
 142             real_data() = part.c_str();
 143          else
 144             real_data() = zero;
 145          s = p;
 146          if (*p && (*p != ')'))
 147          {
 148             ++p;
 149             while (*p && (*p != ')'))
 150                ++p;
 151             part.assign(s + 1, p);
 152          }
 153          else
 154             part.erase();
 155          if (part.size())
 156             imag_data() = part.c_str();
 157          else
 158             imag_data() = zero;
 159
 160          if (eval_fpclassify(imag_data()) == static_cast<int>(FP_NAN))
 161          {
 162             real_data() = imag_data();
 163          }
 164       }
 165       else
 166       {
 167          real_data() = s;
 168          imag_data() = zero;
 169       }
 170       return *this;
 171    }
 172
 173    int compare(const complex_adaptor& o) const
 174    {
 175       // They are either equal or not:
 176       return (m_real.compare(o.real_data()) == 0) && (m_imag.compare(o.imag_data()) == 0) ? 0 : 1;
 177    }
 178    template <class T>
 179    int compare(const T& val) const
 180    {
 181       using default_ops::eval_is_zero;
 182       return (m_real.compare(val) == 0) && eval_is_zero(m_imag) ? 0 : 1;
 183    }
 184    void swap(complex_adaptor& o)
 185    {
 186       real_data().swap(o.real_data());
 187       imag_data().swap(o.imag_data());
 188    }
 189    std::string str(std::streamsize dig, std::ios_base::fmtflags f) const
 190    {
 191       using default_ops::eval_is_zero;
 192       if (eval_is_zero(imag_data()))
 193          return m_real.str(dig, f);
 194       return "(" + m_real.str(dig, f) + "," + m_imag.str(dig, f) + ")";
 195    }
 196    void negate()
 197    {
 198       m_real.negate();
 199       m_imag.negate();
 200    }
 201 };
 202
 203 template <class Backend, class T>
 204 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
 205 {
 206    return a.compare(b) == 0;
 207 }
 208
 209 template <class Backend>
 210 inline void eval_add(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
 211 {
 212    eval_add(result.real_data(), o.real_data());
 213    eval_add(result.imag_data(), o.imag_data());
 214 }
 215 template <class Backend>
 216 inline void eval_subtract(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
 217 {
 218    eval_subtract(result.real_data(), o.real_data());
 219    eval_subtract(result.imag_data(), o.imag_data());
 220 }
 221 template <class Backend>
 222 inline void eval_multiply(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
 223 {
 224    Backend t1, t2, t3;
 225    eval_multiply(t1, result.real_data(), o.real_data());
 226    eval_multiply(t2, result.imag_data(), o.imag_data());
 227    eval_subtract(t3, t1, t2);
 228    eval_multiply(t1, result.real_data(), o.imag_data());
 229    eval_multiply(t2, result.imag_data(), o.real_data());
 230    eval_add(t1, t2);
 231    result.real_data() = std::move(t3);
 232    result.imag_data() = std::move(t1);
 233 }
 234 template <class Backend>
 235 inline void eval_divide(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& z)
 236 {
 237    // (a+bi) / (c + di)
 238    using default_ops::eval_add;
 239    using default_ops::eval_divide;
 240    using default_ops::eval_fabs;
 241    using default_ops::eval_is_zero;
 242    using default_ops::eval_multiply;
 243    using default_ops::eval_subtract;
 244    Backend t1, t2;
 245
 246    //
 247    // Backup sign bits for later, so we can fix up
 248    // signed zeros at the end:
 249    //
 250    int a_sign = eval_signbit(result.real_data());
 251    int b_sign = eval_signbit(result.imag_data());
 252    int c_sign = eval_signbit(z.real_data());
 253    int d_sign = eval_signbit(z.imag_data());
 254
 255    if (eval_is_zero(z.imag_data()))
 256    {
 257       eval_divide(result.real_data(), z.real_data());
 258       eval_divide(result.imag_data(), z.real_data());
 259    }
 260    else
 261    {
 262       eval_fabs(t1, z.real_data());
 263       eval_fabs(t2, z.imag_data());
 264       if (t1.compare(t2) < 0)
 265       {
 266          eval_divide(t1, z.real_data(), z.imag_data()); // t1 = c/d
 267          eval_multiply(t2, z.real_data(), t1);
 268          eval_add(t2, z.imag_data()); // denom = c * (c/d) + d
 269          Backend t_real(result.real_data());
 270          // real = (a * (c/d) + b) / (denom)
 271          eval_multiply(result.real_data(), t1);
 272          eval_add(result.real_data(), result.imag_data());
 273          eval_divide(result.real_data(), t2);
 274          // imag = (b * c/d - a) / denom
 275          eval_multiply(result.imag_data(), t1);
 276          eval_subtract(result.imag_data(), t_real);
 277          eval_divide(result.imag_data(), t2);
 278       }
 279       else
 280       {
 281          eval_divide(t1, z.imag_data(), z.real_data()); // t1 = d/c
 282          eval_multiply(t2, z.imag_data(), t1);
 283          eval_add(t2, z.real_data()); // denom = d * d/c + c
 284
 285          Backend r_t(result.real_data());
 286          Backend i_t(result.imag_data());
 287
 288          // real = (b * d/c + a) / denom
 289          eval_multiply(result.real_data(), result.imag_data(), t1);
 290          eval_add(result.real_data(), r_t);
 291          eval_divide(result.real_data(), t2);
 292          // imag = (-a * d/c + b) / denom
 293          eval_multiply(result.imag_data(), r_t, t1);
 294          result.imag_data().negate();
 295          eval_add(result.imag_data(), i_t);
 296          eval_divide(result.imag_data(), t2);
 297       }
 298    }
 299    //
 300    // Finish off by fixing up signed zeros.
 301    //
 302    // This sets the signs "as if" we had evaluated the result using:
 303    //
 304    // real = (ac + bd) / (c^2 + d^2)
 305    // imag = (bc - ad) / (c^2 + d^2)
 306    //
 307    // ie a zero is negative only if the two parts of the numerator
 308    // are both negative and zero.
 309    //
 310    if (eval_is_zero(result.real_data()))
 311    {
 312       int r_sign = eval_signbit(result.real_data());
 313       int r_required = (a_sign != c_sign) && (b_sign != d_sign);
 314       if (r_required != r_sign)
 315          result.real_data().negate();
 316    }
 317    if (eval_is_zero(result.imag_data()))
 318    {
 319       int i_sign = eval_signbit(result.imag_data());
 320       int i_required = (b_sign != c_sign) && (a_sign == d_sign);
 321       if (i_required != i_sign)
 322          result.imag_data().negate();
 323    }
 324 }
 325 template <class Backend, class T>
 326 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_add(complex_adaptor<Backend>& result, const T& scalar)
 327 {
 328    using default_ops::eval_add;
 329    eval_add(result.real_data(), scalar);
 330 }
 331 template <class Backend, class T>
 332 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_subtract(complex_adaptor<Backend>& result, const T& scalar)
 333 {
 334    using default_ops::eval_subtract;
 335    eval_subtract(result.real_data(), scalar);
 336 }
 337 template <class Backend, class T>
 338 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_multiply(complex_adaptor<Backend>& result, const T& scalar)
 339 {
 340    using default_ops::eval_multiply;
 341    eval_multiply(result.real_data(), scalar);
 342    eval_multiply(result.imag_data(), scalar);
 343 }
 344 template <class Backend, class T>
 345 inline typename std::enable_if< !std::is_same<complex_adaptor<Backend>, T>::value>::type eval_divide(complex_adaptor<Backend>& result, const T& scalar)
 346 {
 347    using default_ops::eval_divide;
 348    eval_divide(result.real_data(), scalar);
 349    eval_divide(result.imag_data(), scalar);
 350 }
 351 // Optimised 3 arg versions:
 352 template <class Backend, class T>
 353 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)
 354 {
 355    using default_ops::eval_add;
 356    eval_add(result.real_data(), a.real_data(), scalar);
 357    result.imag_data() = a.imag_data();
 358 }
 359 template <class Backend, class T>
 360 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)
 361 {
 362    using default_ops::eval_subtract;
 363    eval_subtract(result.real_data(), a.real_data(), scalar);
 364    result.imag_data() = a.imag_data();
 365 }
 366 template <class Backend, class T>
 367 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)
 368 {
 369    using default_ops::eval_multiply;
 370    eval_multiply(result.real_data(), a.real_data(), scalar);
 371    eval_multiply(result.imag_data(), a.imag_data(), scalar);
 372 }
 373 template <class Backend, class T>
 374 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)
 375 {
 376    using default_ops::eval_divide;
 377    eval_divide(result.real_data(), a.real_data(), scalar);
 378    eval_divide(result.imag_data(), a.imag_data(), scalar);
 379 }
 380
 381 template <class Backend>
 382 inline bool eval_is_zero(const complex_adaptor<Backend>& val) noexcept
 383 {
 384    using default_ops::eval_is_zero;
 385    return eval_is_zero(val.real_data()) && eval_is_zero(val.imag_data());
 386 }
 387 template <class Backend>
 388 inline int eval_get_sign(const complex_adaptor<Backend>&)
 389 {
 390    static_assert(sizeof(Backend) == UINT_MAX, "Complex numbers have no sign bit."); // designed to always fail
 391    return 0;
 392 }
 393
 394 template <class Result, class Backend>
 395 inline typename std::enable_if< !boost::multiprecision::detail::is_complex<Result>::value>::type eval_convert_to(Result* result, const complex_adaptor<Backend>& val)
 396 {
 397    using default_ops::eval_convert_to;
 398    using default_ops::eval_is_zero;
 399    if (!eval_is_zero(val.imag_data()))
 400    {
 401       BOOST_MP_THROW_EXCEPTION(std::runtime_error("Could not convert imaginary number to scalar."));
 402    }
 403    eval_convert_to(result, val.real_data());
 404 }
 405
 406 template <class Backend, class T>
 407 inline void assign_components(complex_adaptor<Backend>& result, const T& a, const T& b)
 408 {
 409    result.real_data() = a;
 410    result.imag_data() = b;
 411 }
 412
 413 //
 414 // Native non-member operations:
 415 //
 416 template <class Backend>
 417 inline void eval_sqrt(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& val)
 418 {
 419    // Use the following:
 420    // sqrt(z) = (s, zi / 2s)       for zr >= 0
 421    //           (|zi| / 2s, +-s)   for zr <  0
 422    // where s = sqrt{ [ |zr| + sqrt(zr^2 + zi^2) ] / 2 },
 423    // and the +- sign is the same as the sign of zi.
 424    using default_ops::eval_abs;
 425    using default_ops::eval_add;
 426    using default_ops::eval_divide;
 427    using default_ops::eval_get_sign;
 428    using default_ops::eval_is_zero;
 429
 430    if (eval_is_zero(val.imag_data()) && (eval_get_sign(val.real_data()) >= 0))
 431    {
 432       constexpr const typename std::tuple_element<0, typename Backend::unsigned_types>::type zero = 0u;
 433       eval_sqrt(result.real_data(), val.real_data());
 434       result.imag_data() = zero;
 435       return;
 436    }
 437
 438    const bool __my_real_part_is_neg(eval_get_sign(val.real_data()) < 0);
 439
 440    Backend __my_real_part_fabs(val.real_data());
 441    if (__my_real_part_is_neg)
 442       __my_real_part_fabs.negate();
 443
 444    Backend t, __my_sqrt_part;
 445    eval_abs(__my_sqrt_part, val);
 446    eval_add(__my_sqrt_part, __my_real_part_fabs);
 447    eval_ldexp(t, __my_sqrt_part, -1);
 448    eval_sqrt(__my_sqrt_part, t);
 449
 450    if (__my_real_part_is_neg == false)
 451    {
 452       eval_ldexp(t, __my_sqrt_part, 1);
 453       eval_divide(result.imag_data(), val.imag_data(), t);
 454       result.real_data() = __my_sqrt_part;
 455    }
 456    else
 457    {
 458       const bool __my_imag_part_is_neg(eval_get_sign(val.imag_data()) < 0);
 459
 460       Backend __my_imag_part_fabs(val.imag_data());
 461       if (__my_imag_part_is_neg)
 462          __my_imag_part_fabs.negate();
 463
 464       eval_ldexp(t, __my_sqrt_part, 1);
 465       eval_divide(result.real_data(), __my_imag_part_fabs, t);
 466       if (__my_imag_part_is_neg)
 467          __my_sqrt_part.negate();
 468       result.imag_data() = __my_sqrt_part;
 469    }
 470 }
 471
 472 template <class Backend>
 473 inline void eval_abs(Backend& result, const complex_adaptor<Backend>& val)
 474 {
 475    Backend t1, t2;
 476    eval_multiply(t1, val.real_data(), val.real_data());
 477    eval_multiply(t2, val.imag_data(), val.imag_data());
 478    eval_add(t1, t2);
 479    eval_sqrt(result, t1);
 480 }
 481
 482 template <class Backend>
 483 inline void eval_pow(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& b, const complex_adaptor<Backend>& e)
 484 {
 485    using default_ops::eval_acos;
 486    using default_ops::eval_cos;
 487    using default_ops::eval_exp;
 488    using default_ops::eval_get_sign;
 489    using default_ops::eval_is_zero;
 490    using default_ops::eval_multiply;
 491    using default_ops::eval_sin;
 492
 493    if (eval_is_zero(e))
 494    {
 495       typename std::tuple_element<0, typename Backend::unsigned_types>::type one(1);
 496       result = one;
 497       return;
 498    }
 499    else if (eval_is_zero(b))
 500    {
 501       if (eval_is_zero(e.real_data()))
 502       {
 503          Backend n          = std::numeric_limits<number<Backend> >::quiet_NaN().backend();
 504          result.real_data() = n;
 505          result.imag_data() = n;
 506       }
 507       else if (eval_get_sign(e.real_data()) < 0)
 508       {
 509          Backend n          = std::numeric_limits<number<Backend> >::infinity().backend();
 510          result.real_data() = n;
 511          typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
 512          if (eval_is_zero(e.imag_data()))
 513             result.imag_data() = zero;
 514          else
 515             result.imag_data() = n;
 516       }
 517       else
 518       {
 519          typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
 520          result = zero;
 521       }
 522       return;
 523    }
 524    complex_adaptor<Backend> t;
 525    eval_log(t, b);
 526    eval_multiply(t, e);
 527    eval_exp(result, t);
 528 }
 529
 530 template <class Backend>
 531 inline void eval_exp(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 532 {
 533    using default_ops::eval_cos;
 534    using default_ops::eval_exp;
 535    using default_ops::eval_is_zero;
 536    using default_ops::eval_multiply;
 537    using default_ops::eval_sin;
 538
 539    if (eval_is_zero(arg.imag_data()))
 540    {
 541       eval_exp(result.real_data(), arg.real_data());
 542       typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
 543       result.imag_data() = zero;
 544       return;
 545    }
 546    eval_cos(result.real_data(), arg.imag_data());
 547    eval_sin(result.imag_data(), arg.imag_data());
 548    Backend e;
 549    eval_exp(e, arg.real_data());
 550    if (eval_is_zero(result.real_data()))
 551       eval_multiply(result.imag_data(), e);
 552    else if (eval_is_zero(result.imag_data()))
 553       eval_multiply(result.real_data(), e);
 554    else
 555       eval_multiply(result, e);
 556 }
 557
 558 template <class Backend>
 559 inline void eval_log(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 560 {
 561    using default_ops::eval_add;
 562    using default_ops::eval_atan2;
 563    using default_ops::eval_get_sign;
 564    using default_ops::eval_is_zero;
 565    using default_ops::eval_log;
 566    using default_ops::eval_multiply;
 567
 568    if (eval_is_zero(arg.imag_data()) && (eval_get_sign(arg.real_data()) >= 0))
 569    {
 570       eval_log(result.real_data(), arg.real_data());
 571       typename std::tuple_element<0, typename Backend::unsigned_types>::type zero(0);
 572       result.imag_data() = zero;
 573       return;
 574    }
 575
 576    Backend t1, t2;
 577    eval_multiply(t1, arg.real_data(), arg.real_data());
 578    eval_multiply(t2, arg.imag_data(), arg.imag_data());
 579    eval_add(t1, t2);
 580    eval_log(t2, t1);
 581    eval_ldexp(result.real_data(), t2, -1);
 582    eval_atan2(result.imag_data(), arg.imag_data(), arg.real_data());
 583 }
 584
 585 template <class Backend>
 586 inline void eval_log10(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 587 {
 588    using default_ops::eval_divide;
 589    using default_ops::eval_log;
 590
 591    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
 592
 593    Backend ten;
 594    ten = ui_type(10);
 595    Backend l_ten;
 596    eval_log(l_ten, ten);
 597    eval_log(result, arg);
 598    eval_divide(result, l_ten);
 599 }
 600
 601 template <class Backend>
 602 inline void eval_sin(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 603 {
 604    using default_ops::eval_cos;
 605    using default_ops::eval_cosh;
 606    using default_ops::eval_sin;
 607    using default_ops::eval_sinh;
 608
 609    Backend t1, t2, t3;
 610    eval_sin(t1, arg.real_data());
 611    eval_cosh(t2, arg.imag_data());
 612    eval_multiply(t3, t1, t2);
 613
 614    eval_cos(t1, arg.real_data());
 615    eval_sinh(t2, arg.imag_data());
 616    eval_multiply(result.imag_data(), t1, t2);
 617    result.real_data() = t3;
 618 }
 619
 620 template <class Backend>
 621 inline void eval_cos(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 622 {
 623    using default_ops::eval_cos;
 624    using default_ops::eval_cosh;
 625    using default_ops::eval_sin;
 626    using default_ops::eval_sinh;
 627
 628    Backend t1, t2, t3;
 629    eval_cos(t1, arg.real_data());
 630    eval_cosh(t2, arg.imag_data());
 631    eval_multiply(t3, t1, t2);
 632
 633    eval_sin(t1, arg.real_data());
 634    eval_sinh(t2, arg.imag_data());
 635    eval_multiply(result.imag_data(), t1, t2);
 636    result.imag_data().negate();
 637    result.real_data() = t3;
 638 }
 639
 640 template <class T>
 641 void tanh_imp(const T& r, const T& i, T& r_result, T& i_result)
 642 {
 643    using default_ops::eval_tan;
 644    using default_ops::eval_sinh;
 645    using default_ops::eval_add;
 646    using default_ops::eval_fpclassify;
 647    using default_ops::eval_get_sign;
 648
 649    using ui_type = typename std::tuple_element<0, typename T::unsigned_types>::type;
 650    ui_type one(1);
 651    //
 652    // Set:
 653    // t = tan(i);
 654    // s = sinh(r);
 655    // b = s * (1 + t^2);
 656    // d = 1 + b * s;
 657    //
 658    T t, s, b, d;
 659    eval_tan(t, i);
 660    eval_sinh(s, r);
 661    eval_multiply(d, t, t);
 662    eval_add(d, one);
 663    eval_multiply(b, d, s);
 664    eval_multiply(d, b, s);
 665    eval_add(d, one);
 666
 667    if (eval_fpclassify(d) == FP_INFINITE)
 668    {
 669       r_result = one;
 670       if (eval_get_sign(s) < 0)
 671          r_result.negate();
 672       //
 673       // Imaginary part is a signed zero:
 674       //
 675       ui_type zero(0);
 676       i_result = zero;
 677       if (eval_get_sign(t) < 0)
 678          i_result.negate();
 679    }
 680    //
 681    // Real part is sqrt(1 + s^2) * b / d;
 682    // Imaginary part is t / d;
 683    //
 684    eval_divide(i_result, t, d);
 685    //
 686    // variable t is now spare, as is r_result.
 687    //
 688    eval_multiply(t, s, s);
 689    eval_add(t, one);
 690    eval_sqrt(r_result, t);
 691    eval_multiply(t, r_result, b);
 692    eval_divide(r_result, t, d);
 693 }
 694
 695 template <class Backend>
 696 inline void eval_tanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 697 {
 698    tanh_imp(arg.real_data(), arg.imag_data(), result.real_data(), result.imag_data());
 699 }
 700 template <class Backend>
 701 inline void eval_tan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 702 {
 703    Backend t(arg.imag_data());
 704    t.negate();
 705    tanh_imp(t, arg.real_data(), result.imag_data(), result.real_data());
 706    result.imag_data().negate();
 707 }
 708
 709 template <class Backend>
 710 inline void eval_asin(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 711 {
 712    using default_ops::eval_add;
 713    using default_ops::eval_multiply;
 714
 715    if (eval_is_zero(arg))
 716    {
 717       result = arg;
 718       return;
 719    }
 720
 721    complex_adaptor<Backend> t1, t2;
 722    assign_components(t1, arg.imag_data(), arg.real_data());
 723    t1.real_data().negate();
 724    eval_asinh(t2, t1);
 725
 726    assign_components(result, t2.imag_data(), t2.real_data());
 727    result.imag_data().negate();
 728 }
 729
 730 template <class Backend>
 731 inline void eval_acos(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 732 {
 733    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
 734
 735    using default_ops::eval_asin;
 736
 737    Backend half_pi, t1;
 738    t1 = static_cast<ui_type>(1u);
 739    eval_asin(half_pi, t1);
 740    eval_asin(result, arg);
 741    result.negate();
 742    eval_add(result.real_data(), half_pi);
 743 }
 744
 745 template <class Backend>
 746 inline void eval_atan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 747 {
 748    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
 749    ui_type                                                             one = (ui_type)1u;
 750
 751    using default_ops::eval_add;
 752    using default_ops::eval_is_zero;
 753    using default_ops::eval_log;
 754    using default_ops::eval_subtract;
 755
 756    complex_adaptor<Backend> __my_z_times_i, t1, t2, t3;
 757    assign_components(__my_z_times_i, arg.imag_data(), arg.real_data());
 758    __my_z_times_i.real_data().negate();
 759
 760    eval_add(t1, __my_z_times_i, one);
 761    eval_log(t2, t1);
 762    eval_subtract(t1, one, __my_z_times_i);
 763    eval_log(t3, t1);
 764    eval_subtract(t1, t3, t2);
 765
 766    eval_ldexp(result.real_data(), t1.imag_data(), -1);
 767    eval_ldexp(result.imag_data(), t1.real_data(), -1);
 768    if (!eval_is_zero(result.real_data()))
 769       result.real_data().negate();
 770 }
 771
 772 template <class Backend>
 773 inline void eval_sinh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 774 {
 775    using default_ops::eval_cos;
 776    using default_ops::eval_cosh;
 777    using default_ops::eval_sin;
 778    using default_ops::eval_sinh;
 779
 780    Backend t1, t2, t3;
 781    eval_cos(t1, arg.imag_data());
 782    eval_sinh(t2, arg.real_data());
 783    eval_multiply(t3, t1, t2);
 784
 785    eval_cosh(t1, arg.real_data());
 786    eval_sin(t2, arg.imag_data());
 787    eval_multiply(result.imag_data(), t1, t2);
 788    result.real_data() = t3;
 789 }
 790
 791 template <class Backend>
 792 inline void eval_cosh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 793 {
 794    using default_ops::eval_cos;
 795    using default_ops::eval_cosh;
 796    using default_ops::eval_sin;
 797    using default_ops::eval_sinh;
 798
 799    Backend t1, t2, t3;
 800    eval_cos(t1, arg.imag_data());
 801    eval_cosh(t2, arg.real_data());
 802    eval_multiply(t3, t1, t2);
 803
 804    eval_sin(t1, arg.imag_data());
 805    eval_sinh(t2, arg.real_data());
 806    eval_multiply(result.imag_data(), t1, t2);
 807    result.real_data() = t3;
 808 }
 809
 810 template <class Backend>
 811 inline void eval_asinh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 812 {
 813    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
 814    ui_type                                                             one = (ui_type)1u;
 815
 816    using default_ops::eval_add;
 817    using default_ops::eval_log;
 818    using default_ops::eval_multiply;
 819
 820    complex_adaptor<Backend> t1, t2;
 821    eval_multiply(t1, arg, arg);
 822    eval_add(t1, one);
 823    eval_sqrt(t2, t1);
 824    eval_add(t2, arg);
 825    eval_log(result, t2);
 826 }
 827
 828 template <class Backend>
 829 inline void eval_acosh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 830 {
 831    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
 832    ui_type                                                             one = (ui_type)1u;
 833
 834    using default_ops::eval_add;
 835    using default_ops::eval_divide;
 836    using default_ops::eval_log;
 837    using default_ops::eval_multiply;
 838    using default_ops::eval_subtract;
 839
 840    complex_adaptor<Backend> __my_zp(arg);
 841    eval_add(__my_zp.real_data(), one);
 842    complex_adaptor<Backend> __my_zm(arg);
 843    eval_subtract(__my_zm.real_data(), one);
 844
 845    complex_adaptor<Backend> t1, t2;
 846    eval_divide(t1, __my_zm, __my_zp);
 847    eval_sqrt(t2, t1);
 848    eval_multiply(t2, __my_zp);
 849    eval_add(t2, arg);
 850    eval_log(result, t2);
 851 }
 852
 853 template <class Backend>
 854 inline void eval_atanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 855 {
 856    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
 857    ui_type                                                             one = (ui_type)1u;
 858
 859    using default_ops::eval_add;
 860    using default_ops::eval_divide;
 861    using default_ops::eval_log;
 862    using default_ops::eval_multiply;
 863    using default_ops::eval_subtract;
 864
 865    complex_adaptor<Backend> t1, t2, t3;
 866    eval_add(t1, arg, one);
 867    eval_log(t2, t1);
 868    eval_subtract(t1, one, arg);
 869    eval_log(t3, t1);
 870    eval_subtract(t2, t3);
 871
 872    eval_ldexp(result.real_data(), t2.real_data(), -1);
 873    eval_ldexp(result.imag_data(), t2.imag_data(), -1);
 874 }
 875
 876 template <class Backend>
 877 inline void eval_conj(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 878 {
 879    result = arg;
 880    result.imag_data().negate();
 881 }
 882
 883 template <class Backend>
 884 inline void eval_proj(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 885 {
 886    using default_ops::eval_get_sign;
 887
 888    using ui_type = typename std::tuple_element<0, typename Backend::unsigned_types>::type;
 889    ui_type                                                             zero = (ui_type)0u;
 890
 891    int c1 = eval_fpclassify(arg.real_data());
 892    int c2 = eval_fpclassify(arg.imag_data());
 893    if (c1 == FP_INFINITE)
 894    {
 895       result.real_data() = arg.real_data();
 896       if (eval_get_sign(result.real_data()) < 0)
 897          result.real_data().negate();
 898       result.imag_data() = zero;
 899       if (eval_get_sign(arg.imag_data()) < 0)
 900          result.imag_data().negate();
 901    }
 902    else if (c2 == FP_INFINITE)
 903    {
 904       result.real_data() = arg.imag_data();
 905       if (eval_get_sign(result.real_data()) < 0)
 906          result.real_data().negate();
 907       result.imag_data() = zero;
 908       if (eval_get_sign(arg.imag_data()) < 0)
 909          result.imag_data().negate();
 910    }
 911    else
 912       result = arg;
 913 }
 914
 915 template <class Backend>
 916 inline void eval_real(Backend& result, const complex_adaptor<Backend>& arg)
 917 {
 918    result = arg.real_data();
 919 }
 920 template <class Backend>
 921 inline void eval_imag(Backend& result, const complex_adaptor<Backend>& arg)
 922 {
 923    result = arg.imag_data();
 924 }
 925
 926 template <class Backend, class T>
 927 inline void eval_set_imag(complex_adaptor<Backend>& result, const T& arg)
 928 {
 929    result.imag_data() = arg;
 930 }
 931
 932 template <class Backend, class T>
 933 inline void eval_set_real(complex_adaptor<Backend>& result, const T& arg)
 934 {
 935    result.real_data() = arg;
 936 }
 937
 938 template <class Backend>
 939 inline std::size_t hash_value(const complex_adaptor<Backend>& val)
 940 {
 941    std::size_t result  = hash_value(val.real_data());
 942    std::size_t result2 = hash_value(val.imag_data());
 943    boost::multiprecision::detail::hash_combine(result, result2);
 944    return result;
 945 }
 946
 947 } // namespace backends
 948
 949 using boost::multiprecision::backends::complex_adaptor;
 950
 951 template <class Backend>
 952 struct number_category<complex_adaptor<Backend> > : public std::integral_constant<int, boost::multiprecision::number_kind_complex>
 953 {};
 954
 955 template <class Backend, expression_template_option ExpressionTemplates>
 956 struct component_type<number<complex_adaptor<Backend>, ExpressionTemplates> >
 957 {
 958    using type = number<Backend, ExpressionTemplates>;
 959 };
 960
 961 template <class Backend, expression_template_option ExpressionTemplates>
 962 struct complex_result_from_scalar<number<Backend, ExpressionTemplates> >
 963 {
 964    using type = number<complex_adaptor<Backend>, ExpressionTemplates>;
 965 };
 966
 967 namespace detail {
 968    template <class Backend>
 969    struct is_variable_precision<complex_adaptor<Backend> > : public is_variable_precision<Backend>
 970    {};
 971 #ifdef BOOST_HAS_INT128
 972    template <class Backend>
 973    struct is_convertible_arithmetic<int128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<int128_type, Backend>
 974    {};
 975    template <class Backend>
 976    struct is_convertible_arithmetic<uint128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<uint128_type, Backend>
 977    {};
 978 #endif
 979 #ifdef BOOST_HAS_FLOAT128
 980    template <class Backend>
 981    struct is_convertible_arithmetic<float128_type, complex_adaptor<Backend> > : is_convertible_arithmetic<float128_type, Backend>
 982    {};
 983 #endif
 984    } // namespace detail
 985
 986
 987
 988 template <class Backend, expression_template_option ExpressionTemplates>
 989 struct complex_result_from_scalar<number<backends::debug_adaptor<Backend>, ExpressionTemplates> >
 990 {
 991    using type = number<backends::debug_adaptor<complex_adaptor<Backend> >, ExpressionTemplates>;
 992 };
 993
 994 template <class Backend, expression_template_option ExpressionTemplates>
 995 struct complex_result_from_scalar<number<backends::logged_adaptor<Backend>, ExpressionTemplates> >
 996 {
 997    using type = number<backends::logged_adaptor<complex_adaptor<Backend> >, ExpressionTemplates>;
 998 };
 999
1000 }
1001
1002 } // namespace boost::multiprecision
1003
1004 #endif