HILA
Loading...
Searching...
No Matches
cmplx.h
Go to the documentation of this file.
1#ifndef HILA_CMPLX_H_
2#define HILA_CMPLX_H_
3
4/**
5 * @file cmplx.h
6 * @brief Definition of Complex types
7 * @details This file contains definitions and methods for Complex numbers and Imaginary type.
8 *
9 * > NOTE: All overloads for operators +,-,/,* are not documented separately since there exists a
10 * > function for each combinations of scalar,imaginary and complex number representations. All
11 * > versions are documented in the Complex -- Complex definitions.
12 *
13 */
14
15
16// let's not include the std::complex
17// #include <complex>
18// #include <cmath>
19
20#include "plumbing/defs.h"
21#include "extended.h"
22
23// TEMPORARY location for vector intrinsic analogues -- result obvious
24
25template <typename T, std::enable_if_t<hila::is_arithmetic_or_extended<T>::value, int> = 0>
26inline T mul_add(T a, T b, T c) {
27 return a * b + c;
28}
29
30template <typename T, std::enable_if_t<hila::is_arithmetic_or_extended<T>::value, int> = 0>
31inline T mul_sub(T a, T b, T c) {
32 return a * b - c;
33}
34
35template <typename T, std::enable_if_t<hila::is_arithmetic_or_extended<T>::value, int> = 0>
36inline T nmul_add(T a, T b, T c) {
37 return c - a * b;
38}
39
40
41// forward define Imaginary
42template <typename T>
43class Imaginary_t;
44
45
46/**
47 * @brief Complex definition
48 * @details Define complex type as a class. This allows Hilapp to replace the internal type with
49 * a vector.
50 *
51 * __NOTE:__ T must be arithmetic and integrable. In the following documentation MyType refers to T,
52 * as in an arithmetic and integrable type.
53 * @param re Real part
54 * @param im Imaginary part
55 * @tparam T Arithmetic type
56 */
57template <typename T>
58class Complex {
59
60 static_assert((hila::is_arithmetic_or_extended<T>::value && !std::is_integral<T>::value),
61 "Complex can be used only with floating point type");
62 // This incantation is needed to make Field<Complex<>> vectorized
63
64 public:
65 using base_type = hila::arithmetic_type<T>;
66 using argument_type = T;
67
68 // and the content of the complex number
69 T re, im;
70
71 /**
72 * @brief Default constructo
73 * @details
74 *
75 * \code{.cpp}
76 * Complex<MyType> C;
77 * \endcode
78 *
79 * The default constructor initializes \p Complex#re and \p Complex#im to 0
80 *
81 */
82 Complex<T>() = default;
83 ~Complex<T>() = default;
84
85
86 /**
87 * @brief Copy constructor
88 *
89 * Initialize form already existing Complex number
90 *
91 * \code {.cpp}
92 * MyType a,b;
93 * a = hila::random();
94 * b = hila::random();
95 * Complex<MyType> C(a,b);
96 * Complex<MyType> B = C; // B.re = a, B.im = b
97 * \endcode
98 *
99 * Equivalent initializing is `Complex<MyType> B(C)`
100 *
101 * __Real constructor:__
102 *
103 * Initialize only real element and sets imaginary to 0
104 *
105 * \code {.cpp}
106 * MyType a = hila::random();
107 * Complex<MyType> C(a); // C.re = a, C.im = 0
108 * \endcode
109 *
110 * Not equivalent to `Complex<MyType> C = a`
111 *
112 * __Zero constructor:__
113 *
114 * Initialize to zero with nullpointer trick
115 *
116 * \code {.cpp}
117 * Complex<MyType> C = 0; // C.re = 0, C.im = 0
118 * \endcode
119 *
120 * Equivalent to `Complex<MyType> C` and `Complex<MyType> C(0)`
121 *
122 */
123 Complex<T>(const Complex<T> &a) = default;
124
125 // constructor from single complex --IS THIS NEEDED?
126 template <typename A, std::enable_if_t<hila::is_assignable<T &, A>::value, int> = 0>
127 constexpr Complex<T>(const Complex<A> a) : re(static_cast<T>(a.re)), im(static_cast<T>(a.im)) {}
128
129 // constructor from single scalar value
130 // Remember to mark this explicit, we do not want this to be invoked
131 // in automatic conversions (there should be methods) WHY NOT?
132
133 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
134 explicit constexpr Complex<T>(const S val) : re(val), im(0) {}
135
136 // allow construction from 0 (nullptr)
137 constexpr Complex<T>(const std::nullptr_t n) : re(0), im(0) {}
138
139 // constructor c(a,b)
140 template <
141 typename A, typename B,
142 std::enable_if_t<hila::is_assignable<T &, A>::value && hila::is_assignable<T &, B>::value,
143 int> = 0>
144 explicit constexpr Complex<T>(const A &a, const B &b) : re(a), im(b) {}
145
146 // make also std accessors real() and imag() - don't return reference, because
147 // v.real() would not then work inside loops!
148
149
150 /**
151 * @brief Real part of Complex number
152 * @details Pass by value if Complex number is defined as const
153 *
154 * @return T Real part
155 */
156
157 inline T real() const {
158 return re;
159 }
160 /**
161 * @brief Imaginary part of Complex number
162 * @details Pass by value if Complex number is defined as const
163 * @return T Imaginary part
164 */
165
166 inline T imag() const {
167 return im;
168 }
169
170 /**
171 * @brief Real part of Complex number
172 * @details Pass by reference
173 *
174 * @return T & Real part
175 */
176 inline T &real() const_function {
177 return re;
178 }
179 /**
180 * @brief Imaginary part of Complex number
181 * @details Pass by reference
182 *
183 * @return T & Imaginary part
184 */
185 inline T &imag() const_function {
186 return im;
187 }
188 // automatic casting from Complex<T> -> Complex<A>
189 // TODO: ensure this works if A is vector type!
190 // template <typename A>
191 // operator Complex<A>() const {
192 // return Complex<A>(re, im);
193 // }
194
195 /**
196 * @brief Assignment operator
197 * @details The following ways of assigning a Complex object are
198 *
199 * __Assignment for Complex:__
200 *
201 * \code {.cpp}
202 * Complex<MyType> C(hila::random(),hila::random());
203 * Complex<MyType> B;
204 * B = C; // B.re = C.re, B.im = C.im
205 * \endcode
206 *
207 * __Assignment from real:__
208 *
209 * Assigns real part and imaginary part to zero
210 *
211 * \code {.cpp}
212 * Complex<MyType> B;
213 * MyType a = hila::random;
214 * B = a; // B.re = a, B.im = 0
215 * \endcode
216 *
217 * @param s
218 * @return Complex<T>&
219 */
220 // inline Complex<T> &operator=(Complex<T> & rhs) = default;
221 inline Complex<T> &operator=(const Complex<T> &s) & = default;
222
223 // Assignment from Complex<A>
224 template <typename A, std::enable_if_t<hila::is_assignable<T&,A>::value, int> = 0>
225 inline Complex<T> &operator=(const Complex<A> &s) & {
226 re = s.re;
227 im = s.im;
228 return *this;
229 }
230
231 template <typename S, std::enable_if_t<hila::is_assignable<T&,S>::value, int> = 0>
232 inline Complex<T> &operator=(S s) & {
233 re = s;
234 im = 0;
235 return *this;
236 }
237
238 // delete the assign to rvalue
239 template <typename S>
240 Complex<T> &operator=(const S &s) && = delete;
241
242 /**
243 * @brief Compute square norm of Complex number
244 * @details
245 * \f{align}{ |z|^2 = \Re(z)^2 + \Im(z)^2\f}
246 *
247 * @return T Squared norm
248 */
249 inline T squarenorm() const {
250 return re * re + im * im;
251 }
252
253 /**
254 * @brief Compute absolute value of Complex number
255 * @details Essentially sqrt(squarenorm(z)):
256 * \f{align}{ |z| = \sqrt{\Re(z)^2 + \Im(z)^2}\f}
257 *
258 * @return T
259 */
260 inline T abs() const {
261 return sqrt(squarenorm());
262 }
263
264 /**
265 * @brief Compute argument of Complex number
266 * @details
267 * \f{align}{ \arg(z) = \arctan2(\Im(z)),\Re(z)\f}
268 *
269 * @return T
270 */
271 inline T arg() const {
272 return atan2(im, re);
273 }
274
275 /**
276 * @brief Compute conjugate of Complex number
277 * @details
278 * \f{align}{ z &= x + i\cdot y \\
279 * \Rightarrow z^* &= x - i\cdot y\f}
280 *
281 * @return Complex<T>
282 */
283 inline Complex<T> conj() const {
284 return Complex<T>(re, -im);
285 }
286
287 /**
288 * @brief Compute dagger of Complex number
289 * @details Alias to Complex::conj
290 *
291 * \f{align}{ z^* = z^\dagger \f}
292 * @return Complex<T>
293 */
294 inline Complex<T> dagger() const {
295 return Complex<T>(re, -im);
296 }
297
298 /**
299 * @brief Stores and returns Complex number given in polar coordinates
300 *
301 * \f{align}{ z = r\cdot e^{i\theta} \f}
302 * \code {.cpp}
303 * Complex<double> c;
304 * double r = 1;
305 * double theta = 3.14159265358979/2; // pi/2
306 * c.polar(r,theta); // c.re = 0, c.im = 1
307 * \endcode
308 *
309 *
310 * @param r Radius of Complex number
311 * @param theta Angle of complex number in radians
312 * @return Complex<T> Complex number
313 */
314 inline Complex<T> polar(const T r, const T theta) out_only {
315 re = r * cos(theta);
316 im = r * sin(theta);
317 return *this;
318 }
319
320 /**
321 * @brief Assign random values to Complex real and imaginary part
322 * @details Uses hila::random for both real and imaginary part
323 *
324 * @return Complex<T>&
325 */
326 inline Complex<T> &random() out_only {
327 re = hila::random();
328 im = hila::random();
329 return *this;
330 }
331
332 /**
333 * @brief Produces complex gaussian random values
334 * @details Uses hila::gaussrand2 for both real and imaignary part
335 * Assigns same random value for both real and imaginary part
336 * @param width gaussian_random
337 * @return Complex<T>&
338 */
339 inline Complex<T> &gaussian_random(double width = 1.0) out_only {
340 double d;
341 re = hila::gaussrand2(d) * width;
342 im = d * width;
343 return *this;
344 }
345
346 /**
347 * @brief Unary + operator
348 * @details Identity operation
349 *
350 * \code {.cpp}
351 * Complex<MyType> z;
352 * z == +z;
353 * \endcode
354 *
355 * @return Complex<Tr>
356 */
357 inline Complex<T> operator+() const {
358 return *this;
359 }
360
361 /**
362 * @brief Unary - operator
363 * @details Negates Complex number
364 *
365 * \code {.cpp}
366 * Complex<MyType> z;
367 * z = -z; // Results in -Re.z and -Im.z
368 * \endcode
369 *
370 * @return Complex<T>
371 */
372 inline Complex<T> operator-() const {
373 return Complex<T>(-re, -im);
374 }
375
376 // mark += and *= as loop_functions, because these can be used
377 // in reductions implicitly
378
379 /**
380 * @brief Complex addition assignment operator
381 * @details
382 *
383 * \code{.cpp}
384 * Complex<double> z(0,0);
385 * Complex<double> w(1,1);
386 * z += w; // z.re = 1, z.im = 1
387 * \endcode
388 *
389 * @tparam A Typpe of Complex number to add
390 * @param lhs Complex number to add
391 * @return Complex<T>&
392 */
393 template <typename A>
394 inline Complex<T> &operator+=(const Complex<A> &lhs) {
395 re += lhs.re;
396 im += lhs.im;
397 return *this;
398 }
399
400 /**
401 * @brief Real addition assignment operator
402 *
403 * Add assign only to real part of Complex number
404 *
405 * \code {.cpp}
406 * Complex<double> z(0,0);
407 * z += 1; // z.re = 1, z.im = 0
408 * \endcode
409 *
410 *
411 */
412
413 template <typename A, std::enable_if_t<hila::is_assignable<T &, A>::value, int> = 0>
414 inline Complex<T> &operator+=(const A &a) & {
415 re += a;
416 return *this;
417 }
418
419 /**
420 * @brief Complex subtraction assignment operator
421 * @details
422 *
423 * \code{.cpp}
424 * Complex<double> z(0,0);
425 * Complex<double> w(1,1);
426 * z -= w; // z.re = -1, z.im = -1
427 * \endcode
428 *
429 * @tparam A Type for complex number to subtract
430 * @param lhs Complex number to subtract
431 * @return Complex<T>&
432 */
433 template <typename A>
434 inline Complex<T> &operator-=(const Complex<A> &lhs) & {
435 static_assert(!(hila::is_arithmetic<T>::value && hila::is_extended<A>::value),
436 "Cannot assign ExtendedPrecision value to non-EP variable without casting");
437 re -= lhs.re;
438 im -= lhs.im;
439 return *this;
440 }
441
442 template <typename A, std::enable_if_t<hila::is_arithmetic_or_extended<A>::value, int> = 0>
443 inline Complex<T> &operator-=(const A &a) & {
444 static_assert(!(hila::is_arithmetic<T>::value && hila::is_extended<A>::value),
445 "Cannot assign ExtendedPrecision value to non-EP variable without casting");
446 re -= a;
447 return *this;
448 }
449
450 // inline Complex<T> & operator*= (const Complex<T> & lhs) {
451 // T r = re * lhs.re - im * lhs.im;
452 // im = im * lhs.re + re * lhs.im;
453 // re = r;
454 // return *this;
455 // }
456
457 /**
458 * @brief Complex multiply assignment operator
459 * @details
460 * Standard Complex number multiplication
461 *
462 * \f{align}{z &= x + iy, w = x' + iy' \\
463 * z w &= (x + iy)(x' + iy') = (xx'-yy') + i(xy' + yx')\f}
464 * \code{.cpp}
465 * Complex<double> z,w;
466 * //
467 * // z,w get values
468 * //
469 * z*=w; // z = zw as defined above
470 * \endcode
471 *
472 * @tparam A Type for Complex number multiply
473 * @param lhs Complex number to multiply
474 * @return Complex<T>&
475 */
476 template <typename A>
477 inline Complex<T> &operator*=(const Complex<A> &lhs) & {
478 T r = mul_sub(re, lhs.re, im * lhs.im); // a*b-c
479 im = mul_add(im, lhs.re, re * lhs.im); // a*b+c
480 re = r;
481 return *this;
482 }
483
484 /**
485 * @brief Real multiply assign operator
486 * @details Multiply assign by real number to both components of Complex number
487 *
488 * \code {.cpp}
489 * Complex<double> z(1,1);
490 * z *= 2; // z.re = 2, z.im = 2
491 * \endcode
492 * @tparam A Type for real value a
493 * @param a Real value to multiply
494 * @return Complex<T>&
495 */
496 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
497 inline Complex<T> &operator*=(const A a) & {
498 re *= a;
499 im *= a;
500 return *this;
501 }
502
503 // a/b = a b*/|b|^2 = (a.re*b.re + a.im*b.im + i(a.im*b.re - a.re*b.im))/|b|^2
504 // inline Complex<T> & operator/= (const Complex<T> & lhs) {
505 // T n = lhs.squarenorm();
506 // T r = (re * lhs.re + im * lhs.im)/n;
507 // im = (im * lhs.re - re * lhs.im)/n;
508 // re = r;
509 // return *this;
510 // }
511
512 /**
513 * @brief Complex divide assignment operator
514 * @details
515 *
516 * Standard Complex number division
517 *
518 * \f{align}{z &= x + iy, w = x' + iy' \\
519 * \frac{z}{w} &= \frac{x + iy}{x' + iy'} = \frac{(xx'+ yy') + i( yx' - xy')}{|w|^2}\f}
520 * \code{.cpp}
521 * Complex<double> z,w;
522 * //
523 * // z,w get values
524 * //
525 * z/=w; // z = z/w as defined above
526 * \endcode
527 *
528 * @tparam A Type for Complex number divide
529 * @param lhs Complex number to divide
530 * @return Complex<T>&
531 */
532 template <typename A>
533 inline Complex<T> &operator/=(const Complex<A> &lhs) & {
534 T n = lhs.squarenorm();
535 T r = mul_add(re, lhs.re, im * lhs.im) / n; // a*b+c
536 im = mul_sub(im, lhs.re, re * lhs.im) / n; // a*b-c
537 re = r;
538 return *this;
539 }
540
541 /**
542 * @brief Real divide assignment operator
543 * @details Divide assign by real number to both components of Complex number
544 *
545 * \code {.cpp}
546 * Complex<double> z(2,2);
547 * z /= 2; // z.re = 1, z.im = 1
548 * \endcode
549 * @tparam A Type for real value to divide
550 * @param a Real value to divide
551 * @return Complex<T>&
552 */
553 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
554 inline Complex<T> &operator/=(const A &a) & {
555 re /= a;
556 im /= a;
557 return *this;
558 }
559
560 // define also increment and decrement operators (though not so intuitive for complex)
561
562 /**
563 * @brief Single increment operator
564 *
565 * Increments real part of Complex number
566 *
567 * \code {.cpp}
568 * Complex<double> z(1,1);
569 * z++; // z.re = 2, z.im = 1
570 * \endcode
571 *
572 *
573 * @return Complex<T>&
574 */
576 this->re++;
577 return *this;
578 }
579 /**
580 * @brief Single decrement operator
581 *
582 * Decrement real part of Complex number
583 *
584 * \code {.cpp}
585 * Complex<double> z(1,1);
586 * z--; // z.re = 0, z.im = 1
587 * \endcode
588 *
589 *
590 * @return Complex<T>&
591 */
593 this->re--;
594 return *this;
595 }
596 /**
597 * @brief Multiple increment operator
598 *
599 * Increments real part of Complex number
600 *
601 * \code {.cpp}
602 * Complex<double> z(2,1);
603 * z++2; // z.re = 4, z.im = 1
604 * \endcode
605 *
606 *
607 * @return Complex<T>&
608 */
609 inline Complex<T> operator++(int) {
610 Complex<T> a = *this;
611 this->re++;
612 return a;
613 }
614
615 /**
616 * @brief Multiple decrement operator
617 *
618 * Decrement real part of Complex number
619 *
620 * \code {.cpp}
621 * Complex<double> z(2,1);
622 * z--2; // z.re = 0, z.im = 1
623 * \endcode
624 *
625 *
626 * @return Complex<T>&
627 */
628 inline Complex<T> operator--(int) {
629 Complex<T> a = *this;
630 this->re--;
631 return a;
632 }
633
634 std::string str(int prec = 8, char separator = ' ') const {
635 std::stringstream ss;
636 ss.precision(prec);
637 ss << re << separator << im;
638 return ss.str();
639 }
640
641
642 // Convenience method a.conj_mul(b) == a^* b
643
644 /**
645 * @brief Conjugate multiply method
646 * @details Conjugate a (*this) and multiply with give argument b: `a.conj_mul(b)` \f$ \equiv
647 * a^* \cdot b\f$
648 *
649 * @tparam A Type for Complex number b
650 * @param b number to multiply with
651 * @return Complex<T>
652 */
653 template <typename A>
654 inline Complex<T> conj_mul(const Complex<A> &b) const {
655 return Complex<T>(re * b.re + im * b.im, re * b.im - im * b.re);
656 }
657
658 // Convenience method a.mul_conj(b) == a * b^*
659
660 /**
661 * @brief Multiply conjugate method
662 * @details Multiply a (*this) with conjugate of given argument b: `a.mul_conj(b)` \f$ \equiv a
663 * \cdot b^*\f$
664 *
665 * @tparam A Type for Complex number b
666 * @param b number to conjugate and multiply withv
667 * @return Complex<T>
668 */
669 template <typename A>
670 inline Complex<T> mul_conj(const Complex<A> &b) const {
671 return Complex<T>(re * b.re + im * b.im, im * b.re - re * b.im);
672 }
673
674 // // cast to another number type (IS THIS NEEDED FOR COMPLEX?)
675 // template <typename Ntype>
676 // Complex<Ntype> cast_to() const {
677 // Complex<Ntype> res;
678 // res.re = static_cast<Ntype>(re);
679 // res.im = static_cast<Ntype>(im);
680 // return res;
681 // }
682};
683
684
685namespace hila {
686////////////////////////////////////////////////////////////////////////
687// Section for adding complex utility templates (in namespace hila)
688// hila::is_complex_or_arithmetic<T>::value
689// hila::contains_complex<T>::value
690// hila::number_type<T> returns Complex or arithmetic type
691// hila::ntype_op<A,B> returns the conventionally upgraded complex or scalar number type
692// hila::complex_x_scalar_type<A,B> type of operation Complex<A> * scalar<B>
693
694////////////////////////////////////////////////////////////////////////
695// Define hila::is_complex<T>::value -template, using specialization
696template <typename T>
697struct is_complex : std::integral_constant<bool, false> {};
698
699template <typename T>
700struct is_complex<Complex<T>> : std::integral_constant<bool, true> {};
701
702template <typename T>
703struct is_complex<Imaginary_t<T>> : std::integral_constant<bool, true> {};
704
705/// hila::is_complex_or_arithmetic<T>::value
706template <typename T>
708 : std::integral_constant<bool, hila::is_arithmetic<T>::value || hila::is_complex<T>::value> {};
709
710/////////////////////////////////////////////////////////////////////////
711// Utility to check that the type contains complex numbers
712// Use as contains_complex<T>::value
713
714template <typename T, typename Enable = void>
715struct contains_complex : std::integral_constant<bool, false> {};
716
717template <typename T>
718struct contains_complex<T, typename std::enable_if_t<hila::is_field_class_type<T>::value>>
719 : std::integral_constant<bool,
720 hila::contains_type<T, Complex<hila::arithmetic_type<T>>>::value> {};
721
722/////////////////////////////////////////////////////////////////////////
723// Utility hila::number_type<T> returns complex or arithmetic type
724// depending on the type
725
726template <typename T, typename Enable = void, typename N = hila::arithmetic_type<T>>
727struct complex_or_arithmetic_type_struct {
728 using type = N;
729};
730
731template <typename T>
732struct complex_or_arithmetic_type_struct<
733 T, typename std::enable_if_t<hila::contains_complex<T>::value>> {
735};
736
737template <typename T>
738using number_type = typename complex_or_arithmetic_type_struct<T>::type;
739
740/////////////////////////////////////////////////////////////////////////
741// Utility hila::ntype_op<T1,T2> returns real or complex type,
742// "conventionally" defined double/float etc. upgrading.
743// Note: hila::type_plus<T1,T2> gives e.g. Complex<float>,double -> Complex<float>
744
745template <typename T1, typename T2, typename Enable = void>
746struct ntype_op_s {
747 using type = hila::type_plus<hila::arithmetic_type<T1>, hila::arithmetic_type<T2>>;
748};
749
750// if one type is complex?
751template <typename T1, typename T2>
752struct ntype_op_s<T1, T2,
753 typename std::enable_if_t<(hila::contains_complex<T1>::value ||
754 hila::contains_complex<T2>::value)>> {
755 using type = Complex<hila::type_plus<hila::arithmetic_type<T1>, hila::arithmetic_type<T2>>>;
756};
757
758template <typename T1, typename T2>
759using ntype_op = typename ntype_op_s<T1, T2>::type;
760
761////////////////////////////////////////////////////////////////////////////
762// utility complex_x_scalar, used in Complex<A> + scalar<B> ops
763// - if result is convertible to Complex<A>, return Complex<A>
764// - otherwise return Complex<type_plus<A,B>>
765// This is done in order to keep vector types as vectors
766
767template <typename A, typename B, typename Enable = void>
768struct complex_x_scalar_s {
769 using type = Complex<hila::type_plus<A, B>>;
770};
771
772template <typename A, typename B>
773struct complex_x_scalar_s<
774 A, B, typename std::enable_if_t<std::is_assignable<A &, hila::type_plus<A, B>>::value>> {
775 using type = Complex<A>;
776};
777
778template <typename A, typename B>
779using complex_x_scalar_type = typename complex_x_scalar_s<A, B>::type;
780
781////////////////////////////////////////////////////////////////////////
782// as_complex_array(T var)
783// casts the var to a pointer to complex<scalar_type<T>>
784// assuming var contains complex type. This enables access
785// of complex elements as
786// as_complex_array(var)[i]
787// comment out as hilapp gets confused at the moment
788////////////////////////////////////////////////////////////////////////
789
790// template <typename T, std::enable_if_t<hila::contains_complex<T>::value, int> = 0>
791// inline const Complex<hila::arithmetic_type<T>> * as_complex_array(const T &var) {
792// return (const Complex<hila::arithmetic_type<T>> *)(void *)&var;
793// }
794
795// template <typename T, std::enable_if_t<hila::contains_complex<T>::value, int> = 0>
796// inline Complex<hila::arithmetic_type<T>> * as_complex_array(T &var) {
797// return (Complex<hila::arithmetic_type<T>> *)(void *)&var;
798// }
799
800////////////////////////////////////////////////////////////////////////
801// get_complex_element(var,i)
802// returns the i:th complex number embedded in variable v
803// set_complex_element(var,i,value)
804// sets the i:th element in var
805////////////////////////////////////////////////////////////////////////
806
807template <typename T, std::enable_if_t<hila::contains_complex<T>::value, int> = 0>
808inline Complex<hila::arithmetic_type<T>> get_complex_in_var(const T &var, int i) {
809 return *(reinterpret_cast<const Complex<hila::arithmetic_type<T>> *>(&var) + i);
810}
811
812template <typename T, std::enable_if_t<hila::contains_complex<T>::value, int> = 0>
813inline void set_complex_in_var(T &var, int i, const Complex<hila::arithmetic_type<T>> &val) {
814 *(reinterpret_cast<Complex<hila::arithmetic_type<T>> *>(&var) + i) = val;
815}
816
817} // namespace hila
818
819// generic Complex constructor - type from arguments
820template <typename T, std::enable_if_t<hila::is_floating_point<T>::value, int> = 0>
821Complex<T> complex(const T re, const T im) {
822 return Complex<T>(re, im);
823}
824
825
826/**
827 * @brief Return real value of Complex number
828 * @details Wrapper around Complex::real
829 * @tparam T Arithmetic type of a
830 * @param a Complex number to get real value from
831 * @return T
832 */
833template <typename T>
834inline T real(const Complex<T> &a) {
835 return a.re;
836}
837
838/**
839 * @brief Retrun imaginary value of Complex number
840 * @details Wrapper around Complex::imag
841 * @tparam T Arithmetic type of a
842 * @param a Complex number to get imaginary value from
843 * @return T
844 */
845template <typename T>
846inline T imag(const Complex<T> &a) {
847 return a.im;
848}
849
850/**
851 * @brief Return complex number given by polar representation
852 * @details Same as Complex::polar
853 * @tparam T Arithmetic type r and arg
854 * @param r Radial component of Complex number
855 * @param arg Argument of Complex number
856 * @return Complex<T>
857 */
858template <typename T>
860 Complex<T> res(r * cos(arg), r * sin(arg));
861 return res;
862}
863
864
865// template <typename T>
866// inline Complex<T> operator+(const Complex<T> & a, const Complex<T> & b) {
867// return Complex<T>(a.re + b.re, a.im + b.im);
868// }
869
870
871#pragma hila loop_function
872/**
873 * @brief Addition operator Complex + Complex
874 * @memberof Complex
875 * @details Addition between Complex numbers is defined in the usual way:
876 *
877 * \f{align}{ z_1, z_2 &\in \mathbf{C} \\ z_1 + z_2 &= \Re(z_1) + \Re(z_2) + [\Im(z_1) + \Im(z_2)]i
878 * \f}
879 *
880 * \code {.cpp}
881 * z = z_1 + z_2;
882 * \endcode
883 *
884 * @tparam T1 Type for first complex number
885 * @tparam T2 Type for second complex number
886 * @tparam Tr Conversion type between T1 and T2
887 * @param a First complex number
888 * @param b Second complex number
889 * @return Complex<Tr>
890 */
891template <typename T1, typename T2, typename Tr = hila::type_plus<T1, T2>>
892inline Complex<Tr> operator+(const Complex<T1> &a, const Complex<T2> &b) {
893 return Complex<Tr>(a.re + b.re, a.im + b.im);
894}
895
896// TODO: for avx vector too -- #define new template macro
897/**
898 * @brief Addition operator Complex + Scalar
899 * @memberof Complex
900 * @details Defined in the usual way. \f$z,z_1 \in \mathbb{C}, a \in \mathbb{R}\f$
901 * \code {.cpp}
902 * z = z_1 + a;
903 * \endcode
904 * @tparam T Arithmetic type of c
905 * @tparam A Arithmetic type of a
906 * @param c Complex number to sum
907 * @param a Scalar to sum
908 * @return auto
909 */
910template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
911inline auto operator+(const Complex<T> &c, const A &a) {
912 return hila::complex_x_scalar_type<T, A>(c.re + a, c.im);
913}
914
915/**
916 * @brief Addition operator Scalar + Complex
917 * @memberof Complex
918 * @details Defined in the usual way. \f$z,z_1 \in \mathbb{C}, a \in \mathbb{R}\f$
919 * \code {.cpp}
920 * z = a + z_1;
921 * \endcode
922 * @tparam T Arithmetic type of c
923 * @tparam A Arithmetic type of a
924 * @param a Scalar to sum
925 * @param c Complex number to sum
926 * @return auto
927 */
928template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
929inline auto operator+(const A &a, const Complex<T> &c) {
930 return hila::complex_x_scalar_type<T, A>(c.re + a, c.im);
931}
932
933/**
934 * @brief Subtraction operator Complex - Complex
935 * @memberof Complex
936 * @details Subtraction between Complex numbers is defined in the usual way
937 * \f{align}{ z_1, z_2 &\in \mathbf{C} \\ z_1 - z_2 &= \Re(z_1) - \Re(z_2) + [\Im(z_1) -
938 * \Im(z_2)]i\f}
939 *
940 * \code {.cpp}
941 * z = z_1 - z_2;
942 * \endcode
943 *
944 * @tparam T1 Arithmetic type of a
945 * @tparam T2 Arithmetic type of b
946 * @tparam Tr Resulting type after subtraction
947 * @param a First complex number
948 * @param b Second complex number to subtract
949 * @return Complex<Tr>
950 */
951template <typename T1, typename T2, typename Tr = hila::type_plus<T1, T2>>
952inline Complex<Tr> operator-(const Complex<T1> &a, const Complex<T2> &b) {
953 return Complex<Tr>(a.re - b.re, a.im - b.im);
954}
955
956/**
957 * @brief Subtraction operator Complex - Scalar
958 * @memberof Complex
959 * @details Defined in the usual way. \f$z,z_1 \in \mathbb{C}, a \in \mathbb{R}\f$
960 * \code {.cpp}
961 * z = z_1 - a;
962 * \endcode
963 * @tparam T Arithmetic type of c
964 * @tparam A Arithmetic type of a
965 * @param c Complex number
966 * @param a Scalar to subtract
967 * @return auto
968 */
969template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
970inline auto operator-(const Complex<T> &c, const A &a) {
971 return hila::complex_x_scalar_type<T, A>(c.re - a, c.im);
972}
973
974/**
975 * @memberof Complex
976 * @brief Subtraction operator Scalar - Complex
977 * @details Defined in the usual way. \f$z,z_1 \in \mathbb{C}, a \in \mathbb{R}\f$
978 * \code {.cpp}
979 * z = a - z_1;
980 * \endcode
981 * @tparam T Arithmetic type of c
982 * @tparam A Arithmetic type of a
983 * @param a Scalar
984 * @param c Complex number to subtract
985 * @return auto
986 */
987template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
988inline auto operator-(const A &a, const Complex<T> &c) {
989 return hila::complex_x_scalar_type<T, A>(a - c.re, -c.im);
990}
991
992/**
993 * @brief Multiplication operator Complex * Complex
994 * @memberof Complex
995 * @details Defined in the usual way
996 *
997 * \f{align}{z, z_1 &= x + iy, z_2 = x' + iy' \in \mathbb{C} \\
998 * z_1 z_2 &= (x + iy)(x' + iy') = (xx'-yy') + i(xy' + yx')\f}
999 * \code {.cpp}
1000 * z = z_1 * z_2;
1001 * \endcode
1002 * @tparam T1 Arithmetic type of a
1003 * @tparam T2 Arithmetic type of b
1004 * @tparam Tr Resulting type after multiplication
1005 * @param a
1006 * @param b
1007 * @return Complex<Tr>
1008 */
1009template <typename T1, typename T2, typename Tr = hila::type_mul<T1, T2>>
1010inline Complex<Tr> operator*(const Complex<T1> &a, const Complex<T2> &b) {
1011 return Complex<Tr>(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re);
1012}
1013
1014/**
1015 * @brief Multiplication operator Complex * Scalar
1016 * @memberof Complex
1017 * @details Multiplication between Complex and scalar is defined in the usual way
1018 * \f{align}{z,z_1 &= x + iy \in \mathbb{C}, a \in \mathbb{R} \\
1019 * z_1 * a &= (x\cdot a + iy\cdot a)\f}
1020 * \code {.cpp}
1021 * z = z_1 * a;
1022 * \endcode
1023 * @tparam T Arithmetic type of c
1024 * @tparam A Arithmetic type of a
1025 * @param c Complex number to multiply
1026 * @param a Scalar to multiply
1027 * @return auto
1028 */
1029template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
1030inline auto operator*(const Complex<T> &c, const A &a) {
1031 return hila::complex_x_scalar_type<T, A>(c.re * a, c.im * a);
1032}
1033
1034/**
1035 * @brief Multiplication operator Scalar * Complex
1036 * @memberof Complex
1037 * @details Multiplication between Complex and scalar is defined in the usual way
1038 * \f{align}{z,z_1 &= x + iy \in \mathbb{C}, a \in \mathbb{R} \\
1039 * z * a &= (x\cdot a + iy\cdot a)\f}
1040 * \code {.cpp}
1041 * z = a * z_1;
1042 * \endcode
1043 * @tparam T Arithmetic type of c
1044 * @tparam A Arithmetic type of a
1045 * @param c Complex number to multiply
1046 * @param a Scalar to multiply
1047 * @return auto
1048 */
1049template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
1050inline auto operator*(const A &a, const Complex<T> &c) {
1051 return hila::complex_x_scalar_type<T, A>(a * c.re, a * c.im);
1052}
1053
1054// / a/b = ab*/| b |^2
1055// template <typename T>
1056// inline Complex<T> operator/(const Complex<T> & a, const Complex<T> & b) {
1057// T n = b.squarenorm();
1058// return Complex<T>( (a.re*b.re + a.im*b.im)/n, (a.im*b.re - a.re*b.im)/n );
1059// }
1060
1061/**
1062 * @brief Division operator Complex / Complex
1063 * @memberof Complex
1064 * @details Defined in the usual way
1065 * \f{align}{ z,z_1,z_2 &\in \mathbb{C} \\
1066 * \frac{z_1}{z_2} &= \frac{z_1\cdot z_2^{*}}{|z_2|^2} \f}
1067 * \code {.cpp}
1068 * z =z_1 / z_2;
1069 * \endcode
1070 * @tparam T1 Arithmetic type of a
1071 * @tparam T2 Arithmetic type of b
1072 * @tparam Tr Resulting type after division
1073 * @param a Complex number to divide
1074 * @param b Complex number to divide with
1075 * @return Complex<Tr>
1076 */
1077template <typename T1, typename T2, typename Tr = hila::type_mul<T1, T2>>
1078inline Complex<Tr> operator/(const Complex<T1> &a, const Complex<T2> &b) {
1079 T2 n = 1.0 / b.squarenorm();
1080 return Complex<Tr>((a.re * b.re + a.im * b.im) * n, (a.im * b.re - a.re * b.im) * n);
1081}
1082
1083/**
1084 * @memberof Complex
1085 * @brief Division operator Complex / Scalar
1086 * @details Defined in the usual way
1087 * \f{align}{z, z_1=x + iy &\in \mathbb{C}, a \in \mathbb{R} \\
1088 * \frac{z_1}{a} &= \frac{x}{a} + i\cdot \frac{y}{a} \f}
1089 * \code {.cpp}
1090 * z = z_1 / a;
1091 * \endcode
1092 * @tparam T Arithmetic type of c
1093 * @tparam A Arithmetic type of a
1094 * @param c Complex number to divide
1095 * @param a Scalar to divide with
1096 * @return auto
1097 */
1098template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
1099inline auto operator/(const Complex<T> &c, const A &a) {
1100 return hila::complex_x_scalar_type<T, A>(c.re / a, c.im / a);
1101}
1102
1103// a/c = ac*/|c|^2
1104/**
1105 * @memberof Complex
1106 * @brief Division operator Scalar / Complex
1107 * @details Defined in the usual way
1108 * \f{align}{ z,z_1=x + iy &\in \mathbb{C}, a \in \mathbb{R} \\
1109 * \frac{a}{z_1} &= \frac{az^*}{|z|^2} \f}
1110 * \code {.cpp}
1111 * z = a / z_1;
1112 * \endcode
1113 * @tparam T Arithmetic type of c
1114 * @tparam A Arithmetic type of a
1115 * @param c Complex number to divide
1116 * @param a Scalar to divide with
1117 * @return auto
1118 */
1119template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
1120inline auto operator/(const A &a, const Complex<T> &c) {
1121 T n = c.squarenorm();
1122 return hila::complex_x_scalar_type<T, A>((a * c.re) / n, -(a * c.im) / n);
1123}
1124
1125/**
1126 * @brief Multiply add with Complex numbers
1127 * @memberof Complex
1128 * @details Defined as
1129 *
1130 * \code{.cpp}
1131 * mul_add(a,b,c) = a*b + c;
1132 * \endcode
1133 * @tparam T Arithmetic type of a,b and c
1134 * @param a Complex number to multiply
1135 * @param b Complex number to multiply
1136 * @param c Complex number to add to result of multiplication
1137 * @return Complex<T>
1138 */
1139template <typename T>
1140inline Complex<T> mul_add(const Complex<T> &a, const Complex<T> &b, const Complex<T> &c) {
1141 // a*b + c
1142 Complex<T> r;
1143 T t1 = mul_add(a.re, b.re, c.re);
1144 T t2 = mul_add(a.re, b.im, c.im);
1145 r.re = mul_add(a.im, b.im, t1); // -a.im*b.im + a.re*b.re + c.re
1146 r.im = mul_add(a.im, b.re, t2); // a.im*b.re + a.re*b.im + c.im
1147 return r;
1148}
1149
1150/**
1151 * @brief Compare equality of two complex numbers
1152 * @memberof Complex
1153 * Two numbers are equal, if the real and imaginary components are respectively equal
1154 *
1155 * @tparam A Arithmetic type of a
1156 * @tparam B Arithmetic type of b
1157 * @param a Complex number to compare
1158 * @param b Complex number to compare
1159 * @return true if values compare to equal
1160 */
1161
1162template <typename A, typename B>
1163inline bool operator==(const Complex<A> &a, const Complex<B> &b) {
1164 return (a.re == b.re && a.im == b.im);
1165}
1166
1167/**
1168 * @memberof Complex
1169 * @brief Compare equality of Complex and scalar
1170 *
1171 * Two numbers are equal, if the arithmetic values are equal: thus,
1172 * complex and real comparison (a + i b) == a
1173 * is true if b == 0.
1174 *
1175 * @tparam A Arithmetic type of a
1176 * @tparam B Arithmetic type of b
1177 * @param a Complex number to compare
1178 * @param b Scalar to compare
1179 * @return true if values compare to equal
1180 */
1181template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<B>::value, int> = 0>
1182inline bool operator==(const Complex<A> &a, const B b) {
1183 return (a.re == b && a.im == 0);
1184}
1185
1186/**
1187 * @memberof Complex
1188 * @brief Compare equality of Scalar and Complex
1189 *
1190 * Two numbers are equal, if the arithmetic values are equal: thus,
1191 * complex and real comparison (a + i b) == a
1192 * is true if b == 0.
1193 *
1194 * @tparam A Arithmetic type of a
1195 * @tparam B Arithmetic type of b
1196 * @param b Complex number to compare
1197 * @param a Scalar to compare
1198 * @return true if values compare to equal
1199 */
1200template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
1201inline bool operator==(const A a, const Complex<B> &b) {
1202 return b == a;
1203}
1204
1205/**
1206 * @brief Compare non-equality of two complex numbers
1207 * @memberof Complex
1208 * @details Negation of operator==()
1209 *
1210 * @tparam A Arithmetic type of a
1211 * @tparam B Arithmetic type of b
1212 * @param a Complex number to compare
1213 * @param b Complex number to compare
1214 * @return true if values are not arithmetically equal
1215 */
1216template <typename A, typename B>
1217inline bool operator!=(const Complex<A> &a, const Complex<B> &b) {
1218 return (a.re != b.re || a.im != b.im);
1219}
1220
1221/**
1222 * @brief Compare non-equality of Complex number and Scalar
1223 * @memberof Complex
1224 * @details Negation of operator==()
1225 *
1226 * @tparam A Arithmetic type of a
1227 * @tparam B Arithmetic type of b
1228 * @param a Complex number to compare
1229 * @param b Scalar to compare
1230 * @return true if values are not arithmetically equal
1231 */
1232template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<B>::value, int> = 0>
1233inline bool operator!=(const Complex<A> &a, const B b) {
1234 return (a.re != b || a.im != 0);
1235}
1236
1237/**
1238 * @brief Compare non-equality of Scalar and Complex number
1239 * @memberof Complex
1240 * @details Negation of operator==()
1241 *
1242 * @tparam A Arithmetic type of a
1243 * @tparam B Arithmetic type of b
1244 * @param a Scalar to compare
1245 * @param b Complex number to compare
1246 * @return true if values are not arithmetically equal
1247 */
1248template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
1249inline bool operator!=(const A a, const Complex<B> &b) {
1250 return b != a;
1251}
1252
1253
1254
1255//////////////////////////////////////////////////////////////////////////////////
1256// Some operations in function form. Useful in templates when the arg type is not known
1257
1258/**
1259 * @brief Return absolute value of Complex number
1260 * @details Wrapper around Complex::abs
1261 * @tparam T Arithmetic type of a
1262 * @param a Complex number to get abs from
1263 * @return T
1264 */
1265template <typename T>
1266inline T abs(const Complex<T> &a) {
1267 return a.abs();
1268}
1269
1270/**
1271 * @brief Return argument of Complex number
1272 * @details Wrapper around Complex::arg
1273 * @tparam T Arithmetic type of a
1274 * @param a Complex number to get arg from
1275 * @return T
1276 */
1277template <typename T>
1278inline T arg(const Complex<T> &a) {
1279 return a.arg();
1280}
1281
1282/**
1283 * @brief Return conjugate of Complex number
1284 * @details Wrapper around Complex::conj
1285 * @tparam T Arithmetic type of a
1286 * @param a Complex number to get conjugate from
1287 * @return Complex<T>
1288 */
1289template <typename T>
1290inline Complex<T> conj(const Complex<T> &val) {
1291 return val.conj();
1292}
1293
1294/**
1295 * @brief Return dagger of Complex number
1296 * @details Wrapper around Complex::conj
1297 * @tparam T Arithmetic type of a
1298 * @param a Complex number to get conjugate from
1299 * @return Complex<T>
1300 */
1301template <typename T>
1302inline Complex<T> dagger(const Complex<T> &val) {
1303 return val.conj();
1304}
1305
1306/**
1307 * @brief Return Squarenorm of Complex number
1308 * @details Wrapper around Complex::squarenorm
1309 * @tparam T Arithmetic type of val
1310 * @param val Complex number to compute squarenorm of
1311 * @return T
1312 */
1313template <typename T>
1314inline auto squarenorm(const Complex<T> &val) {
1315 return val.squarenorm();
1316}
1317
1318
1319//////////////////////////////////////////////////////////////////////////////////
1320/// Print a complex value as (re,im)
1321//////////////////////////////////////////////////////////////////////////////////
1322
1323/**
1324 * @brief Stream operator
1325 * @details Print Complex number as "A.re A.im"
1326 * @tparam T Arithmetic type of A
1327 * @param strm
1328 * @param A Complex number to print
1329 * @return std::ostream&
1330 */
1331template <typename T>
1332std::ostream &operator<<(std::ostream &strm, const Complex<T> &A) {
1333 return strm << A.real() << ' ' << A.imag();
1334}
1335
1336/////////////////////////////////////////////////////////////////////////////////
1337// Function hila::to_string
1338
1339namespace hila {
1340
1341/**
1342 * @brief Cast to different basic number type
1343 *
1344 * hila::cast_to<double>(a); - cast complex a to Complex<double>
1345 */
1346
1347template <typename Ntype, typename T>
1349 Complex<Ntype> res;
1350 res.re = static_cast<Ntype>(m.re);
1351 res.im = static_cast<Ntype>(m.im);
1352 return res;
1353}
1354
1355/**
1356 * @brief Return Complex number as std::string
1357 *
1358 * @tparam T Arithmetic type of A
1359 * @param A Complex number to convert
1360 * @param prec Precision to represent Complex number at
1361 * @param separator Separator for Complex number
1362 * @return std::string
1363 */
1364template <typename T>
1365std::string to_string(const Complex<T> &A, int prec = 8, char separator = ' ') {
1366 return A.str(prec, separator);
1367}
1368
1369/**
1370 * @brief Return well formatted Complex number as std::string
1371 * @details Output is as "(A.re, A.im)"
1372 *
1373 * @tparam T Arithmetic type of A
1374 * @param A Complex number to print
1375 * @param prec Precision to print Complex number as
1376 * @return std::string
1377 */
1378template <typename T>
1379std::string prettyprint(const Complex<T> &A, int prec = 8) {
1380 std::stringstream ss;
1381 ss.precision(prec);
1382 ss << "( " << A.real() << ", " << A.imag() << " )";
1383 return ss.str();
1384}
1385
1386
1387
1388
1389} // namespace hila
1390
1391
1392/**
1393 * @brief Imaginary type, used to represent purely imaginary numbers
1394 *
1395 * Useful for reducing multiply operations in im * complex or im * real -ops
1396 * Derived from Complex class, so generic complex ops should remain valid
1397 * Defines only operators * and /, others go via Complex class
1398 *
1399 * Note: Imaginary_t should NOT be used in Field variables
1400 *
1401 * @tparam T type of imaginary (float/double)
1402 */
1403template <typename T>
1404class Imaginary_t : public Complex<T> {
1405 public:
1406 constexpr Imaginary_t() = default;
1407 ~Imaginary_t() = default;
1408 constexpr Imaginary_t(const Imaginary_t &i) = default;
1409
1410 // construct from scalar
1411 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value, int> = 0>
1412 explicit constexpr Imaginary_t(const A v) : Complex<T>(0, v) {}
1413
1414 constexpr Imaginary_t operator-() const {
1415 return Imaginary_t<T>(-this->im);
1416 }
1417
1418 constexpr Imaginary_t operator+() const {
1419 return *this;
1420 }
1421};
1422
1423///////////////////////////////////////////////////////////////////////////////////////
1424/// @brief Imaginary unit I - global variable
1425///
1426///
1427// separate definition for GPU kernel code and host (and generic) code
1428
1429#if defined(_GPU_DEVICE_COMPILE_)
1430__device__ constexpr Imaginary_t<double> I(1.0);
1431#else
1432constexpr Imaginary_t<double> I(1.0);
1433#endif
1434
1435
1436///////////////////////////////////////////////////////////////////////////////////////
1437
1438// Imaginary * object containing complex
1439template <typename T, std::enable_if_t<hila::contains_complex<T>::value, int> = 0>
1440inline auto operator*(const Imaginary_t<double> &i, const T &c) {
1442 T res;
1443 constexpr int n_cmplx = sizeof(T) / sizeof(Complex<hila::arithmetic_type<T>>);
1444 for (int k = 0; k < n_cmplx; k++) {
1445 ca = hila::get_complex_in_var(c, k);
1446 cb.re = -ca.im * i.imag();
1447 cb.im = ca.re * i.imag();
1448 hila::set_complex_in_var(res, k, cb);
1449 }
1450 return res;
1451}
1452
1453// object containing complex * imaginary
1454template <typename T, std::enable_if_t<hila::contains_complex<T>::value, int> = 0>
1455inline auto operator*(const T &c, const Imaginary_t<double> &i) {
1456 return i * c;
1457}
1458
1459// Imag * scalar, returns imag
1460// note: using std::is_arithmetic, not done for vector types
1461template <typename A, typename T, std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>
1462inline Imaginary_t<A> operator*(Imaginary_t<A> i, const T &c) {
1463 i.imag() *= c;
1464 return i;
1465}
1466
1467// scalar * imag, returns imag
1468template <typename A, typename T, std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>
1469inline Imaginary_t<A> operator*(const T &c, Imaginary_t<A> i) {
1470 i.imag() *= c;
1471 return i;
1472}
1473
1474// Imaginary * Imaginary, returns scalar
1475template <typename A, typename B>
1476inline auto operator*(const Imaginary_t<A> &a, const Imaginary_t<B> &b) {
1477 return -(a.imag() * b.imag());
1478}
1479
1480////////////////////////////
1481
1482
1483// @brief Imaginary / real value, returning imaginary
1484// Note: for vectorized types the generic version is used
1485template <typename T, typename A, std::enable_if_t<std::is_arithmetic<A>::value, int> = 0>
1486inline auto operator/(Imaginary_t<T> i, const A &a) {
1487 i.imag() /= a;
1488 return i;
1489}
1490
1491// @brief Imaginary / Imaginary, returning real
1492template <typename A, typename B>
1493inline auto operator/(const Imaginary_t<A> &a, const Imaginary_t<B> &b) {
1494 return a.imag() / b.imag();
1495}
1496
1497
1498///////////////////////////////////////////////////////////////////////////////
1499// Set of complex functions
1500///////////////////////////////////////////////////////////////////////////////
1501/**
1502 * @name Mathematical functions
1503 * @return Complex T
1504 * @memberof Complex
1505 *
1506 * @{
1507 */
1508/// \f$\exp(z)\f$
1509template <typename T>
1510inline Complex<T> exp(const Complex<T> z) {
1511 return exp(z.re) * Complex<T>(cos(z.im), sin(z.im));
1512}
1513
1514/// \f$\exp(i\cdot x)\f$
1515template <typename T, std::enable_if_t<hila::is_arithmetic<T>::value, int> = 0>
1516inline Complex<T> expi(T a) {
1517 return Complex<T>(cos(a), sin(a));
1518}
1519
1520/// \f$\exp(\Im(z))\f$
1521template <typename T>
1523 return expi(im.imag());
1524}
1525
1526/// \f$\log{z}\f$
1527template <typename T>
1529 return Complex<T>(static_cast<T>(0.5) * log(z.squarenorm()), z.arg());
1530}
1531
1532// sqrt(z) branch cut at -x axis
1533
1534/// \f$\sqrt{z}\f$
1535template <typename T>
1537 T r = z.squarenorm();
1538 T a = z.arg();
1539 return pow(r, 0.25) * expi(0.5 * a);
1540}
1541
1542/// \f$\sqrt[3]{z}\f$
1543template <typename T>
1545 T r = z.squarenorm();
1546 T a = z.arg();
1547 return pow(r, (1.0 / 6.0)) * expi((1.0 / 3.0) * a);
1548}
1549
1550/// pow(z.p) = \f$z^p\f$ = \f$exp(p*log(z))\f$
1551template <typename A, typename B>
1552inline auto pow(Complex<A> z, Complex<B> p) {
1553 return exp(p * log(z));
1554}
1555
1556/// pow(z.p) with scalar power
1557template <typename T, typename S, std::enable_if_t<hila::is_arithmetic<S>::value, int> = 0>
1558inline Complex<T> pow(Complex<T> z, S p) {
1559 return exp(p * log(z));
1560}
1561
1562/// pow(z.p) with scalar base
1563template <typename T, typename S, std::enable_if_t<hila::is_arithmetic<S>::value, int> = 0>
1564inline Complex<T> pow(S z, Complex<T> p) {
1565 return exp(p * log(z));
1566}
1567
1568/// \f{align}{sin(z) &= \sin(\Re(z) + i \Im(z)) \\
1569/// &= \sin(\Re(z))\cos(i \Im(z)) + \cos(\Re(z))\sin(i \Im(z)) \\
1570/// &= \sin(\Re(z)) \cosh(\Im(z)) + i \cos(\Re(z)) \sinh(\Im(z)) \f}
1571template <typename T>
1573 return Complex<T>(sin(z.re) * cosh(z.im), cos(z.re) * sinh(z.im));
1574}
1575
1576/// \f{align}{\cos(z) &= \cos(\Re{z})\cos(i \Im(z)) - \sin(\Re{z})\sin(i \Im{z}) \\
1577/// &= \cos(\Re{z})\cosh(\Im(z)) - i \sin(\Re{z})\sinh(\Im(z))\f}
1578template <typename T>
1580 return Complex<T>(cos(z.re) * cosh(z.im), -sin(z.re) * sinh(z.im));
1581}
1582
1583/// \f$\tan(z) = \frac{\sin(z)}{\cos(z)}\f$ - rely on optimizer to simplify
1584template <typename T>
1586 return sin(z) / cos(z);
1587}
1588
1589/// sinh(z) = sinh(re)cosh(i im) + cosh(re)sinh(i im)
1590/// = sinh(re)cos(im) + i cosh(re)sin(im)
1591template <typename T>
1593 return Complex<T>(sinh(z.re) * cos(z.im), cosh(z.re) * sin(z.im));
1594}
1595
1596/// cosh(z) = cosh(re)cosh(i im) - sinh(re)sinh(i im)
1597/// = cosh(re)cos(im) - i sinh(re)sin(im)
1598template <typename T>
1600 return Complex<T>(cosh(z.re) * cos(z.im), sinh(z.re) * sin(z.im));
1601}
1602
1603/// \f$\tanh(z)\f$
1604template <typename T>
1606 return sinh(z) / cosh(z);
1607}
1608
1609/// \f$\arctan(z)\f$
1610template <typename T>
1612 return -0.5 * I * (log((I - z) / (I + z)));
1613}
1614
1615/// \f$\arcsin(z)\f$
1616template <typename T>
1618 return -I * (log(I * z + sqrt(1 - z * z)));
1619}
1620
1621/// \f$\arccos(z)\f$
1622template <typename T>
1624 return -I * (log(z + I * (sqrt(1 - z * z))));
1625}
1626
1627/// \f$\text{artanh}(z)\f$
1628template <typename T>
1630 return 0.5 * log((1 + z) / (1 - z));
1631}
1632
1633/// \f$\text{arsinh}(z)\f$
1634template <typename T>
1636 return log(z + sqrt(1 + z * z));
1637}
1638
1639/// \f$\text{arcosh}(z)\f$
1640template <typename T>
1642 return log(z + sqrt(z * z - 1));
1643}
1644
1645/** @} */
1646
1647
1648//////////////////////////////////////////////////////////////////////////////////
1649/// Operators to implement imaginary unit 1_i, enablig expressions 3 + 2_i etc.
1650/// This is defined as an user-defined literal, which requires an underscore.
1651////////////////////////////////////////////////////////////////////////////////
1652
1653constexpr Imaginary_t<double> operator""_i(long double a) {
1654 return Imaginary_t<double>{a};
1655}
1656
1657constexpr Imaginary_t<double> operator""_i(unsigned long long a) {
1658 return Imaginary_t<double>(a);
1659}
1660
1661
1662#endif
Complex definition.
Definition cmplx.h:58
Complex< T > dagger() const
Compute dagger of Complex number.
Definition cmplx.h:294
auto operator*(const A &a, const Complex< T > &c)
Multiplication operator Scalar * Complex.
Definition cmplx.h:1050
Complex< T > asin(Complex< T > z)
Definition cmplx.h:1617
auto operator+(const A &a, const Complex< T > &c)
Addition operator Scalar + Complex.
Definition cmplx.h:929
Complex< T > & operator*=(const Complex< A > &lhs) &
Complex multiply assignment operator.
Definition cmplx.h:477
Complex< T > cbrt(Complex< T > z)
Definition cmplx.h:1544
bool operator!=(const A a, const Complex< B > &b)
Compare non-equality of Scalar and Complex number.
Definition cmplx.h:1249
Complex< T > exp(const Imaginary_t< T > im)
Definition cmplx.h:1522
Complex< Tr > operator-(const Complex< T1 > &a, const Complex< T2 > &b)
Subtraction operator Complex - Complex.
Definition cmplx.h:952
Complex< T > & operator*=(const A a) &
Real multiply assign operator.
Definition cmplx.h:497
auto operator-(const Complex< T > &c, const A &a)
Subtraction operator Complex - Scalar.
Definition cmplx.h:970
bool operator==(const Complex< A > &a, const Complex< B > &b)
Compare equality of two complex numbersTwo numbers are equal, if the real and imaginary components ar...
Definition cmplx.h:1163
Complex< T > conj() const
Compute conjugate of Complex number.
Definition cmplx.h:283
auto operator*(const Complex< T > &c, const A &a)
Multiplication operator Complex * Scalar.
Definition cmplx.h:1030
Complex< T > cosh(Complex< T > z)
Definition cmplx.h:1599
Complex< T > acosh(Complex< T > z)
Definition cmplx.h:1641
bool operator==(const A a, const Complex< B > &b)
Compare equality of Scalar and Complex.
Definition cmplx.h:1201
T abs() const
Compute absolute value of Complex number.
Definition cmplx.h:260
auto operator/(const Complex< T > &c, const A &a)
Division operator Complex / Scalar.
Definition cmplx.h:1099
Complex< T > polar(const T r, const T theta)
Stores and returns Complex number given in polar coordinates.
Definition cmplx.h:314
T real() const
Real part of Complex number.
Definition cmplx.h:157
Complex< Tr > operator*(const Complex< T1 > &a, const Complex< T2 > &b)
Multiplication operator Complex * Complex.
Definition cmplx.h:1010
Complex< T > & operator++()
Single increment operator.
Definition cmplx.h:575
T arg() const
Compute argument of Complex number.
Definition cmplx.h:271
Complex< T > pow(Complex< T > z, S p)
pow(z.p) with scalar power
Definition cmplx.h:1558
Complex< Tr > operator/(const Complex< T1 > &a, const Complex< T2 > &b)
Division operator Complex / Complex.
Definition cmplx.h:1078
Complex< Tr > operator+(const Complex< T1 > &a, const Complex< T2 > &b)
Addition operator Complex + Complex.
Definition cmplx.h:892
bool operator!=(const Complex< A > &a, const Complex< B > &b)
Compare non-equality of two complex numbers.
Definition cmplx.h:1217
Complex< T > operator++(int)
Multiple increment operator.
Definition cmplx.h:609
Complex< T > & operator+=(const Complex< A > &lhs)
Complex addition assignment operator.
Definition cmplx.h:394
T & imag()
Imaginary part of Complex number.
Definition cmplx.h:185
Complex< T > tanh(Complex< T > z)
Definition cmplx.h:1605
Complex< T > sin(Complex< T > z)
Definition cmplx.h:1572
Complex< T > atan(Complex< T > z)
Definition cmplx.h:1611
Complex< T > & operator+=(const A &a) &
Real addition assignment operator.
Definition cmplx.h:414
Complex< T > tan(Complex< T > z)
- rely on optimizer to simplify
Definition cmplx.h:1585
Complex< T > & operator/=(const Complex< A > &lhs) &
Complex divide assignment operator.
Definition cmplx.h:533
Complex< T > operator--(int)
Multiple decrement operator.
Definition cmplx.h:628
T squarenorm() const
Compute square norm of Complex number.
Definition cmplx.h:249
Complex< T > log(Complex< T > z)
Definition cmplx.h:1528
Complex< T > & gaussian_random(double width=1.0)
Produces complex gaussian random values.
Definition cmplx.h:339
Complex< T > & operator/=(const A &a) &
Real divide assignment operator.
Definition cmplx.h:554
Complex< T > & operator--()
Single decrement operator.
Definition cmplx.h:592
T & real()
Real part of Complex number.
Definition cmplx.h:176
Complex< T > operator+() const
Unary + operator.
Definition cmplx.h:357
T imag() const
Imaginary part of Complex number.
Definition cmplx.h:166
Complex< T > mul_add(const Complex< T > &a, const Complex< T > &b, const Complex< T > &c)
Multiply add with Complex numbers.
Definition cmplx.h:1140
Complex< T > & random()
Assign random values to Complex real and imaginary part.
Definition cmplx.h:326
Complex< T > pow(S z, Complex< T > p)
pow(z.p) with scalar base
Definition cmplx.h:1564
Complex< T > operator-() const
Unary - operator.
Definition cmplx.h:372
Complex< T > mul_conj(const Complex< A > &b) const
Multiply conjugate method.
Definition cmplx.h:670
Complex< T > expi(T a)
Definition cmplx.h:1516
auto operator+(const Complex< T > &c, const A &a)
Addition operator Complex + Scalar.
Definition cmplx.h:911
bool operator==(const Complex< A > &a, const B b)
Compare equality of Complex and scalar.
Definition cmplx.h:1182
Complex< T > & operator-=(const Complex< A > &lhs) &
Complex subtraction assignment operator.
Definition cmplx.h:434
Complex< T > conj_mul(const Complex< A > &b) const
Conjugate multiply method.
Definition cmplx.h:654
Complex< T > cos(Complex< T > z)
Definition cmplx.h:1579
Complex< T > sinh(Complex< T > z)
Definition cmplx.h:1592
auto pow(Complex< A > z, Complex< B > p)
pow(z.p) = =
Definition cmplx.h:1552
Complex< T > atanh(Complex< T > z)
Definition cmplx.h:1629
Complex< T > acos(Complex< T > z)
Definition cmplx.h:1623
bool operator!=(const Complex< A > &a, const B b)
Compare non-equality of Complex number and Scalar.
Definition cmplx.h:1233
auto operator-(const A &a, const Complex< T > &c)
Subtraction operator Scalar - Complex.
Definition cmplx.h:988
auto operator/(const A &a, const Complex< T > &c)
Division operator Scalar / Complex.
Definition cmplx.h:1120
Complex< T > & operator=(const Complex< T > &s) &=default
Assignment operator.
Complex< T > sqrt(Complex< T > z)
Definition cmplx.h:1536
Complex< T > asinh(Complex< T > z)
Definition cmplx.h:1635
Complex< T > exp(const Complex< T > z)
Definition cmplx.h:1510
Imaginary type, used to represent purely imaginary numbers.
Definition cmplx.h:1404
T imag(const Complex< T > &a)
Retrun imaginary value of Complex number.
Definition cmplx.h:846
Complex< T > conj(const Complex< T > &val)
Return conjugate of Complex number.
Definition cmplx.h:1290
auto squarenorm(const Complex< T > &val)
Return Squarenorm of Complex number.
Definition cmplx.h:1314
Complex< T > dagger(const Complex< T > &val)
Return dagger of Complex number.
Definition cmplx.h:1302
Complex< T > polar(T r, T arg)
Return complex number given by polar representation.
Definition cmplx.h:859
T real(const Complex< T > &a)
Return real value of Complex number.
Definition cmplx.h:834
constexpr Imaginary_t< double > I(1.0)
Imaginary unit I - global variable.
T abs(const Complex< T > &a)
Return absolute value of Complex number.
Definition cmplx.h:1266
std::ostream & operator<<(std::ostream &strm, const Complex< T > &A)
Print a complex value as (re,im)
Definition cmplx.h:1332
T arg(const Complex< T > &a)
Return argument of Complex number.
Definition cmplx.h:1278
This file defines all includes for HILA.
This files containts definitions for the extended precision class that allows for high precision redu...
Implement hila::swap for gauge fields.
Definition array.h:982
double random()
Real valued uniform random number generator.
Definition hila_gpu.cpp:121
Array< n, m, Ntype > cast_to(const Array< n, m, T > &mat)
Array casting operation.
Definition array.h:1292
double gaussrand2(double &out2)
hila::gaussrand2 returns 2 Gaussian distributed random numbers with variance .
Definition random.cpp:188
decltype(std::declval< A >()+std::declval< B >()) type_plus
Definition type_tools.h:103
std::string to_string(const Array< n, m, T > &A, int prec=8, char separator=' ')
Converts Array object to string.
Definition array.h:996
hila::is_complex_or_arithmetic<T>::value
Definition cmplx.h:708