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 <boost/cstdint.hpp>
  11 #include <boost/multiprecision/detail/digits.hpp>
  12 #include <boost/functional/hash_fwd.hpp>
  13 #include <boost/type_traits/is_complex.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 complex_adaptor
  24 {
  25  protected:
  26    Backend m_real, m_imag;
  27
  28  public:
  29    Backend& real_data()
  30    {
  31       return m_real;
  32    }
  33    const Backend& real_data() const
  34    {
  35       return m_real;
  36    }
  37    Backend& imag_data()
  38    {
  39       return m_imag;
  40    }
  41    const Backend& imag_data() const
  42    {
  43       return m_imag;
  44    }
  45
  46    typedef typename Backend::signed_types   signed_types;
  47    typedef typename Backend::unsigned_types unsigned_types;
  48    typedef typename Backend::float_types    float_types;
  49    typedef typename Backend::exponent_type  exponent_type;
  50
  51    complex_adaptor() {}
  52    complex_adaptor(const complex_adaptor& o) : m_real(o.real_data()), m_imag(o.imag_data()) {}
  53 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
  54    complex_adaptor(complex_adaptor&& o) : m_real(std::move(o.real_data())), m_imag(std::move(o.imag_data()))
  55    {}
  56 #endif
  57    complex_adaptor(const Backend& val)
  58        : m_real(val)
  59    {}
  60
  61    complex_adaptor(const std::complex<float>& val)
  62    {
  63       m_real = (long double)val.real();
  64       m_imag = (long double)val.imag();
  65    }
  66    complex_adaptor(const std::complex<double>& val)
  67    {
  68       m_real = (long double)val.real();
  69       m_imag = (long double)val.imag();
  70    }
  71    complex_adaptor(const std::complex<long double>& val)
  72    {
  73       m_real = val.real();
  74       m_imag = val.imag();
  75    }
  76
  77    complex_adaptor& operator=(const complex_adaptor& o)
  78    {
  79       m_real = o.real_data();
  80       m_imag = o.imag_data();
  81       return *this;
  82    }
  83 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
  84    complex_adaptor& operator=(complex_adaptor&& o) BOOST_NOEXCEPT
  85    {
  86       m_real = std::move(o.real_data());
  87       m_imag = std::move(o.imag_data());
  88       return *this;
  89    }
  90 #endif
  91    template <class V>
  92    complex_adaptor& operator=(const V& v)
  93    {
  94       typedef typename mpl::front<unsigned_types>::type ui_type;
  95       m_real = v;
  96       m_imag = ui_type(0u);
  97       return *this;
  98    }
  99    template <class T>
 100    complex_adaptor& operator=(const std::complex<T>& val)
 101    {
 102       m_real = (long double)val.real();
 103       m_imag = (long double)val.imag();
 104       return *this;
 105    }
 106    complex_adaptor& operator=(const char* s)
 107    {
 108       typedef typename mpl::front<unsigned_types>::type ui_type;
 109       ui_type                                           zero = 0u;
 110
 111       using default_ops::eval_fpclassify;
 112
 113       if (s && (*s == '('))
 114       {
 115          std::string part;
 116          const char* p = ++s;
 117          while (*p && (*p != ',') && (*p != ')'))
 118             ++p;
 119          part.assign(s, p);
 120          if (part.size())
 121             real_data() = part.c_str();
 122          else
 123             real_data() = zero;
 124          s = p;
 125          if (*p && (*p != ')'))
 126          {
 127             ++p;
 128             while (*p && (*p != ')'))
 129                ++p;
 130             part.assign(s + 1, p);
 131          }
 132          else
 133             part.erase();
 134          if (part.size())
 135             imag_data() = part.c_str();
 136          else
 137             imag_data() = zero;
 138
 139          if (eval_fpclassify(imag_data()) == (int)FP_NAN)
 140          {
 141             real_data() = imag_data();
 142          }
 143       }
 144       else
 145       {
 146          real_data() = s;
 147          imag_data() = zero;
 148       }
 149       return *this;
 150    }
 151
 152    int compare(const complex_adaptor& o) const
 153    {
 154       // They are either equal or not:
 155       return (m_real.compare(o.real_data()) == 0) && (m_imag.compare(o.imag_data()) == 0) ? 0 : 1;
 156    }
 157    template <class T>
 158    int compare(const T& val) const
 159    {
 160       using default_ops::eval_is_zero;
 161       return (m_real.compare(val) == 0) && eval_is_zero(m_imag) ? 0 : 1;
 162    }
 163    void swap(complex_adaptor& o)
 164    {
 165       real_data().swap(o.real_data());
 166       imag_data().swap(o.imag_data());
 167    }
 168    std::string str(std::streamsize dig, std::ios_base::fmtflags f) const
 169    {
 170       using default_ops::eval_is_zero;
 171       if (eval_is_zero(imag_data()))
 172          return m_real.str(dig, f);
 173       return "(" + m_real.str(dig, f) + "," + m_imag.str(dig, f) + ")";
 174    }
 175    void negate()
 176    {
 177       m_real.negate();
 178       m_imag.negate();
 179    }
 180 };
 181
 182 template <class Backend, class T>
 183 inline typename enable_if<is_arithmetic<T>, bool>::type eval_eq(const complex_adaptor<Backend>& a, const T& b) BOOST_NOEXCEPT
 184 {
 185    return a.compare(b) == 0;
 186 }
 187
 188 template <class Backend>
 189 inline void eval_add(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
 190 {
 191    eval_add(result.real_data(), o.real_data());
 192    eval_add(result.imag_data(), o.imag_data());
 193 }
 194 template <class Backend>
 195 inline void eval_subtract(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
 196 {
 197    eval_subtract(result.real_data(), o.real_data());
 198    eval_subtract(result.imag_data(), o.imag_data());
 199 }
 200 template <class Backend>
 201 inline void eval_multiply(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& o)
 202 {
 203    Backend t1, t2, t3;
 204    eval_multiply(t1, result.real_data(), o.real_data());
 205    eval_multiply(t2, result.imag_data(), o.imag_data());
 206    eval_subtract(t3, t1, t2);
 207    eval_multiply(t1, result.real_data(), o.imag_data());
 208    eval_multiply(t2, result.imag_data(), o.real_data());
 209    eval_add(t1, t2);
 210    result.real_data() = BOOST_MP_MOVE(t3);
 211    result.imag_data() = BOOST_MP_MOVE(t1);
 212 }
 213 template <class Backend>
 214 inline void eval_divide(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& z)
 215 {
 216    // (a+bi) / (c + di)
 217    using default_ops::eval_add;
 218    using default_ops::eval_divide;
 219    using default_ops::eval_fabs;
 220    using default_ops::eval_is_zero;
 221    using default_ops::eval_multiply;
 222    using default_ops::eval_subtract;
 223    Backend t1, t2;
 224
 225    if (eval_is_zero(z.imag_data()))
 226    {
 227       eval_divide(result.real_data(), z.real_data());
 228       eval_divide(result.imag_data(), z.real_data());
 229       return;
 230    }
 231
 232    eval_fabs(t1, z.real_data());
 233    eval_fabs(t2, z.imag_data());
 234    if (t1.compare(t2) < 0)
 235    {
 236       eval_divide(t1, z.real_data(), z.imag_data()); // t1 = c/d
 237       eval_multiply(t2, z.real_data(), t1);
 238       eval_add(t2, z.imag_data()); // denom = c * (c/d) + d
 239       Backend t_real(result.real_data());
 240       // real = (a * (c/d) + b) / (denom)
 241       eval_multiply(result.real_data(), t1);
 242       eval_add(result.real_data(), result.imag_data());
 243       eval_divide(result.real_data(), t2);
 244       // imag = (b * c/d - a) / denom
 245       eval_multiply(result.imag_data(), t1);
 246       eval_subtract(result.imag_data(), t_real);
 247       eval_divide(result.imag_data(), t2);
 248    }
 249    else
 250    {
 251       eval_divide(t1, z.imag_data(), z.real_data()); // t1 = d/c
 252       eval_multiply(t2, z.imag_data(), t1);
 253       eval_add(t2, z.real_data()); // denom = d * d/c + c
 254
 255       Backend r_t(result.real_data());
 256       Backend i_t(result.imag_data());
 257
 258       // real = (b * d/c + a) / denom
 259       eval_multiply(result.real_data(), result.imag_data(), t1);
 260       eval_add(result.real_data(), r_t);
 261       eval_divide(result.real_data(), t2);
 262       // imag = (-a * d/c + b) / denom
 263       eval_multiply(result.imag_data(), r_t, t1);
 264       result.imag_data().negate();
 265       eval_add(result.imag_data(), i_t);
 266       eval_divide(result.imag_data(), t2);
 267    }
 268 }
 269 template <class Backend, class T>
 270 inline typename boost::disable_if_c<boost::is_same<complex_adaptor<Backend>, T>::value>::type eval_add(complex_adaptor<Backend>& result, const T& scalar)
 271 {
 272    using default_ops::eval_add;
 273    eval_add(result.real_data(), scalar);
 274 }
 275 template <class Backend, class T>
 276 inline typename boost::disable_if_c<boost::is_same<complex_adaptor<Backend>, T>::value>::type eval_subtract(complex_adaptor<Backend>& result, const T& scalar)
 277 {
 278    using default_ops::eval_subtract;
 279    eval_subtract(result.real_data(), scalar);
 280 }
 281 template <class Backend, class T>
 282 inline typename boost::disable_if_c<boost::is_same<complex_adaptor<Backend>, T>::value>::type eval_multiply(complex_adaptor<Backend>& result, const T& scalar)
 283 {
 284    using default_ops::eval_multiply;
 285    eval_multiply(result.real_data(), scalar);
 286    eval_multiply(result.imag_data(), scalar);
 287 }
 288 template <class Backend, class T>
 289 inline typename boost::disable_if_c<boost::is_same<complex_adaptor<Backend>, T>::value>::type eval_divide(complex_adaptor<Backend>& result, const T& scalar)
 290 {
 291    using default_ops::eval_divide;
 292    eval_divide(result.real_data(), scalar);
 293    eval_divide(result.imag_data(), scalar);
 294 }
 295 // Optimised 3 arg versions:
 296 template <class Backend, class T>
 297 inline typename boost::disable_if_c<boost::is_same<complex_adaptor<Backend>, T>::value>::type eval_add(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& a, const T& scalar)
 298 {
 299    using default_ops::eval_add;
 300    eval_add(result.real_data(), a.real_data(), scalar);
 301    result.imag_data() = a.imag_data();
 302 }
 303 template <class Backend, class T>
 304 inline typename boost::disable_if_c<boost::is_same<complex_adaptor<Backend>, T>::value>::type eval_subtract(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& a, const T& scalar)
 305 {
 306    using default_ops::eval_subtract;
 307    eval_subtract(result.real_data(), a.real_data(), scalar);
 308    result.imag_data() = a.imag_data();
 309 }
 310 template <class Backend, class T>
 311 inline typename boost::disable_if_c<boost::is_same<complex_adaptor<Backend>, T>::value>::type eval_multiply(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& a, const T& scalar)
 312 {
 313    using default_ops::eval_multiply;
 314    eval_multiply(result.real_data(), a.real_data(), scalar);
 315    eval_multiply(result.imag_data(), a.imag_data(), scalar);
 316 }
 317 template <class Backend, class T>
 318 inline typename boost::disable_if_c<boost::is_same<complex_adaptor<Backend>, T>::value>::type eval_divide(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& a, const T& scalar)
 319 {
 320    using default_ops::eval_divide;
 321    eval_divide(result.real_data(), a.real_data(), scalar);
 322    eval_divide(result.imag_data(), a.imag_data(), scalar);
 323 }
 324
 325 template <class Backend>
 326 inline bool eval_is_zero(const complex_adaptor<Backend>& val) BOOST_NOEXCEPT
 327 {
 328    using default_ops::eval_is_zero;
 329    return eval_is_zero(val.real_data()) && eval_is_zero(val.imag_data());
 330 }
 331 template <class Backend>
 332 inline int eval_get_sign(const complex_adaptor<Backend>&)
 333 {
 334    BOOST_STATIC_ASSERT_MSG(sizeof(Backend) == UINT_MAX, "Complex numbers have no sign bit."); // designed to always fail
 335    return 0;
 336 }
 337
 338 template <class Result, class Backend>
 339 inline typename disable_if_c<boost::is_complex<Result>::value>::type eval_convert_to(Result* result, const complex_adaptor<Backend>& val)
 340 {
 341    using default_ops::eval_convert_to;
 342    using default_ops::eval_is_zero;
 343    if (!eval_is_zero(val.imag_data()))
 344    {
 345       BOOST_THROW_EXCEPTION(std::runtime_error("Could not convert imaginary number to scalar."));
 346    }
 347    eval_convert_to(result, val.real_data());
 348 }
 349
 350 template <class Backend, class T>
 351 inline void assign_components(complex_adaptor<Backend>& result, const T& a, const T& b)
 352 {
 353    result.real_data() = a;
 354    result.imag_data() = b;
 355 }
 356
 357 //
 358 // Native non-member operations:
 359 //
 360 template <class Backend>
 361 inline void eval_sqrt(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& val)
 362 {
 363    // Use the following:
 364    // sqrt(z) = (s, zi / 2s)       for zr >= 0
 365    //           (|zi| / 2s, +-s)   for zr <  0
 366    // where s = sqrt{ [ |zr| + sqrt(zr^2 + zi^2) ] / 2 },
 367    // and the +- sign is the same as the sign of zi.
 368    using default_ops::eval_abs;
 369    using default_ops::eval_add;
 370    using default_ops::eval_divide;
 371    using default_ops::eval_get_sign;
 372    using default_ops::eval_is_zero;
 373
 374    if (eval_is_zero(val.imag_data()) && (eval_get_sign(val.real_data()) >= 0))
 375    {
 376       static const typename mpl::front<typename Backend::unsigned_types>::type zero = 0u;
 377       eval_sqrt(result.real_data(), val.real_data());
 378       result.imag_data() = zero;
 379       return;
 380    }
 381
 382    const bool __my_real_part_is_neg(eval_get_sign(val.real_data()) < 0);
 383
 384    Backend __my_real_part_fabs(val.real_data());
 385    if (__my_real_part_is_neg)
 386       __my_real_part_fabs.negate();
 387
 388    Backend t, __my_sqrt_part;
 389    eval_abs(__my_sqrt_part, val);
 390    eval_add(__my_sqrt_part, __my_real_part_fabs);
 391    eval_ldexp(t, __my_sqrt_part, -1);
 392    eval_sqrt(__my_sqrt_part, t);
 393
 394    if (__my_real_part_is_neg == false)
 395    {
 396       eval_ldexp(t, __my_sqrt_part, 1);
 397       eval_divide(result.imag_data(), val.imag_data(), t);
 398       result.real_data() = __my_sqrt_part;
 399    }
 400    else
 401    {
 402       const bool __my_imag_part_is_neg(eval_get_sign(val.imag_data()) < 0);
 403
 404       Backend __my_imag_part_fabs(val.imag_data());
 405       if (__my_imag_part_is_neg)
 406          __my_imag_part_fabs.negate();
 407
 408       eval_ldexp(t, __my_sqrt_part, 1);
 409       eval_divide(result.real_data(), __my_imag_part_fabs, t);
 410       if (__my_imag_part_is_neg)
 411          __my_sqrt_part.negate();
 412       result.imag_data() = __my_sqrt_part;
 413    }
 414 }
 415
 416 template <class Backend>
 417 inline void eval_abs(Backend& result, const complex_adaptor<Backend>& val)
 418 {
 419    Backend t1, t2;
 420    eval_multiply(t1, val.real_data(), val.real_data());
 421    eval_multiply(t2, val.imag_data(), val.imag_data());
 422    eval_add(t1, t2);
 423    eval_sqrt(result, t1);
 424 }
 425
 426 template <class Backend>
 427 inline void eval_pow(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& b, const complex_adaptor<Backend>& e)
 428 {
 429    using default_ops::eval_acos;
 430    using default_ops::eval_cos;
 431    using default_ops::eval_exp;
 432    using default_ops::eval_get_sign;
 433    using default_ops::eval_is_zero;
 434    using default_ops::eval_multiply;
 435    using default_ops::eval_sin;
 436
 437    if (eval_is_zero(e))
 438    {
 439       typename mpl::front<typename Backend::unsigned_types>::type one(1);
 440       result = one;
 441       return;
 442    }
 443    else if (eval_is_zero(b))
 444    {
 445       if (eval_is_zero(e.real_data()))
 446       {
 447          Backend n          = std::numeric_limits<number<Backend> >::quiet_NaN().backend();
 448          result.real_data() = n;
 449          result.imag_data() = n;
 450       }
 451       else if (eval_get_sign(e.real_data()) < 0)
 452       {
 453          Backend n          = std::numeric_limits<number<Backend> >::infinity().backend();
 454          result.real_data() = n;
 455          typename mpl::front<typename Backend::unsigned_types>::type zero(0);
 456          if (eval_is_zero(e.imag_data()))
 457             result.imag_data() = zero;
 458          else
 459             result.imag_data() = n;
 460       }
 461       else
 462       {
 463          typename mpl::front<typename Backend::unsigned_types>::type zero(0);
 464          result = zero;
 465       }
 466       return;
 467    }
 468    complex_adaptor<Backend> t;
 469    eval_log(t, b);
 470    eval_multiply(t, e);
 471    eval_exp(result, t);
 472 }
 473
 474 template <class Backend>
 475 inline void eval_exp(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 476 {
 477    using default_ops::eval_cos;
 478    using default_ops::eval_exp;
 479    using default_ops::eval_is_zero;
 480    using default_ops::eval_multiply;
 481    using default_ops::eval_sin;
 482
 483    if (eval_is_zero(arg.imag_data()))
 484    {
 485       eval_exp(result.real_data(), arg.real_data());
 486       typename mpl::front<typename Backend::unsigned_types>::type zero(0);
 487       result.imag_data() = zero;
 488       return;
 489    }
 490    eval_cos(result.real_data(), arg.imag_data());
 491    eval_sin(result.imag_data(), arg.imag_data());
 492    Backend e;
 493    eval_exp(e, arg.real_data());
 494    if (eval_is_zero(result.real_data()))
 495       eval_multiply(result.imag_data(), e);
 496    else if (eval_is_zero(result.imag_data()))
 497       eval_multiply(result.real_data(), e);
 498    else
 499       eval_multiply(result, e);
 500 }
 501
 502 template <class Backend>
 503 inline void eval_log(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 504 {
 505    using default_ops::eval_add;
 506    using default_ops::eval_atan2;
 507    using default_ops::eval_get_sign;
 508    using default_ops::eval_is_zero;
 509    using default_ops::eval_log;
 510    using default_ops::eval_multiply;
 511
 512    if (eval_is_zero(arg.imag_data()) && (eval_get_sign(arg.real_data()) >= 0))
 513    {
 514       eval_log(result.real_data(), arg.real_data());
 515       typename mpl::front<typename Backend::unsigned_types>::type zero(0);
 516       result.imag_data() = zero;
 517       return;
 518    }
 519
 520    Backend t1, t2;
 521    eval_multiply(t1, arg.real_data(), arg.real_data());
 522    eval_multiply(t2, arg.imag_data(), arg.imag_data());
 523    eval_add(t1, t2);
 524    eval_log(t2, t1);
 525    eval_ldexp(result.real_data(), t2, -1);
 526    eval_atan2(result.imag_data(), arg.imag_data(), arg.real_data());
 527 }
 528
 529 template <class Backend>
 530 inline void eval_log10(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 531 {
 532    using default_ops::eval_divide;
 533    using default_ops::eval_log;
 534
 535    typedef typename mpl::front<typename Backend::unsigned_types>::type ui_type;
 536
 537    Backend ten;
 538    ten = ui_type(10);
 539    Backend l_ten;
 540    eval_log(l_ten, ten);
 541    eval_log(result, arg);
 542    eval_divide(result, l_ten);
 543 }
 544
 545 template <class Backend>
 546 inline void eval_sin(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 547 {
 548    using default_ops::eval_cos;
 549    using default_ops::eval_cosh;
 550    using default_ops::eval_sin;
 551    using default_ops::eval_sinh;
 552
 553    Backend t1, t2;
 554    eval_sin(t1, arg.real_data());
 555    eval_cosh(t2, arg.imag_data());
 556    eval_multiply(result.real_data(), t1, t2);
 557
 558    eval_cos(t1, arg.real_data());
 559    eval_sinh(t2, arg.imag_data());
 560    eval_multiply(result.imag_data(), t1, t2);
 561 }
 562
 563 template <class Backend>
 564 inline void eval_cos(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 565 {
 566    using default_ops::eval_cos;
 567    using default_ops::eval_cosh;
 568    using default_ops::eval_sin;
 569    using default_ops::eval_sinh;
 570
 571    Backend t1, t2;
 572    eval_cos(t1, arg.real_data());
 573    eval_cosh(t2, arg.imag_data());
 574    eval_multiply(result.real_data(), t1, t2);
 575
 576    eval_sin(t1, arg.real_data());
 577    eval_sinh(t2, arg.imag_data());
 578    eval_multiply(result.imag_data(), t1, t2);
 579    result.imag_data().negate();
 580 }
 581
 582 template <class Backend>
 583 inline void eval_tan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 584 {
 585    complex_adaptor<Backend> c;
 586    eval_cos(c, arg);
 587    eval_sin(result, arg);
 588    eval_divide(result, c);
 589 }
 590
 591 template <class Backend>
 592 inline void eval_asin(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 593 {
 594    using default_ops::eval_add;
 595    using default_ops::eval_multiply;
 596
 597    if (eval_is_zero(arg))
 598    {
 599       result = arg;
 600       return;
 601    }
 602
 603    complex_adaptor<Backend> t1, t2;
 604    assign_components(t1, arg.imag_data(), arg.real_data());
 605    t1.real_data().negate();
 606    eval_asinh(t2, t1);
 607
 608    assign_components(result, t2.imag_data(), t2.real_data());
 609    result.imag_data().negate();
 610 }
 611
 612 template <class Backend>
 613 inline void eval_acos(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 614 {
 615    typedef typename mpl::front<typename Backend::unsigned_types>::type ui_type;
 616
 617    using default_ops::eval_asin;
 618
 619    Backend half_pi, t1;
 620    t1 = static_cast<ui_type>(1u);
 621    eval_asin(half_pi, t1);
 622    eval_asin(result, arg);
 623    result.negate();
 624    eval_add(result.real_data(), half_pi);
 625 }
 626
 627 template <class Backend>
 628 inline void eval_atan(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 629 {
 630    typedef typename mpl::front<typename Backend::unsigned_types>::type ui_type;
 631    ui_type                                                             one = (ui_type)1u;
 632
 633    using default_ops::eval_add;
 634    using default_ops::eval_is_zero;
 635    using default_ops::eval_log;
 636    using default_ops::eval_subtract;
 637
 638    complex_adaptor<Backend> __my_z_times_i, t1, t2, t3;
 639    assign_components(__my_z_times_i, arg.imag_data(), arg.real_data());
 640    __my_z_times_i.real_data().negate();
 641
 642    eval_add(t1, __my_z_times_i, one);
 643    eval_log(t2, t1);
 644    eval_subtract(t1, one, __my_z_times_i);
 645    eval_log(t3, t1);
 646    eval_subtract(t1, t3, t2);
 647
 648    eval_ldexp(result.real_data(), t1.imag_data(), -1);
 649    eval_ldexp(result.imag_data(), t1.real_data(), -1);
 650    if (!eval_is_zero(result.real_data()))
 651       result.real_data().negate();
 652 }
 653
 654 template <class Backend>
 655 inline void eval_sinh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 656 {
 657    using default_ops::eval_cos;
 658    using default_ops::eval_cosh;
 659    using default_ops::eval_sin;
 660    using default_ops::eval_sinh;
 661
 662    Backend t1, t2;
 663    eval_cos(t1, arg.imag_data());
 664    eval_sinh(t2, arg.real_data());
 665    eval_multiply(result.real_data(), t1, t2);
 666
 667    eval_cosh(t1, arg.real_data());
 668    eval_sin(t2, arg.imag_data());
 669    eval_multiply(result.imag_data(), t1, t2);
 670 }
 671
 672 template <class Backend>
 673 inline void eval_cosh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 674 {
 675    using default_ops::eval_cos;
 676    using default_ops::eval_cosh;
 677    using default_ops::eval_sin;
 678    using default_ops::eval_sinh;
 679
 680    Backend t1, t2;
 681    eval_cos(t1, arg.imag_data());
 682    eval_cosh(t2, arg.real_data());
 683    eval_multiply(result.real_data(), t1, t2);
 684
 685    eval_sin(t1, arg.imag_data());
 686    eval_sinh(t2, arg.real_data());
 687    eval_multiply(result.imag_data(), t1, t2);
 688 }
 689
 690 template <class Backend>
 691 inline void eval_tanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 692 {
 693    using default_ops::eval_divide;
 694    complex_adaptor<Backend> s, c;
 695    eval_sinh(s, arg);
 696    eval_cosh(c, arg);
 697    eval_divide(result, s, c);
 698 }
 699
 700 template <class Backend>
 701 inline void eval_asinh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 702 {
 703    typedef typename mpl::front<typename Backend::unsigned_types>::type ui_type;
 704    ui_type                                                             one = (ui_type)1u;
 705
 706    using default_ops::eval_add;
 707    using default_ops::eval_log;
 708    using default_ops::eval_multiply;
 709
 710    complex_adaptor<Backend> t1, t2;
 711    eval_multiply(t1, arg, arg);
 712    eval_add(t1, one);
 713    eval_sqrt(t2, t1);
 714    eval_add(t2, arg);
 715    eval_log(result, t2);
 716 }
 717
 718 template <class Backend>
 719 inline void eval_acosh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 720 {
 721    typedef typename mpl::front<typename Backend::unsigned_types>::type ui_type;
 722    ui_type                                                             one = (ui_type)1u;
 723
 724    using default_ops::eval_add;
 725    using default_ops::eval_divide;
 726    using default_ops::eval_log;
 727    using default_ops::eval_multiply;
 728    using default_ops::eval_subtract;
 729
 730    complex_adaptor<Backend> __my_zp(arg);
 731    eval_add(__my_zp.real_data(), one);
 732    complex_adaptor<Backend> __my_zm(arg);
 733    eval_subtract(__my_zm.real_data(), one);
 734
 735    complex_adaptor<Backend> t1, t2;
 736    eval_divide(t1, __my_zm, __my_zp);
 737    eval_sqrt(t2, t1);
 738    eval_multiply(t2, __my_zp);
 739    eval_add(t2, arg);
 740    eval_log(result, t2);
 741 }
 742
 743 template <class Backend>
 744 inline void eval_atanh(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 745 {
 746    typedef typename mpl::front<typename Backend::unsigned_types>::type ui_type;
 747    ui_type                                                             one = (ui_type)1u;
 748
 749    using default_ops::eval_add;
 750    using default_ops::eval_divide;
 751    using default_ops::eval_log;
 752    using default_ops::eval_multiply;
 753    using default_ops::eval_subtract;
 754
 755    complex_adaptor<Backend> t1, t2, t3;
 756    eval_add(t1, arg, one);
 757    eval_log(t2, t1);
 758    eval_subtract(t1, one, arg);
 759    eval_log(t3, t1);
 760    eval_subtract(t2, t3);
 761
 762    eval_ldexp(result.real_data(), t2.real_data(), -1);
 763    eval_ldexp(result.imag_data(), t2.imag_data(), -1);
 764 }
 765
 766 template <class Backend>
 767 inline void eval_conj(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 768 {
 769    result = arg;
 770    result.imag_data().negate();
 771 }
 772
 773 template <class Backend>
 774 inline void eval_proj(complex_adaptor<Backend>& result, const complex_adaptor<Backend>& arg)
 775 {
 776    using default_ops::eval_get_sign;
 777
 778    typedef typename mpl::front<typename Backend::unsigned_types>::type ui_type;
 779    ui_type                                                             zero = (ui_type)0u;
 780
 781    int c1 = eval_fpclassify(arg.real_data());
 782    int c2 = eval_fpclassify(arg.imag_data());
 783    if (c1 == FP_INFINITE)
 784    {
 785       result.real_data() = arg.real_data();
 786       if (eval_get_sign(result.real_data()) < 0)
 787          result.real_data().negate();
 788       result.imag_data() = zero;
 789       if (eval_get_sign(arg.imag_data()) < 0)
 790          result.imag_data().negate();
 791    }
 792    else if (c2 == FP_INFINITE)
 793    {
 794       result.real_data() = arg.imag_data();
 795       if (eval_get_sign(result.real_data()) < 0)
 796          result.real_data().negate();
 797       result.imag_data() = zero;
 798       if (eval_get_sign(arg.imag_data()) < 0)
 799          result.imag_data().negate();
 800    }
 801    else
 802       result = arg;
 803 }
 804
 805 template <class Backend>
 806 inline void eval_real(Backend& result, const complex_adaptor<Backend>& arg)
 807 {
 808    result = arg.real_data();
 809 }
 810 template <class Backend>
 811 inline void eval_imag(Backend& result, const complex_adaptor<Backend>& arg)
 812 {
 813    result = arg.imag_data();
 814 }
 815
 816 template <class Backend, class T>
 817 inline void eval_set_imag(complex_adaptor<Backend>& result, const T& arg)
 818 {
 819    result.imag_data() = arg;
 820 }
 821
 822 template <class Backend, class T>
 823 inline void eval_set_real(complex_adaptor<Backend>& result, const T& arg)
 824 {
 825    result.real_data() = arg;
 826 }
 827
 828 template <class Backend>
 829 inline std::size_t hash_value(const complex_adaptor<Backend>& val)
 830 {
 831    std::size_t result  = hash_value(val.real_data());
 832    std::size_t result2 = hash_value(val.imag_data());
 833    boost::hash_combine(result, result2);
 834    return result;
 835 }
 836
 837 } // namespace backends
 838
 839 using boost::multiprecision::backends::complex_adaptor;
 840
 841 template <class Backend>
 842 struct number_category<complex_adaptor<Backend> > : public boost::mpl::int_<boost::multiprecision::number_kind_complex>
 843 {};
 844
 845 template <class Backend, expression_template_option ExpressionTemplates>
 846 struct component_type<number<complex_adaptor<Backend>, ExpressionTemplates> >
 847 {
 848    typedef number<Backend, ExpressionTemplates> type;
 849 };
 850
 851 template <class Backend, expression_template_option ExpressionTemplates>
 852 struct complex_result_from_scalar<number<Backend, ExpressionTemplates> >
 853 {
 854    typedef number<complex_adaptor<Backend>, ExpressionTemplates> type;
 855 };
 856
 857 }
 858
 859 } // namespace boost::multiprecision
 860
 861 #endif