23template <typename T, std::enable_if_t<hila::is_arithmetic<T>::value,
int> = 0>
24inline T mul_add(T a, T b, T c) {
28template <typename T, std::enable_if_t<hila::is_arithmetic<T>::value,
int> = 0>
29inline T mul_sub(T a, T b, T c) {
33template <typename T, std::enable_if_t<hila::is_arithmetic<T>::value,
int> = 0>
34inline T nmul_add(T a, T b, T c) {
53 "Complex can be used only with floating point type");
57 using base_type = hila::arithmetic_type<T>;
58 using argument_type = T;
120 template <typename S, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
121 explicit constexpr Complex<T>(
const S val) : re(val), im(0) {}
135 constexpr Complex<T>(
const std::nullptr_t n) : re(0), im(0) {}
153 template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0,
154 std::enable_if_t<hila::is_arithmetic<B>::value,
int> = 0>
155 explicit constexpr Complex<T>(
const A &a,
const B &b) : re(a), im(b) {}
185 inline T &
real() const_function {
194 inline T &
imag() const_function {
199 template <
typename A>
200#pragma hila loop_function
245 template <
typename A>
267 template <typename S, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
282 return re * re + im * im;
304 return atan2(im, re);
425 template <
typename A>
446 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
466 template <
typename A>
486 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
520 template <
typename A>
522 T r = mul_sub(re, lhs.re, im * lhs.im);
523 im =
mul_add(im, lhs.re, re * lhs.im);
540 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
576 template <
typename A>
579 T r =
mul_add(re, lhs.re, im * lhs.im) / n;
580 im = mul_sub(im, lhs.re, re * lhs.im) / n;
597 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
678 std::string str(
int prec = 8,
char separator =
' ')
const {
679 std::stringstream ss;
681 ss << re << separator << im;
697 template <
typename A>
699 return Complex<T>(re * b.re + im * b.im, re * b.im - im * b.re);
713 template <
typename A>
715 return Complex<T>(re * b.re + im * b.im, im * b.re - re * b.im);
719 template <
typename Ntype>
742struct is_complex : std::integral_constant<bool, false> {};
745struct is_complex<
Complex<T>> : std::integral_constant<bool, true> {};
750 : std::integral_constant<bool, hila::is_arithmetic<T>::value || hila::is_complex<T>::value> {};
756template <
typename T,
typename Enable =
void>
757struct contains_complex : std::integral_constant<bool, false> {};
760struct contains_complex<T, typename
std::enable_if_t<hila::is_field_class_type<T>::value>>
761 : std::integral_constant<bool,
762 hila::contains_type<T, Complex<hila::arithmetic_type<T>>>::value> {};
768template <
typename T,
typename Enable =
void,
typename N = hila::arithmetic_type<T>>
769struct complex_or_arithmetic_type_struct {
774struct complex_or_arithmetic_type_struct<
775 T, typename
std::enable_if_t<hila::contains_complex<T>::value>> {
780using number_type =
typename complex_or_arithmetic_type_struct<T>::type;
787template <
typename T1,
typename T2,
typename Enable =
void>
793template <
typename T1,
typename T2>
794struct ntype_op_s<T1, T2,
795 typename
std::enable_if_t<(hila::contains_complex<T1>::value ||
796 hila::contains_complex<T2>::value)>> {
800template <
typename T1,
typename T2>
801using ntype_op =
typename ntype_op_s<T1, T2>::type;
809template <
typename A,
typename B,
typename Enable =
void>
810struct complex_x_scalar_s {
814template <
typename A,
typename B>
815struct complex_x_scalar_s<
816 A, B, typename
std::enable_if_t<std::is_assignable<A &, hila::type_plus<A, B>>::value>> {
820template <
typename A,
typename B>
821using complex_x_scalar_type =
typename complex_x_scalar_s<A, B>::type;
851template <typename T, std::enable_if_t<hila::contains_complex<T>::value,
int> = 0>
856template <typename T, std::enable_if_t<hila::contains_complex<T>::value,
int> = 0>
857inline void set_complex_in_var(T &var,
int i,
const Complex<hila::arithmetic_type<T>> &val) {
864template <typename T, std::enable_if_t<hila::is_floating_point<T>::value,
int> = 0>
915#pragma hila loop_function
935template <
typename T1,
typename T2,
typename Tr = hila::type_plus<T1, T2>>
954template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
956 return hila::complex_x_scalar_type<T, A>(c.re + a, c.im);
972template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
974 return hila::complex_x_scalar_type<T, A>(c.re + a, c.im);
995template <
typename T1,
typename T2,
typename Tr = hila::type_plus<T1, T2>>
1013template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1015 return hila::complex_x_scalar_type<T, A>(c.re - a, c.im);
1031template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1033 return hila::complex_x_scalar_type<T, A>(a - c.re, -c.im);
1053template <
typename T1,
typename T2,
typename Tr = hila::type_mul<T1, T2>>
1055 return Complex<Tr>(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re);
1073template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1075 return hila::complex_x_scalar_type<T, A>(c.re * a, c.im * a);
1093template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1095 return hila::complex_x_scalar_type<T, A>(a * c.re, a * c.im);
1121template <
typename T1,
typename T2,
typename Tr = hila::type_mul<T1, T2>>
1124 return Complex<Tr>((a.re * b.re + a.im * b.im) * n, (a.im * b.re - a.re * b.im) * n);
1142template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1144 return hila::complex_x_scalar_type<T, A>(c.re / a, c.im / a);
1163template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1166 return hila::complex_x_scalar_type<T, A>((a * c.re) / n, -(a * c.im) / n);
1183template <
typename T>
1187 T t1 = mul_add(a.re, b.re, c.re);
1188 T t2 = mul_add(a.re, b.im, c.im);
1189 r.re = mul_add(a.im, b.im, t1);
1190 r.im = mul_add(a.im, b.re, t2);
1206template <
typename A,
typename B>
1208 return (a.re == b.re && a.im == b.im);
1225template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<B>::value,
int> = 0>
1227 return (a.re == b && a.im == 0);
1244template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1260template <
typename A,
typename B>
1262 return (a.re != b.re || a.im != b.im);
1276template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<B>::value,
int> = 0>
1278 return (a.re != b || a.im != 0);
1292template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1304template <typename Ntype, typename T, std::enable_if_t<hila::is_arithmetic<Ntype>::value,
int> = 0>
1321template <
typename T>
1333template <
typename T>
1345template <
typename T>
1357template <
typename T>
1369template <
typename T>
1387template <
typename T>
1389 return strm << A.
real() <<
' ' << A.
imag();
1405template <
typename T>
1407 return A.str(prec, separator);
1419template <
typename T>
1421 std::stringstream ss;
1423 ss <<
"( " << A.
real() <<
", " << A.
imag() <<
" )";
1441template <
typename T>
1449#pragma hila loop_function
1450 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1468#if defined(CUDA) || defined(HIP)
1480template <typename A, typename T, std::enable_if_t<hila::contains_complex<T>::value,
int> = 0>
1484 constexpr int n_cmplx =
sizeof(T) /
sizeof(
Complex<hila::arithmetic_type<T>>);
1485 for (
int k = 0; k < n_cmplx; k++) {
1486 ca = hila::get_complex_in_var(c, k);
1487 cb.re = -ca.im * i.
imag();
1488 cb.im = ca.re * i.
imag();
1489 hila::set_complex_in_var(res, k, cb);
1495template <typename A, typename T, std::enable_if_t<hila::contains_complex<T>::value,
int> = 0>
1503template <typename A, typename T, std::enable_if_t<std::is_arithmetic<T>::value,
int> = 0>
1510template <typename A, typename T, std::enable_if_t<std::is_arithmetic<T>::value,
int> = 0>
1517template <
typename A,
typename B>
1527template <typename T, typename A, std::enable_if_t<std::is_arithmetic<A>::value,
int> = 0>
1534template <
typename A,
typename B>
1551template <
typename T>
1553 return exp(z.re) *
Complex<T>(cos(z.im), sin(z.im));
1557template <typename T, std::enable_if_t<hila::is_arithmetic<T>::value,
int> = 0>
1563template <
typename T>
1565 return expi(im.
imag());
1569template <
typename T>
1577template <
typename T>
1581 return pow(r, 0.25) * expi(0.5 * a);
1585template <
typename T>
1589 return pow(r, (1.0 / 6.0)) * expi((1.0 / 3.0) * a);
1593template <
typename A,
typename B>
1595 return exp(p * log(z));
1599template <typename T, typename S, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
1601 return exp(p * log(z));
1605template <typename T, typename S, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
1607 return exp(p * log(z));
1613template <
typename T>
1615 return Complex<T>(sin(z.re) * cosh(z.im), cos(z.re) * sinh(z.im));
1620template <
typename T>
1622 return Complex<T>(cos(z.re) * cosh(z.im), -sin(z.re) * sinh(z.im));
1626template <
typename T>
1628 return sin(z) / cos(z);
1633template <
typename T>
1635 return Complex<T>(sinh(z.re) * cos(z.im), cosh(z.re) * sin(z.im));
1640template <
typename T>
1642 return Complex<T>(cosh(z.re) * cos(z.im), sinh(z.re) * sin(z.im));
1646template <
typename T>
1648 return sinh(z) / cosh(z);
1652template <
typename T>
1654 return -0.5 *
I * (log((
I - z) / (
I + z)));
1658template <
typename T>
1660 return -
I * (log(
I * z + sqrt(1 - z * z)));
1664template <
typename T>
1666 return -
I * (log(z +
I * (sqrt(1 - z * z))));
1670template <
typename T>
1672 return 0.5 * log((1 + z) / (1 - z));
1676template <
typename T>
1678 return log(z + sqrt(1 + z * z));
1682template <
typename T>
1684 return log(z + sqrt(z * z - 1));
Complex< T > dagger() const
Compute dagger of Complex number.
auto operator*(const A &a, const Complex< T > &c)
Multiplication operator Scalar * Complex.
Complex< T > asin(Complex< T > z)
Complex< T > & operator/=(const Complex< A > &lhs)
Complex divide assignment operator.
auto operator+(const A &a, const Complex< T > &c)
Addition operator Scalar + Complex.
Complex< T > & operator-=(const Complex< A > &lhs)
Complex subtraction assignment operator.
Complex< T > cbrt(Complex< T > z)
bool operator!=(const A a, const Complex< B > &b)
Compare non-equality of Scalar and Complex number.
Complex< T > exp(const Imaginary_t< T > im)
Complex< Tr > operator-(const Complex< T1 > &a, const Complex< T2 > &b)
Subtraction operator Complex - Complex.
auto operator-(const Complex< T > &c, const A &a)
Subtraction operator Complex - Scalar.
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...
Complex< T > & operator*=(const A a)
Real multiply assign operator.
Complex< T > conj() const
Compute conjugate of Complex number.
auto operator*(const Complex< T > &c, const A &a)
Multiplication operator Complex * Scalar.
Complex< T > cosh(Complex< T > z)
Complex< T > acosh(Complex< T > z)
bool operator==(const A a, const Complex< B > &b)
Compare equality of Scalar and Complex.
T abs() const
Compute absolute value of Complex number.
auto operator/(const Complex< T > &c, const A &a)
Division operator Complex / Scalar.
Complex< T > polar(const T r, const T theta)
Stores and returns Complex number given in polar coordinates.
T real() const
Real part of Complex number.
Complex< Tr > operator*(const Complex< T1 > &a, const Complex< T2 > &b)
Multiplication operator Complex * Complex.
Complex< T > & operator++()
Single increment operator.
T arg() const
Compute argument of Complex number.
Complex< T > pow(Complex< T > z, S p)
pow(z.p) with scalar power
Complex< Tr > operator/(const Complex< T1 > &a, const Complex< T2 > &b)
Division operator Complex / Complex.
Complex< Tr > operator+(const Complex< T1 > &a, const Complex< T2 > &b)
Addition operator Complex + Complex.
bool operator!=(const Complex< A > &a, const Complex< B > &b)
Compare non-equality of two complex numbers.
Complex< T > operator++(int)
Multiple increment operator.
Complex< T > & operator+=(const Complex< A > &lhs)
Complex addition assignment operator.
T & imag()
Imaginary part of Complex number.
Complex< T > tanh(Complex< T > z)
Complex< T > sin(Complex< T > z)
Complex< T > atan(Complex< T > z)
Complex< T > tan(Complex< T > z)
- rely on optimizer to simplify
Complex< T > operator--(int)
Multiple decrement operator.
T squarenorm() const
Compute square norm of Complex number.
Complex< T > log(Complex< T > z)
Complex< T > & operator*=(const Complex< A > &lhs)
Complex multiply assignment operator.
Complex< T > & gaussian_random(double width=1.0)
Produces complex gaussian random values.
Complex< T > & operator/=(const A &a)
Real divide assignment operator.
Complex< T > & operator--()
Single decrement operator.
T & real()
Real part of Complex number.
Complex< T > operator+() const
Unary + operator.
Complex< T > & operator=(S s)
Assignment from real.
Complex< T > & operator-=(const A &a)
Real subtraction assignment operator Subtract assign only to real part of Complex number.
T imag() const
Imaginary part of Complex number.
Complex< T > mul_add(const Complex< T > &a, const Complex< T > &b, const Complex< T > &c)
Multiply add with Complex numbers.
Complex< T > & random()
Assign random values to Complex real and imaginary part.
Complex< T > pow(S z, Complex< T > p)
pow(z.p) with scalar base
Complex< T > operator-() const
Unary - operator.
Complex< T > mul_conj(const Complex< A > &b) const
Multiply conjugate method.
auto operator+(const Complex< T > &c, const A &a)
Addition operator Complex + Scalar.
bool operator==(const Complex< A > &a, const B b)
Compare equality of Complex and scalar.
Complex< T > & operator=(const Complex< A > &s)
Assignment from complex.
Complex< T > & operator=(const Complex< T > &s)=default
Assignment operator.
Complex< T > conj_mul(const Complex< A > &b) const
Conjugate multiply method.
Complex< T > cos(Complex< T > z)
Complex< T > sinh(Complex< T > z)
auto pow(Complex< A > z, Complex< B > p)
pow(z.p) = =
Complex< T > atanh(Complex< T > z)
Complex< T > acos(Complex< T > z)
bool operator!=(const Complex< A > &a, const B b)
Compare non-equality of Complex number and Scalar.
auto operator-(const A &a, const Complex< T > &c)
Subtraction operator Scalar - Complex.
auto operator/(const A &a, const Complex< T > &c)
Division operator Scalar / Complex.
Complex< T > sqrt(Complex< T > z)
Complex< T > asinh(Complex< T > z)
Complex< T > exp(const Complex< T > z)
Complex< T > & operator+=(const A &a)
Real addition assignment operator.
Imaginary type, used to represent purely imaginary numbers.
T imag(const Complex< T > &a)
Retrun imaginary value of Complex number.
Complex< T > conj(const Complex< T > &val)
Return conjugate of Complex number.
auto squarenorm(const Complex< T > &val)
Return Squarenorm of Complex number.
Complex< T > dagger(const Complex< T > &val)
Return dagger of Complex number.
Complex< T > polar(T r, T arg)
Return complex number given by polar representation.
T real(const Complex< T > &a)
Return real value of Complex number.
constexpr Imaginary_t< double > I(1.0)
Imaginary unit I - global variable.
T abs(const Complex< T > &a)
Return absolute value of Complex number.
std::ostream & operator<<(std::ostream &strm, const Complex< T > &A)
Print a complex value as (re,im)
T arg(const Complex< T > &a)
Return argument of Complex number.
This file defines all includes for HILA.
Invert diagonal + const. matrix using Sherman-Morrison formula.
double random()
Real valued uniform random number generator.
double gaussrand2(double &out2)
hila::gaussrand2 returns 2 Gaussian distributed random numbers with variance .
decltype(std::declval< A >()+std::declval< B >()) type_plus
std::string to_string(const Array< n, m, T > &A, int prec=8, char separator=' ')
Converts Array object to string.
hila::is_complex_or_arithmetic<T>::value