backend.hpp

Go to the documentation of this file.

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
 
#ifndef backend_h
#define backend_h
 
#include <algorithm>
#include <vector>
#include <cmath>
 
#include "special_functions.hpp"
#include "register.hpp"
 
namespace backend {
//******************************************************************************
//  Data buffer.
//******************************************************************************
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    class buffer : public std::vector<T> {
    public:
        using std::vector<T>::size;
        using std::vector<T>::data;
        using std::vector<T>::assign;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer() :
        std::vector<T> () {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const size_t s) :
        std::vector<T> (s) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const size_t s, const T d) :
        std::vector<T> (s, d) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const std::vector<T> &d) :
        std::vector<T> (d) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const buffer &d) :
        std::vector<T> (d) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const T min, const T dx, const size_t num) : std::vector<T> (num) {
            for (size_t i = 0; i < num; i++) {
                (*this)[i] = dx*i + min;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void set(const T d) {
            assign(size(), d);
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void set(const std::vector<T> &d) {
            assign(d.cbegin(), d.cend());
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_same() const {
            const T same = (*this)[0];
            for (size_t i = 1, ie = size(); i < ie; i++) {
                if ((*this)[i] != same) {
                    return false;
                }
            }
 
            return true;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_zero() const {
            for (const T &d : *this) {
                if (d != static_cast<T> (0.0)) {
                    return false;
                }
            }
 
            return true;
        }
        
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool has_zero() const {
            for (const T &d : *this) {
                if (d == static_cast<T> (0.0)) {
                    return true;
                }
            }
 
            return false;
        }
        
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_negative() const {
            for (const T &d : *this) {
                if (std::real(d) > std::real(static_cast<T> (0.0))) {
                    return false;
                }
            }
 
            return true;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_even() const {
            for (const T &d : *this) {
                if (std::fmod(std::real(d), std::real(static_cast<T> (2.0)))) {
                    return false;
                }
            }
 
            return true;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_none() const {
            for (const T &d : *this) {
                if (d != static_cast<T> (-1.0)) {
                    return false;
                }
            }
 
            return true;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define apply_op(op) \
for (T &d : *this) { \
    d = op(d);       \
}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void sqrt() {
            apply_op(std::sqrt)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void exp() {
            apply_op(std::exp)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void log() {
            apply_op(std::log)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void sin() {
            apply_op(std::sin)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void cos() {
            apply_op(std::cos)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void real() {
            apply_op(std::real)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void erf() requires(std::floating_point<T>) {
            apply_op(std::erf)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void erfi() requires(jit::complex_scalar<T>) {
            apply_op(special::erfi)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_normal() const {
            for (const T &x : *this) {
                if constexpr (jit::complex_scalar<T>) {
                    if (std::isnan(std::real(x)) || std::isinf(std::real(x)) ||
                        std::isnan(std::imag(x)) || std::isinf(std::imag(x))) {
                        return false;
                    }
                } else {
                    if (std::isnan(x) || std::isinf(x)) {
                        return false;
                    }
                }
            }
            return true;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer<T> index_row(const size_t index, const size_t num_columns) {
            buffer<T> b(num_columns);
            const size_t num_rows = size()/num_columns;
            for (size_t j = 0; j < num_columns; j++) {
                b[j] = (*this)[index*num_rows + j];
            }
            return b;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer<T> index_column(const size_t index, const size_t num_columns) {
            const size_t num_rows = size()/num_columns;
            buffer<T> b(num_rows);
            for (size_t i = 0; i < num_rows; i++) {
                b[i] = (*this)[i*num_rows + index];
            }
            return b;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define row_op(opp, oppeq)                                                  \
if (size() > x.size()) {                                                    \
    assert(size()%x.size() == 0 &&                                          \
           "Vector operand size is not a multiple of matrix operand size"); \
                                                                            \
    const size_t num_columns = size()/x.size();                             \
    const size_t num_rows = x.size();                                       \
    for (size_t i = 0; i < num_rows; i++) {                                 \
        for (size_t j = 0; j < num_columns; j++) {                          \
            (*this)[i*num_columns + j] oppeq x[i];                          \
        }                                                                   \
    }                                                                       \
} else {                                                                    \
    assert(x.size()%size() == 0 &&                                          \
           "Vector operand size is not a multiple of matrix operand size"); \
                                                                            \
    std::vector<T> m(x.size());                                             \
    const size_t num_columns = x.size()/size();                             \
    const size_t num_rows = size();                                         \
    for (size_t i = 0; i < num_rows; i++) {                                 \
        for (size_t j = 0; j < num_columns; j++) {                          \
            m[i*num_columns + j] = (*this)[i] opp x[i*num_columns + j];     \
        }                                                                   \
    }                                                                       \
    *this = m;                                                              \
}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void add_row(const buffer<T> &x) {
            row_op(+, +=)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define col_op(opp, oppeq)                                                  \
if (size() > x.size()) {                                                    \
    assert(size()%x.size() == 0 &&                                          \
           "Vector operand size is not a multiple of matrix operand size"); \
                                                                            \
    const size_t num_columns = size()/x.size();                             \
    const size_t num_rows = x.size();                                       \
    for (size_t i = 0; i < num_rows; i++) {                                 \
        for (size_t j = 0; j < num_columns; j++) {                          \
            (*this)[i*num_columns + j] oppeq x[j];                          \
        }                                                                   \
    }                                                                       \
} else {                                                                    \
    assert(x.size()%size() == 0 &&                                          \
           "Vector operand size is not a multiple of matrix operand size"); \
                                                                            \
    std::vector<T> m(x.size());                                             \
    const size_t num_columns = x.size()/size();                             \
    const size_t num_rows = size();                                         \
    for (size_t i = 0; i < num_rows; i++) {                                 \
        for (size_t j = 0; j < num_columns; j++) {                          \
            m[i*num_columns + j] = (*this)[j] opp x[i*num_columns + j];     \
        }                                                                   \
    }                                                                       \
    *this = m;                                                              \
}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void add_col(const buffer<T> &x) {
            col_op(+, +=)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void subtract_row(const buffer<T> &x) {
            row_op(-, -=)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void subtract_col(const buffer<T> &x) {
            col_op(-, -=)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void multiply_row(const buffer<T> &x) {
            row_op(*, *=)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void multiply_col(const buffer<T> &x) {
            col_op(*, *=)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void divide_row(const buffer<T> &x) {
            row_op(/, /=)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void divide_col(const buffer<T> &x) {
            col_op(/, /=)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        #define row_fn(fn)                                                             \
        if (size() > x.size()) {                                                       \
            assert(size()%x.size() == 0 &&                                             \
                   "Vector operand size is not a multiple of matrix operand size");    \
                                                                                       \
            const size_t num_columns = size()/x.size();                                \
            const size_t num_rows = x.size();                                          \
            for (size_t i = 0; i < num_rows; i++) {                                    \
                for (size_t j = 0; j < num_columns; j++) {                             \
                    (*this)[i*num_columns + j] = fn((*this)[i*num_columns + j], x[i]); \
                }                                                                      \
            }                                                                          \
        } else {                                                                       \
            assert(x.size()%size() == 0 &&                                             \
                   "Vector operand size is not a multiple of matrix operand size");    \
                                                                                       \
            std::vector<T> m(x.size());                                                \
            const size_t num_columns = x.size()/size();                                \
            const size_t num_rows = size();                                            \
            for (size_t i = 0; i < num_rows; i++) {                                    \
                for (size_t j = 0; j < num_columns; j++) {                             \
                    m[i*num_columns + j] = fn((*this)[i], x[i*num_columns + j]);       \
                }                                                                      \
            }                                                                          \
            *this = m;                                                                 \
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void min_row(const buffer<T> &x) {
            row_fn(std::min)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        #define col_fn(fn)                                                             \
        if (size() > x.size()) {                                                       \
            assert(size()%x.size() == 0 &&                                             \
                   "Vector operand size is not a multiple of matrix operand size");    \
                                                                                       \
            const size_t num_columns = size()/x.size();                                \
            const size_t num_rows = x.size();                                          \
            for (size_t i = 0; i < num_rows; i++) {                                    \
                for (size_t j = 0; j < num_columns; j++) {                             \
                    (*this)[i*num_columns + j] = fn((*this)[i*num_columns + j], x[j]); \
                }                                                                      \
            }                                                                          \
        } else {                                                                       \
            assert(x.size()%size() == 0 &&                                             \
                   "Vector operand size is not a multiple of matrix operand size");    \
                                                                                       \
            std::vector<T> m(x.size());                                                \
            const size_t num_columns = x.size()/size();                                \
            const size_t num_rows = size();                                            \
            for (size_t i = 0; i < num_rows; i++) {                                    \
                for (size_t j = 0; j < num_columns; j++) {                             \
                    m[i*num_columns + j] = fn((*this)[j], x[i*num_columns + j]);       \
                }                                                                      \
            }                                                                          \
            *this = m;                                                                 \
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void min_col(const buffer<T> &x) {
            col_fn(std::min)
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void atan_row(const buffer<T> &x) {
            if (size() > x.size()) {
                assert(size()%x.size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                const size_t num_columns = size()/x.size();
                const size_t num_rows = x.size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        if constexpr (jit::complex_scalar<T>) {
                            (*this)[i*num_columns + j] = std::atan(x[i]/(*this)[i*num_columns + j]);
                        } else {
                            (*this)[i*num_columns + j] = std::atan2(x[i], (*this)[i*num_columns + j]);
                        }
                    }
                }
            } else {
                assert(x.size()%size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                std::vector<T> m(x.size());
                const size_t num_columns = x.size()/size();
                const size_t num_rows = size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        if constexpr (jit::complex_scalar<T>) {
                            m[i*num_columns + j] = std::atan(x[i*num_columns + j]/(*this)[i]);
                        } else {
                            m[i*num_columns + j] = std::atan2(x[i*num_columns + j], (*this)[i]);
                        }
                    }
                }
                *this = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void atan_col(const buffer<T> &x) {
            if (size() > x.size()) {
                assert(size()%x.size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                const size_t num_columns = size()/x.size();
                const size_t num_rows = x.size();
                for (size_t i = 0; i < num_columns; i++) {
                    for (size_t j = 0; j < num_rows; j++) {
                        if constexpr (jit::complex_scalar<T>) {
                            (*this)[i*num_columns + j] = std::atan(x[j]/(*this)[i*num_columns + j]);
                        } else {
                            (*this)[i*num_columns + j] = std::atan2(x[j], (*this)[i*num_columns + j]);
                        }
                    }
                }
            } else {
                assert(x.size()%size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                std::vector<T> m(x.size());
                const size_t num_columns = x.size()/size();
                const size_t num_rows = size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        if constexpr (jit::complex_scalar<T>) {
                            m[i*num_columns + j] = std::atan(x[i*num_columns + j]/(*this)[j]);
                        } else {
                            m[i*num_columns + j] = std::atan2(x[i*num_columns + j], (*this)[j]);
                        }
                    }
                }
                *this = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void pow_row(const buffer<T> &x) {
            if (size() > x.size()) {
                assert(size()%x.size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                const size_t num_columns = size()/x.size();
                const size_t num_rows = x.size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        (*this)[i*num_columns + j] = std::pow((*this)[i*num_columns + j], x[i]);
                    }
                }
            } else {
                assert(x.size()%size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                std::vector<T> m(x.size());
                const size_t num_columns = x.size()/size();
                const size_t num_rows = size();
                for (size_t i = 0; i < num_columns; i++) {
                    for (size_t j = 0; j < num_rows; j++) {
                        m[i*num_columns + j] = std::pow((*this)[i], x[i*num_columns + j]);
                    }
                }
                *this = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void pow_col(const buffer<T> &x) {
            if (size() > x.size()) {
                assert(size()%x.size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                const size_t num_columns = size()/x.size();
                const size_t num_rows = x.size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        (*this)[i*num_columns + j] = std::pow((*this)[i*num_columns + j], x[j]);
                    }
                }
            } else {
                assert(x.size()%size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                std::vector<T> m(x.size());
                const size_t num_columns = x.size()/size();
                const size_t num_rows = size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        m[i*num_columns + j] = std::pow((*this)[j], x[i*num_columns + j]);
                    }
                }
                *this = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void hypot_row(const buffer<T> &x) requires(std::floating_point<T>) {
            if (size() > x.size()) {
                assert(size()%x.size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                const size_t num_columns = size()/x.size();
                const size_t num_rows = x.size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        (*this)[i*num_columns + j] = std::hypot((*this)[i*num_columns + j], x[i]);
                    }
                }
            } else {
                assert(x.size()%size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                std::vector<T> m(x.size());
                const size_t num_columns = x.size()/size();
                const size_t num_rows = size();
                for (size_t i = 0; i < num_columns; i++) {
                    for (size_t j = 0; j < num_rows; j++) {
                        m[i*num_columns + j] = std::hypot((*this)[i], x[i*num_columns + j]);
                    }
                }
                *this = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void hypot_col(const buffer<T> &x) requires(std::floating_point<T>) {
            if (size() > x.size()) {
                assert(size()%x.size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                const size_t num_columns = size()/x.size();
                const size_t num_rows = x.size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        (*this)[i*num_columns + j] = std::hypot((*this)[i*num_columns + j], x[j]);
                    }
                }
            } else {
                assert(x.size()%size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                std::vector<T> m(x.size());
                const size_t num_columns = x.size()/size();
                const size_t num_rows = size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        m[i*num_columns + j] = std::hypot((*this)[j], x[i*num_columns + j]);
                    }
                }
                *this = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void copysign_row(const buffer<T> &x) requires(std::floating_point<T>) {
            if (size() > x.size()) {
                assert(size()%x.size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                const size_t num_columns = size()/x.size();
                const size_t num_rows = x.size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        (*this)[i*num_columns + j] = std::copysign((*this)[i*num_columns + j], x[i]);
                    }
                }
            } else {
                assert(x.size()%size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                std::vector<T> m(x.size());
                const size_t num_columns = x.size()/size();
                const size_t num_rows = size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        m[i*num_columns + j] = std::copysign((*this)[i], x[i*num_columns + j]);
                    }
                }
                *this = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void copysign_col(const buffer<T> &x) {
            if (size() > x.size()) {
                assert(size()%x.size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                const size_t num_columns = size()/x.size();
                const size_t num_rows = x.size();
                for (size_t i = 0; i < num_columns; i++) {
                    for (size_t j = 0; j < num_rows; j++) {
                        (*this)[i*num_columns + j] = std::copysign((*this)[i*num_columns + j], x[j]);
                    }
                }
            } else {
                assert(x.size()%size() == 0 &&
                       "Vector operand size is not a multiple of matrix operand size");
 
                std::vector<T> m(x.size());
                const size_t num_columns = x.size()/size();
                const size_t num_rows = size();
                for (size_t i = 0; i < num_rows; i++) {
                    for (size_t j = 0; j < num_columns; j++) {
                        m[i*num_columns + j] = std::copysign((*this)[j], x[i*num_columns + j]);
                    }
                }
                *this = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer<T> operator!() {
            for (T &d : *this) {
                if constexpr (jit::complex_scalar<T>) {
                    assert(d.imag() == 0.0 && "Imaginary part not zero.");
                    d = static_cast<T> (!d.real());
                } else {
                    d = !d;
                }
            }
            return *this;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer<T> if_(const buffer<T> &t, const buffer<T> &f) {
            if (size() == 1) {
                if constexpr (std::floating_point<T>) {
                    return (*this)[0] ? t : f;
                } else {
                    return (*this)[0] != static_cast<T> (0) ? t : f;
                }
            } else {
                if (t.size() == 1) {
                    if (f.size() == 1) {
                        for (T &d : *this) {
                            if constexpr (std::floating_point<T>) {
                                d = d ? t[0] : f[0];
                            } else {
                                d = d != static_cast<T> (0) ? t[0] : f[0];
                            }
                        }
                        return *this;
                    } else {
                        assert(size() == f.size() && "Incompatable buffersize.");
                        for (size_t i = 0, ie = size(); i < ie; i++) {
                            if constexpr (std::floating_point<T>) {
                                (*this)[i] = (*this)[i] ? t[0] : f[i];
                            } else {
                                (*this)[i] = (*this)[i] != static_cast<T> (0) ? t[0] : f[i];
                            }
                        }
                        return *this;
                    }
                } else {
                    assert(size() == t.size() && "Incompatable buffersize.");
                    if (f.size() == 1) {
                        for (size_t i = 0, ie = size(); i < ie; i++) {
                            if constexpr (std::floating_point<T>) {
                                (*this)[i] = (*this)[i] ? t[i] : f[0];
                            } else {
                                (*this)[i] = (*this)[i] != static_cast<T> (0) ? t[i] : f[0];
                            }
                        }
                        return *this;
                    } else {
                        assert(size() == f.size() && "Incompatable buffersize.");
                        for (size_t i = 0, ie = size(); i < ie; i++) {
                            if constexpr (std::floating_point<T>) {
                                (*this)[i] = (*this)[i] ? t[i] : f[i];
                            } else {
                                (*this)[i] = (*this)[i] != static_cast<T> (0) ? t[i] : f[i];
                            }
                        }
                        return *this;
                    }
                }
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define logic_is(op) \
for (T &d : *this) { \
    d = std::op(d);  \
}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    void isinf() {
        logic_is(isinf)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    void isnan() {
        logic_is(isnan)
    }
 
        typedef T base;
    };
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline bool operator==(const buffer<T> &a, const buffer<T> &b) {
        if (a.size() != b.size()) {
            return false;
        }
 
        for (size_t i = 0, ie = a.size(); i < ie; i++) {
            if (a[i] != b[i]) {
                return false;
            }
        }
        return true;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define build_assoc_func(func)                    \
if (b.size() == 1) {                              \
    const T right = b[0];                         \
    for (T &l : a) {                              \
        l = func(std::real(l),                    \
                 std::real(right));               \
    }                                             \
    return a;                                     \
} else if (a.size() == 1) {                       \
    const T left = a[0];                          \
    for (T &r : b) {                              \
        r = func(std::real(r),                    \
                 std::real(left));                \
    }                                             \
    return b;                                     \
}                                                 \
                                                  \
assert(a.size() == b.size() &&                    \
       "Left and right sizes are incompatible."); \
for (size_t i = 0, ie = a.size(); i < ie; i++) {  \
    a[i] = func(std::real(a[i]),                  \
                std::real(b[i]));                 \
}                                                 \
return a;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> max(buffer<T> &a, buffer<T> &b) {
        build_assoc_func(std::max);
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> min(buffer<T> &a, buffer<T> &b) {
        build_assoc_func(std::min);
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define build_assoc_op(op)                        \
if (b.size() == 1) {                              \
    const T right = b[0];                         \
    for (T &l : a) {                              \
        l op right;                               \
    }                                             \
    return a;                                     \
} else if (a.size() == 1) {                       \
    const T left = a[0];                          \
    for (T &r : b) {                              \
        r op left;                                \
    }                                             \
    return b;                                     \
}                                                 \
                                                  \
assert(a.size() == b.size() &&                    \
       "Left and right sizes are incompatible."); \
for (size_t i = 0, ie = a.size(); i < ie; i++) {  \
    a[i] op b[i];                                 \
}                                                 \
return a;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator+(buffer<T> &a, buffer<T> &b) {
        build_assoc_op(+=)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define build_non_assoc_op(op, opeq)              \
if (b.size() == 1) {                              \
    const T right = b[0];                         \
    for (T &l : a) {                              \
        l opeq right;                             \
    }                                             \
    return a;                                     \
} else if (a.size() == 1) {                       \
    const T left = a[0];                          \
    for (T &r : b) {                              \
        r = left op r;                            \
    }                                             \
    return b;                                     \
}                                                 \
                                                  \
assert(a.size() == b.size() &&                    \
       "Left and right sizes are incompatible."); \
for (size_t i = 0, ie = a.size(); i < ie; i++) {  \
    a[i] opeq b[i];                               \
}                                                 \
return a;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator-(buffer<T> &a, buffer<T> &b) {
        build_non_assoc_op(-, -=)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator*(buffer<T> &a, buffer<T> &b) {
        build_assoc_op(*=)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator/(buffer<T> &a, buffer<T> &b) {
        build_non_assoc_op(/, /=)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define build_assoc_func(func)                    \
if (b.size() == 1) {                              \
    const T right = b[0];                         \
    for (T &l : a) {                              \
        l = func(l, right);                       \
    }                                             \
    return a;                                     \
} else if (a.size() == 1) {                       \
    const T left = a[0];                          \
    for (T &r : b) {                              \
        r = func(r, left);                        \
    }                                             \
    return b;                                     \
}                                                 \
                                                  \
assert(a.size() == b.size() &&                    \
       "Left and right sizes are incompatible."); \
for (size_t i = 0, ie = a.size(); i < ie; i++) {  \
    a[i] = func(a[i], b[i]);                      \
}                                                 \
return a;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> hypot(buffer<T> &a, buffer<T> &b) {
        build_assoc_func(std::hypot)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> fma(buffer<T> &a, buffer<T> &b, buffer<T> &c) {
        constexpr bool use_fma = !jit::complex_scalar<T> &&
#ifdef FP_FAST_FMA
                                 true;
#else
                                 false;
#endif
        
        if (a.size() == 1) {
            const T left = a[0];
 
            if (b.size() == 1) {
                const T middle = b[0];
                for (size_t i = 0, ie = c.size(); i < ie; i++) {
                    if constexpr (use_fma) {
                        c[i] = std::fma(left, middle, c[i]);
                    } else {
                        c[i] = left*middle + c[i];
                    }
                }
                return c;
            } else if (c.size() == 1) {
                const T right = c[0];
                for (size_t i = 0, ie = b.size(); i < ie; i++) {
                    if constexpr (use_fma) {
                        b[i] = std::fma(left, b[i], right);
                    } else {
                        b[i] = left*b[i] + right;
                    }
                }
                return b;
            }
 
            assert(b.size() == c.size() &&
                   "Size mismatch between middle and right.");
            for (size_t i = 0, ie = b.size(); i < ie; i++) {
                if constexpr (use_fma) {
                    b[i] = std::fma(left, b[i], c[i]);
                } else {
                    b[i] = left*b[i] + c[i];
                }
            }
            return b;
        } else if (b.size() == 1) {
            const T middle = b[0];
            if (c.size() == 1) {
                const T right = c[0];
                for (size_t i = 0, ie = a.size(); i < ie; i++) {
                    if constexpr (use_fma) {
                        a[i] = std::fma(a[i], middle, right);
                    } else {
                        a[i] = a[i]*middle + right;
                    }
                }
                return a;
            }
 
            assert(a.size() == c.size() &&
                   "Size mismatch between left and right.");
            for (size_t i = 0, ie = a.size(); i < ie; i++) {
                if constexpr (use_fma) {
                    a[i] = std::fma(a[i], middle, c[i]);
                } else {
                    a[i] = a[i]*middle + c[i];
                }
            }
            return a;
        } else if (c.size() == 1) {
            assert(a.size() == b.size() &&
                   "Size mismatch between left and middle.");
            const T right = c[0];
            for (size_t i = 0, ie = a.size(); i < ie; i++) {
                if constexpr (use_fma) {
                    a[i] = std::fma(a[i], b[i], right);
                } else {
                    a[i] = a[i]*b[i] + right;
                }
            }
            return a;
        }
 
        assert(a.size() == b.size() &&
               b.size() == c.size() &&
               a.size() == c.size() &&
               "Left, middle and right sizes are incompatible.");
        for (size_t i = 0, ie = a.size(); i < ie; i++) {
            if constexpr (use_fma) {
                a[i] = std::fma(a[i], b[i], c[i]);
            } else {
                a[i] = a[i]*b[i] + c[i];
            }
        }
        return a;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> operator%(buffer<T> &a, buffer<T> &b) {
        if (b.size() == 1) {
            const T right = b[0];
            for (size_t i = 0, ie = a.size(); i < ie; i++) {
                a[i] = std::fmod(a[i], right);
            }
            return a;
        } else if (a.size() == 1) {
            const T left = a[0];
            for (size_t i = 0, ie = b.size(); i < ie; i++) {
                b[i] = std::fmod(left, b[i]);
            }
            return b;
        }
 
        assert(a.size() == b.size() &&
               "Left and right sizes are incompatible.");
        for (size_t i = 0, ie = a.size(); i < ie; i++) {
            a[i] = std::fmod(a[i], b[i]);
        }
        return a;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define logic_op(op)                                 \
if (b.size() == 1) {                                 \
    const T right = b[0];                            \
    for (size_t i = 0, ie = a.size(); i < ie; i++) { \
        a[i] = static_cast<T> (a[i] op right);       \
    }                                                \
    return a;                                        \
} else if (a.size() == 1) {                          \
    const T left = a[0];                             \
    for (size_t i = 0, ie = b.size(); i < ie; i++) { \
        b[i] = static_cast<T> (left op b[i]);        \
    }                                                \
    return b;                                        \
}                                                    \
                                                     \
assert(a.size() == b.size() &&                       \
       "Left and right sizes are incompatible.");    \
for (size_t i = 0, ie = a.size(); i < ie; i++) {     \
    a[i] = static_cast<T> (a[i] op b[i]);            \
}                                                    \
return a;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator==(buffer<T> &a, buffer<T> &b) {
        logic_op(==)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator!=(buffer<T> &a, buffer<T> &b) {
        logic_op(!=)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> operator>(buffer<T> &a, buffer<T> &b) {
        logic_op(>)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> operator<(buffer<T> &a, buffer<T> &b) {
        logic_op(<)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> operator>=(buffer<T> &a, buffer<T> &b) {
        logic_op(>=)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> operator<=(buffer<T> &a, buffer<T> &b) {
        logic_op(<=)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> operator&&(buffer<T> &a, buffer<T> &b) {
        logic_op(&&)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> operator||(buffer<T> &a, buffer<T> &b) {
        logic_op(||)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
#define branch_fn(fn)                                \
if (b.size() == 1) {                                 \
    const T right = b[0];                            \
    for (size_t i = 0, ie = a.size(); i < ie; i++) { \
        a[i] = fn(a[i], right);                      \
    }                                                \
    return a;                                        \
} else if (a.size() == 1) {                          \
    const T left = a[0];                             \
    for (size_t i = 0, ie = b.size(); i < ie; i++) { \
        b[i] = fn(left, b[i]);                       \
    }                                                \
    return b;                                        \
}                                                    \
                                                     \
assert(a.size() == b.size() &&                       \
       "Left and right sizes are incompatible.");    \
for (size_t i = 0, ie = a.size(); i < ie; i++) {     \
    a[i] = fn(a[i], b[i]);                           \
}                                                    \
return a;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> min(buffer<T> &a, buffer<T> &b) {
        branch_fn(std::min)
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> pow(buffer<T> &base, buffer<T> &exponent) {
        if (exponent.size() == 1) {
            const T right = exponent[0];
            if (std::imag(right) == 0) {
                const int64_t right_int = static_cast<int64_t> (std::real(right));
                if (std::real(right) - right_int) {
                    if (right == static_cast<T> (0.5)) {
                        base.sqrt();
                        return base;
                    }
 
                    for (size_t i = 0, ie = base.size(); i < ie; i++) {
                        base[i] = std::pow(base[i], right);
                    }
                    return base;
                }
 
                if (right_int > 0) {
                    for (size_t i = 0, ie = base.size(); i < ie; i++) {
                        const T left = base[i];
                        for (size_t j = 0, je = right_int - 1; j < je; j++) {
                            base[i] *= left;
                        }
                    }
                    return base;
                } else if (right_int == 0) {
                    for (size_t i = 0, ie = base.size(); i < ie; i++) {
                        base[i] = 1.0;
                    }
                    return base;
                } else {
                    for (size_t i = 0, ie = base.size(); i < ie; i++) {
                        const T left = static_cast<T> (1.0)/base[i];
                        base[i] = left;
                        for (size_t j = 0, je = std::abs(right_int) - 1; j < je; j++) {
                            base[i] *= left;
                        }
                    }
                    return base;
                }
            } else {
                for (size_t i = 0, ie = base.size(); i < ie; i++) {
                    base[i] = std::pow(base[i], right);
                }
                return base;
            }
        } else if (base.size() == 1) {
            const T left = base[0];
            for (size_t i = 0, ie = exponent.size(); i < ie; i++) {
                exponent[i] = std::pow(left, exponent[i]);
            }
            return exponent;
        }
 
        assert(base.size() == exponent.size() &&
               "Left and right sizes are incompatible.");
        for (size_t i = 0, ie = base.size(); i < ie; i++) {
            base[i] = std::pow(base[i], exponent[i]);
        }
        return base;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> atan(buffer<T> &x, buffer<T> &y) {
        if (y.size() == 1) {
            const T right = y[0];
            for (size_t i = 0, ie = x.size(); i < ie; i++) {
                if constexpr (jit::complex_scalar<T>) {
                    x[i] = std::atan(right/x[i]);
                } else {
                    x[i] = std::atan2(right, x[i]);
                }
            }
            return x;
        } else if (x.size() == 1) {
            const T left = x[0];
            for (size_t i = 0, ie = y.size(); i < ie; i++) {
                if constexpr (jit::complex_scalar<T>) {
                    y[i] = std::atan(y[i]/left);
                } else {
                    y[i] = std::atan2(y[i], left);
                }
            }
            return y;
        }
 
        assert(x.size() == y.size() &&
               "Left and right sizes are incompatible.");
        for (size_t i = 0, ie = x.size(); i < ie; i++) {
            if constexpr (jit::complex_scalar<T>) {
                x[i] = std::atan(y[i]/x[i]);
            } else {
                x[i] = std::atan2(y[i], x[i]);
            }
        }
        return x;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<std::floating_point T>
    inline buffer<T> copysign(buffer<T> &x,
                              buffer<T> &y) {
        if (y.size() == 1) {
            const T right = y[0];
            for (size_t i = 0, ie = x.size(); i < ie; i++) {
                x[i] = std::copysign(x[i], right);
            }
            return x;
        } else if (x.size() == 1) {
            const T left = x[0];
            for (size_t i = 0, ie = y.size(); i < ie; i++) {
                y[i] = std::copysign(left, y[i]);
            }
            return y;
        }
 
        assert(x.size() == y.size() &&
               "Left and right sizes are incompatible.");
        for (size_t i = 0, ie = x.size(); i < ie; i++) {
            x[i] = std::copysign(x[i], y[i]);
        }
        return x;
    }
}
 
#endif /* backend_h */