complex_adaptor.hpp 26 KB

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