complex_adaptor.hpp 33 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049
  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