complex_adaptor.hpp 27 KB

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