complex_adaptor.hpp 25 KB

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