complex_adaptor.hpp 33 KB

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