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) {
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) {
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) {
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");
65 using base_type = hila::arithmetic_type<T>;
66 using argument_type = T;
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)) {}
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) {}
137 constexpr Complex<T>(
const std::nullptr_t n) : re(0), im(0) {}
141 typename A,
typename B,
142 std::enable_if_t<hila::is_assignable<T &, A>::value && hila::is_assignable<T &, B>::value,
144 explicit constexpr Complex<T>(
const A &a,
const B &b) : re(a), im(b) {}
176 inline T &
real() const_function {
185 inline T &
imag() const_function {
224 template <typename A, std::enable_if_t<hila::is_assignable<T&,A>::value,
int> = 0>
231 template <typename S, std::enable_if_t<hila::is_assignable<T&,S>::value,
int> = 0>
239 template <
typename S>
250 return re * re + im * im;
272 return atan2(im, re);
393 template <
typename A>
413 template <typename A, std::enable_if_t<hila::is_assignable<T &, A>::value,
int> = 0>
433 template <
typename A>
436 "Cannot assign ExtendedPrecision value to non-EP variable without casting");
442 template <typename A, std::enable_if_t<hila::is_arithmetic_or_extended<A>::value,
int> = 0>
445 "Cannot assign ExtendedPrecision value to non-EP variable without casting");
476 template <
typename A>
478 T r = mul_sub(re, lhs.re, im * lhs.im);
479 im =
mul_add(im, lhs.re, re * lhs.im);
496 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
532 template <
typename A>
534 T n = lhs.squarenorm();
535 T r =
mul_add(re, lhs.re, im * lhs.im) / n;
536 im = mul_sub(im, lhs.re, re * lhs.im) / n;
553 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
634 std::string str(
int prec = 8,
char separator =
' ')
const {
635 std::stringstream ss;
637 ss << re << separator << im;
653 template <
typename A>
655 return Complex<T>(re * b.re + im * b.im, re * b.im - im * b.re);
669 template <
typename A>
671 return Complex<T>(re * b.re + im * b.im, im * b.re - re * b.im);
697struct is_complex : std::integral_constant<bool, false> {};
700struct is_complex<
Complex<T>> : std::integral_constant<bool, true> {};
703struct is_complex<
Imaginary_t<T>> : std::integral_constant<bool, true> {};
708 : std::integral_constant<bool, hila::is_arithmetic<T>::value || hila::is_complex<T>::value> {};
714template <
typename T,
typename Enable =
void>
715struct contains_complex : std::integral_constant<bool, false> {};
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> {};
726template <
typename T,
typename Enable =
void,
typename N = hila::arithmetic_type<T>>
727struct complex_or_arithmetic_type_struct {
732struct complex_or_arithmetic_type_struct<
733 T, typename std::enable_if_t<hila::contains_complex<T>::value>> {
738using number_type =
typename complex_or_arithmetic_type_struct<T>::type;
745template <
typename T1,
typename T2,
typename Enable =
void>
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)>> {
758template <
typename T1,
typename T2>
759using ntype_op =
typename ntype_op_s<T1, T2>::type;
767template <
typename A,
typename B,
typename Enable =
void>
768struct complex_x_scalar_s {
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>> {
778template <
typename A,
typename B>
779using complex_x_scalar_type =
typename complex_x_scalar_s<A, B>::type;
807template <typename T, std::enable_if_t<hila::contains_complex<T>::value,
int> = 0>
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) {
820template <typename T, std::enable_if_t<hila::is_floating_point<T>::value,
int> = 0>
871#pragma hila loop_function
891template <
typename T1,
typename T2,
typename Tr = hila::type_plus<T1, T2>>
910template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
912 return hila::complex_x_scalar_type<T, A>(c.re + a, c.im);
928template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
930 return hila::complex_x_scalar_type<T, A>(c.re + a, c.im);
951template <
typename T1,
typename T2,
typename Tr = hila::type_plus<T1, T2>>
969template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
971 return hila::complex_x_scalar_type<T, A>(c.re - a, c.im);
987template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
989 return hila::complex_x_scalar_type<T, A>(a - c.re, -c.im);
1009template <
typename T1,
typename T2,
typename Tr = hila::type_mul<T1, T2>>
1011 return Complex<Tr>(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re);
1029template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1031 return hila::complex_x_scalar_type<T, A>(c.re * a, c.im * a);
1049template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1051 return hila::complex_x_scalar_type<T, A>(a * c.re, a * c.im);
1077template <
typename T1,
typename T2,
typename Tr = hila::type_mul<T1, T2>>
1080 return Complex<Tr>((a.re * b.re + a.im * b.im) * n, (a.im * b.re - a.re * b.im) * n);
1098template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1100 return hila::complex_x_scalar_type<T, A>(c.re / a, c.im / a);
1119template <typename T, typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1122 return hila::complex_x_scalar_type<T, A>((a * c.re) / n, -(a * c.im) / n);
1139template <
typename T>
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);
1146 r.im = mul_add(a.im, b.re, t2);
1162template <
typename A,
typename B>
1164 return (a.re == b.re && a.im == b.im);
1181template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<B>::value,
int> = 0>
1183 return (a.re == b && a.im == 0);
1200template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1216template <
typename A,
typename B>
1218 return (a.re != b.re || a.im != b.im);
1232template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<B>::value,
int> = 0>
1234 return (a.re != b || a.im != 0);
1248template <typename A, typename B, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1265template <
typename T>
1277template <
typename T>
1289template <
typename T>
1301template <
typename T>
1313template <
typename T>
1331template <
typename T>
1333 return strm << A.
real() <<
' ' << A.
imag();
1347template <
typename Ntype,
typename T>
1350 res.re =
static_cast<Ntype
>(m.re);
1351 res.im =
static_cast<Ntype
>(m.im);
1364template <
typename T>
1366 return A.str(prec, separator);
1378template <
typename T>
1380 std::stringstream ss;
1382 ss <<
"( " << A.
real() <<
", " << A.
imag() <<
" )";
1403template <
typename T>
1411 template <typename A, std::enable_if_t<hila::is_arithmetic<A>::value,
int> = 0>
1429#if defined(_GPU_DEVICE_COMPILE_)
1439template <typename T, std::enable_if_t<hila::contains_complex<T>::value,
int> = 0>
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);
1454template <typename T, std::enable_if_t<hila::contains_complex<T>::value,
int> = 0>
1461template <typename A, typename T, std::enable_if_t<std::is_arithmetic<T>::value,
int> = 0>
1468template <typename A, typename T, std::enable_if_t<std::is_arithmetic<T>::value,
int> = 0>
1475template <
typename A,
typename B>
1485template <typename T, typename A, std::enable_if_t<std::is_arithmetic<A>::value,
int> = 0>
1492template <
typename A,
typename B>
1509template <
typename T>
1511 return exp(z.re) *
Complex<T>(cos(z.im), sin(z.im));
1515template <typename T, std::enable_if_t<hila::is_arithmetic<T>::value,
int> = 0>
1521template <
typename T>
1523 return expi(im.
imag());
1527template <
typename T>
1535template <
typename T>
1539 return pow(r, 0.25) * expi(0.5 * a);
1543template <
typename T>
1547 return pow(r, (1.0 / 6.0)) * expi((1.0 / 3.0) * a);
1551template <
typename A,
typename B>
1553 return exp(p * log(z));
1557template <typename T, typename S, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
1559 return exp(p * log(z));
1563template <typename T, typename S, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
1565 return exp(p * log(z));
1571template <
typename T>
1573 return Complex<T>(sin(z.re) * cosh(z.im), cos(z.re) * sinh(z.im));
1578template <
typename T>
1580 return Complex<T>(cos(z.re) * cosh(z.im), -sin(z.re) * sinh(z.im));
1584template <
typename T>
1586 return sin(z) / cos(z);
1591template <
typename T>
1593 return Complex<T>(sinh(z.re) * cos(z.im), cosh(z.re) * sin(z.im));
1598template <
typename T>
1600 return Complex<T>(cosh(z.re) * cos(z.im), sinh(z.re) * sin(z.im));
1604template <
typename T>
1606 return sinh(z) / cosh(z);
1610template <
typename T>
1612 return -0.5 *
I * (log((
I - z) / (
I + z)));
1616template <
typename T>
1618 return -
I * (log(
I * z + sqrt(1 - z * z)));
1622template <
typename T>
1624 return -
I * (log(z +
I * (sqrt(1 - z * z))));
1628template <
typename T>
1630 return 0.5 * log((1 + z) / (1 - z));
1634template <
typename T>
1636 return log(z + sqrt(1 + z * z));
1640template <
typename T>
1642 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)
auto operator+(const A &a, const Complex< T > &c)
Addition operator Scalar + Complex.
Complex< T > & operator*=(const Complex< A > &lhs) &
Complex multiply 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.
Complex< T > & operator*=(const A a) &
Real multiply assign operator.
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 > 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 > & operator+=(const A &a) &
Real addition assignment operator.
Complex< T > tan(Complex< T > z)
- rely on optimizer to simplify
Complex< T > & operator/=(const Complex< A > &lhs) &
Complex divide assignment operator.
Complex< T > operator--(int)
Multiple decrement operator.
T squarenorm() const
Compute square norm of Complex number.
Complex< T > log(Complex< T > z)
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.
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 > &lhs) &
Complex subtraction 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 > & operator=(const Complex< T > &s) &=default
Assignment operator.
Complex< T > sqrt(Complex< T > z)
Complex< T > asinh(Complex< T > z)
Complex< T > exp(const Complex< T > z)
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.
This files containts definitions for the extended precision class that allows for high precision redu...
Implement hila::swap for gauge fields.
double random()
Real valued uniform random number generator.
Array< n, m, Ntype > cast_to(const Array< n, m, T > &mat)
Array casting operation.
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