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