gpu_complex.hpp

/***************************************************************************\
* Copyright (c) 2008-2024, Sofie Martins, Erik Kjellgren, Claudio Pica      *
* All rights reserved.                                                      *
\***************************************************************************/
 
#ifndef GPU_COMPLEX_H
#define GPU_COMPLEX_H
 
#include "Core/gpu.h"
 
template <class T> struct __align__(sizeof(T)) hr_complex_t {
    T re, im;
 
    constexpr visible hr_complex_t(const T a = (T)0, const T b = (T)0)
        : re(a)
        , im(b) {
    }
 
    template <class L>
    constexpr visible hr_complex_t(const hr_complex_t<L> a)
        : re(a.re)
        , im(a.im) {
    }
 
    template <class L> hr_complex_t visible inline __attribute__((always_inline)) &operator=(const L x) {
        re = x;
        im = 0;
        return *this;
    }
 
    template <class L> hr_complex_t visible inline __attribute__((always_inline)) &operator=(const hr_complex_t<L> x) {
        re = x.re;
        im = x.im;
        return *this;
    }
 
    hr_complex_t visible inline __attribute__((always_inline)) operator-(void) const {
        return hr_complex_t(-re, -im);
    }
 
    template <class L> hr_complex_t visible inline __attribute__((always_inline)) &operator+=(const L x) {
        re += x;
        return *this;
    }
 
    template <class L> hr_complex_t visible inline __attribute__((always_inline)) &operator+=(const hr_complex_t<L> x) {
        re += x.re;
        im += x.im;
        return *this;
    }
 
    template <class L> hr_complex_t visible inline __attribute__((always_inline)) &operator-=(const L x) {
        re -= x;
        return *this;
    }
 
    template <class L> auto visible inline __attribute__((always_inline)) &operator-=(const hr_complex_t<L> x) {
        re -= x.re;
        im -= x.im;
        return *this;
    }
 
    template <class L> auto visible inline __attribute__((always_inline)) &operator*=(const L x) {
        re *= x;
        im *= x;
        return *this;
    }
 
    template <class L> auto visible inline __attribute__((always_inline)) &operator*=(const hr_complex_t<L> x) {
        T re_new = re * x.re - im * x.im;
        T im_new = re * x.im + im * x.re;
        re = re_new;
        im = im_new;
        return *this;
    }
 
    template <class L> auto visible inline __attribute__((always_inline)) &operator/=(const L x) {
        re /= x;
        im /= x;
        return *this;
    }
 
    template <class L> auto visible inline __attribute__((always_inline)) &operator/=(const hr_complex_t<L> x) {
        T re_new = (re * x.re + im * x.im) / (x.re * x.re + x.im * x.im);
        T im_new = (im * x.re - re * x.im) / (x.re * x.re + x.im * x.im);
        re = re_new;
        im = im_new;
        return *this;
    }
 
    hr_complex_t visible inline __attribute__((always_inline)) conj() {
        return hr_complex_t(re, -im);
    }
};
 
template <class L, class T> auto visible inline __attribute__((always_inline)) operator+(const hr_complex_t<L> &c, const T x) {
    return hr_complex_t(c.re + x, c.im + (T)0);
}
 
template <class L, class T> auto visible inline __attribute__((always_inline)) operator+(const L x, const hr_complex_t<T> &c) {
    return hr_complex_t(c.re + x, c.im + (L)0);
}
 
template <class L, class T>
auto visible inline __attribute__((always_inline)) operator+(const hr_complex_t<L> &x, const hr_complex_t<T> &c) {
    return hr_complex_t(c.re + x.re, c.im + x.im);
}
 
template <class T, class L> auto visible inline __attribute__((always_inline)) operator-(const hr_complex_t<T> &c, const L x) {
    return hr_complex_t(c.re - x, c.im - (L)0);
}
 
template <class T, class L> auto visible inline __attribute__((always_inline)) operator-(const T x, const hr_complex_t<L> &c) {
    return hr_complex_t(x - c.re, (T)0 - c.im);
}
 
template <class T, class L>
auto visible inline __attribute__((always_inline)) operator-(const hr_complex_t<L> &x, const hr_complex_t<T> &c) {
    return hr_complex_t(x.re - c.re, x.im - c.im);
}
 
template <class T, class L> auto visible inline __attribute__((always_inline)) operator*(const hr_complex_t<T> &c, const L x) {
    return hr_complex_t(c.re * x, c.im * x);
}
 
template <class T, class L> auto visible inline __attribute__((always_inline)) operator*(const T x, const hr_complex_t<L> &c) {
    return hr_complex_t(c.re * x, c.im * x);
}
 
template <class T, class L>
auto visible inline __attribute__((always_inline)) operator*(const hr_complex_t<T> &a, const hr_complex_t<L> &b) {
    return hr_complex_t(a.re * b.re - b.im * a.im, a.im * b.re + b.im * a.re);
}
 
template <class T, class L> auto visible inline __attribute__((always_inline)) operator/(const hr_complex_t<T> &c, const L x) {
    return hr_complex_t(c.re / x, c.im / x);
}
 
template <class T, class L> auto visible inline __attribute__((always_inline)) operator/(const T x, const hr_complex_t<L> &c) {
    return hr_complex_t(x * c.re / (c.re * c.re + c.im * c.im), -x * c.im / (c.re * c.re + c.im * c.im));
}
 
template <class T, class L>
auto visible inline __attribute__((always_inline)) operator/(const hr_complex_t<T> x, const hr_complex_t<L> c) {
    return hr_complex_t((x.re * c.re + x.im * c.im) / (c.re * c.re + c.im * c.im),
                        (x.im * c.re - x.re * c.im) / (c.re * c.re + c.im * c.im));
}
 
typedef struct hr_complex_t<double> hr_complex;
typedef struct hr_complex_t<float> hr_complex_flt;
 
#ifndef HIP
#define I (hr_complex_flt(0.0f, 1.0f))
#else
constexpr hr_complex I = hr_complex_flt(0.0f, 1.0f);
#endif
#define creal(a) ((a).re)
#define cimag(a) ((a).im)
#define conj(a) ((a).conj())
 
visible double carg(hr_complex c);
 
visible hr_complex cpow(hr_complex c, double pow);
 
#endif