ComplexDouble.h

Go to the documentation of this file.

 
#pragma once
 
 
#include <ogdf/basic/internal/intrinsics.h> // IWYU pragma: keep
 
namespace ogdf {
 
namespace sse {
 
#ifdef OGDF_FME_KERNEL_USE_SSE
class ComplexDouble {
public:
    __m128d reg;
 
    inline ComplexDouble() { reg = _mm_setzero_pd(); }
 
    inline ComplexDouble(const ComplexDouble& other) { reg = other.reg; }
 
    inline ComplexDouble(double x) { reg = _mm_setr_pd((x), (0)); }
 
    inline ComplexDouble(double x, double y) { reg = _mm_setr_pd((x), (y)); }
 
    inline ComplexDouble(const double* ptr) { reg = _mm_load_pd(ptr); }
 
    inline ComplexDouble(__m128d r) : reg(r) { }
 
    inline ComplexDouble(float x, float y) { reg = _mm_cvtps_pd(_mm_setr_ps((x), (y), 0, 0)); }
 
 
    inline ComplexDouble operator+(const ComplexDouble& other) const {
        return ComplexDouble(_mm_add_pd(reg, other.reg));
    }
 
    inline ComplexDouble operator-(const ComplexDouble& other) const {
        return ComplexDouble(_mm_sub_pd(reg, other.reg));
    }
 
    inline ComplexDouble operator-(void) const {
        return ComplexDouble(_mm_sub_pd(_mm_setzero_pd(), reg));
    }
 
    inline ComplexDouble operator*(const ComplexDouble& other) const {
        // ---------------------------------
        // | a0*b0 - a1*b1 | a0*b1 + a1*b0 |
        // ---------------------------------
        // bt = | b1 | b0 |
        __m128d b_t = _mm_shuffle_pd(other.reg, other.reg, _MM_SHUFFLE2(0, 1));
        // left = | a0*b0 | a1*b1 |
        __m128d left = _mm_mul_pd(reg, other.reg);
        // right = | a0*b1 | a1*b0 |
        __m128d right = _mm_mul_pd(reg, b_t);
        // left = | a0*b0 | -a1*b1 |
        left = _mm_mul_pd(left, _mm_setr_pd(1.0, -1.0));
        // left = | a0*b0 + (-a1*b1) | a0*b1 + a1*b0 |
        return ComplexDouble(_mm_hadd_pd(left, right));
    }
 
    inline ComplexDouble operator/(const ComplexDouble& other) const {
        // 1/(length(other)^2 * this * other.conj;
        // bt = | b0 | -b1 |
        __m128d conj_reg = _mm_mul_pd(other.reg, _mm_setr_pd(1.0, -1.0));
        // bt = | b1 | b0 |
        __m128d b_t = _mm_shuffle_pd(conj_reg, conj_reg, _MM_SHUFFLE2(0, 1));
        // left = | a0*b0 | a1*b1 |
        __m128d left = _mm_mul_pd(reg, conj_reg);
        // right = | a0*b1 | a1*b0 |
        __m128d right = _mm_mul_pd(reg, b_t);
        // left = | a0*b0 | -a1*b1 |
        left = _mm_mul_pd(left, _mm_setr_pd(1.0, -1.0));
        // left = | a0*b0 + (-a1*b1) | a0*b1 + a1*b0 |
        __m128d product = _mm_hadd_pd(left, right);
        // product = reg*other.reg.conj
        // ell = b0*b0 | b1*b1
        __m128d ell = _mm_mul_pd(conj_reg, conj_reg);
        // ell = b0*b0 + b1*b1 | b0*b0 + b1*b1
        ell = _mm_hadd_pd(ell, ell);
        // ell = length^2 | length^2
        return ComplexDouble(_mm_div_pd(product, ell));
    }
 
    inline ComplexDouble operator*(double scalar) const {
        return ComplexDouble(_mm_mul_pd(reg, _mm_setr_pd(scalar, scalar)));
    }
 
    inline ComplexDouble operator/(double scalar) const {
        //double rcp = 1.0/scalar;
        return ComplexDouble(_mm_div_pd(reg, _mm_setr_pd(scalar, scalar)));
    }
 
    inline ComplexDouble operator*(unsigned int scalar) const {
        return ComplexDouble(_mm_mul_pd(reg, _mm_setr_pd((double)scalar, (double)scalar)));
    }
 
    inline void operator+=(const ComplexDouble& other) { reg = _mm_add_pd(reg, other.reg); }
 
    inline void operator-=(const ComplexDouble& other) { reg = _mm_sub_pd(reg, other.reg); }
 
    inline void operator*=(const ComplexDouble& other) {
        // bt = | b1 | b0 |
        __m128d b_t = _mm_shuffle_pd(other.reg, other.reg, _MM_SHUFFLE2(0, 1));
        // left = | a0*b0 | a1*b1 |
        __m128d left = _mm_mul_pd(reg, other.reg);
        // right = | a0*b1 | a1*b0 |
        __m128d right = _mm_mul_pd(reg, b_t);
        // left = | a0*b0 | -a1*b1 |
        left = _mm_mul_pd(left, _mm_setr_pd(1.0, -1.0));
        // left = | a0*b0 + (-a1*b1) | a0*b1 + a1*b0 |
        reg = _mm_hadd_pd(left, right);
    }
 
    inline void operator*=(double scalar) {
        // (real*scalar, imag*scalar)
        reg = _mm_mul_pd(reg, _mm_setr_pd(scalar, scalar));
    }
 
    inline void operator/=(const ComplexDouble& other) {
        // 1/(length(other)^2 * this * other.conj;
        // bt = | b0 | -b1 |
        __m128d conj_reg = _mm_mul_pd(other.reg, _mm_setr_pd(1.0, -1.0));
        // bt = | b1 | b0 |
        __m128d b_t = _mm_shuffle_pd(conj_reg, conj_reg, _MM_SHUFFLE2(0, 1));
        // left = | a0*b0 | a1*b1 |
        __m128d left = _mm_mul_pd(reg, conj_reg);
        // right = | a0*b1 | a1*b0 |
        __m128d right = _mm_mul_pd(reg, b_t);
        // left = | a0*b0 | -a1*b1 |
        left = _mm_mul_pd(left, _mm_setr_pd(1.0, -1.0));
        // left = | a0*b0 + (-a1*b1) | a0*b1 + a1*b0 |
        __m128d product = _mm_hadd_pd(left, right);
        // ell = b0*b0 | b1*b1
        __m128d ell = _mm_mul_pd(conj_reg, conj_reg);
        // ell = b0*b0 + b1*b1 | b0*b0 + b1*b1
        ell = _mm_hadd_pd(ell, ell);
        // ell = length^2 | length^2
        reg = _mm_div_pd(product, ell);
    }
 
 
    inline double length() const {
        // sqrt(real*real + imag*imag)
        double res;
        __m128d r = _mm_mul_pd(reg, reg);
        r = _mm_hadd_pd(r, _mm_setzero_pd());
        r = _mm_sqrt_sd(r, r);
        _mm_store_sd(&res, r);
        return res;
    }
 
    inline ComplexDouble conj() const {
        // (real, -imag)
        return ComplexDouble(_mm_mul_pd(reg, _mm_setr_pd(1.0, -1.0)));
    }
 
 
    inline void operator=(const ComplexDouble& other) { reg = other.reg; }
 
    inline void operator=(double* ptr) { reg = _mm_load_pd(ptr); }
 
 
    inline void load(const double* ptr) { reg = _mm_load_pd(ptr); }
 
    inline void load_unaligned(const double* ptr) { reg = _mm_loadu_pd(ptr); }
 
    inline void store(double* ptr) const { _mm_store_pd(ptr, reg); }
 
    inline void store_unaligned(double* ptr) const { _mm_storeu_pd(ptr, reg); }
 
};
 
#else
class ComplexDouble {
public:
    double reg[2];
 
    inline ComplexDouble() {
        reg[0] = 0.0;
        reg[1] = 0.0;
    }
 
    inline ComplexDouble(const ComplexDouble& other) {
        reg[0] = other.reg[0];
        reg[1] = other.reg[1];
    }
 
    inline ComplexDouble(double x) {
        reg[0] = x;
        reg[1] = 0;
    }
 
    inline ComplexDouble(double x, double y) {
        reg[0] = x;
        reg[1] = y;
    }
 
    inline ComplexDouble(double* ptr) {
        reg[0] = ptr[0];
        reg[1] = ptr[1];
    }
 
    inline ComplexDouble operator+(const ComplexDouble& other) const {
        return ComplexDouble(reg[0] + other.reg[0], reg[1] + other.reg[1]);
    }
 
    inline ComplexDouble operator-(const ComplexDouble& other) const {
        return ComplexDouble(reg[0] - other.reg[0], reg[1] - other.reg[1]);
    }
 
    inline ComplexDouble operator-(void) const { return ComplexDouble(-reg[0], -reg[1]); }
 
    inline ComplexDouble operator*(const ComplexDouble& other) const {
        return ComplexDouble(reg[0] * other.reg[0] - reg[1] * other.reg[1],
                reg[0] * other.reg[1] + reg[1] * other.reg[0]);
    }
 
    inline ComplexDouble operator/(const ComplexDouble& other) const {
        return (*this) * other.conj() / (other.reg[0] * other.reg[0] + other.reg[1] * other.reg[1]);
    }
 
    inline ComplexDouble operator*(double scalar) const {
        return ComplexDouble(reg[0] * scalar, reg[1] * scalar);
    }
 
    inline ComplexDouble operator/(double scalar) const {
        return ComplexDouble(reg[0] / scalar, reg[1] / scalar);
    }
 
    inline ComplexDouble operator*(unsigned int scalar) const {
        return ComplexDouble(reg[0] * (double)scalar, reg[1] * (double)scalar);
    }
 
    inline void operator+=(const ComplexDouble& other) {
        reg[0] += other.reg[0];
        reg[1] += other.reg[1];
    }
 
    inline void operator-=(const ComplexDouble& other) {
        reg[0] -= other.reg[0];
        reg[1] -= other.reg[1];
    }
 
    inline void operator*=(const ComplexDouble& other) {
        double t[2];
        t[0] = reg[0] * other.reg[0] - reg[1] * other.reg[1];
        t[1] = reg[0] * other.reg[1] + reg[1] * other.reg[0];
        reg[0] = t[0];
        reg[1] = t[1];
    }
 
    inline void operator*=(double scalar) {
        reg[0] *= scalar;
        reg[1] *= scalar;
    }
 
    inline void operator/=(const ComplexDouble& other) {
        ComplexDouble t = other.conj() / (other.reg[0] * other.reg[0] + other.reg[1] * other.reg[1]);
        double r[2];
        r[0] = reg[0] * t.reg[0] - reg[1] * t.reg[1];
        r[1] = reg[0] * t.reg[1] + reg[1] * t.reg[0];
        reg[0] = r[0];
        reg[1] = r[1];
    }
 
 
    inline double length() const {
        // sqrt(real*real + imag*imag)
        return sqrt(reg[0] * reg[0] + reg[1] * reg[1]);
    }
 
    inline ComplexDouble conj() const {
        // (real, -imag)
        return ComplexDouble(reg[0], -reg[1]);
    }
 
 
    inline ComplexDouble& operator=(const ComplexDouble& other) {
        reg[0] = other.reg[0];
        reg[1] = other.reg[1];
        return *this;
    }
 
    inline ComplexDouble& operator=(double* ptr) {
        reg[0] = ptr[0];
        reg[1] = ptr[1];
        return *this;
    }
 
 
    inline void load(const double* ptr) {
        reg[0] = ptr[0];
        reg[1] = ptr[1];
    }
 
    inline void load_unaligned(const double* ptr) {
        reg[0] = ptr[0];
        reg[1] = ptr[1];
    }
 
    inline void store(double* ptr) const {
        ptr[0] = reg[0];
        ptr[1] = reg[1];
    }
 
    inline void store_unaligned(double* ptr) const {
        ptr[0] = reg[0];
        ptr[1] = reg[1];
    }
 
};
 
#endif
}
 
}