complex_adaptor.hpp 32 KB

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