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 {
    private:
        std::vector<T> memory;
 
    public:
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer() :
        memory() {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const size_t s) :
        memory(s) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const size_t s, const T d) :
        memory(s, d) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const std::vector<T> &d) :
        memory(d) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        buffer(const buffer &d) :
        memory(d.memory) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        T &operator[] (const size_t index) {
            return memory[index];
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        const T &operator[] (const size_t index) const {
            return memory[index];
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        const T at(const size_t index) const {
            return memory.at(index);
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void set(const T d) {
            memory.assign(memory.size(), d);
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void set(const std::vector<T> &d) {
            memory.assign(d.cbegin(), d.cend());
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        size_t size() const {
            return memory.size();
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_same() const {
            const T same = memory.at(0);
            for (size_t i = 1, ie = memory.size(); i < ie; i++) {
                if (memory.at(i) != same) {
                    return false;
                }
            }
 
            return true;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_zero() const {
            for (const T &d : memory) {
                if (d != static_cast<T> (0.0)) {
                    return false;
                }
            }
 
            return true;
        }
        
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool has_zero() const {
            for (const T &d : memory) {
                if (d == static_cast<T> (0.0)) {
                    return true;
                }
            }
 
            return false;
        }
        
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_negative() const {
            for (const T &d : memory) {
                if (std::real(d) > std::real(static_cast<T> (0.0))) {
                    return false;
                }
            }
 
            return true;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_even() const {
            for (const T &d : memory) {
                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 : memory) {
                if (d != static_cast<T> (-1.0)) {
                    return false;
                }
            }
 
            return true;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void sqrt() {
            for (T &d : memory) {
                d = std::sqrt(d);
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void exp() {
            for (T &d : memory) {
                d = std::exp(d);
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void log() {
            for (T &d : memory) {
                d = std::log(d);
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void sin() {
            for (T &d : memory) {
                d = std::sin(d);
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void cos() {
            for (T &d : memory) {
                d = std::cos(d);
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void erfi() requires(jit::complex_scalar<T>) {
            for (T &d : memory) {
                d = special::erfi(d);
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        T *data() {
            return memory.data();
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        bool is_normal() const {
            for (const T &x : memory) {
                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] = memory[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] = memory[i*num_rows + index];
            }
            return b;
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void add_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++) {
                        memory[i*num_rows + 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] = memory[i] + x[i*num_columns + j];
                    }
                }
                memory = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void add_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++) {
                        memory[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] = memory[j] + x[i*num_columns + j];
                    }
                }
                memory = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void subtract_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++) {
                        memory[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] = memory[i] - x[i*num_columns + j];
                    }
                }
                memory = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void subtract_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++) {
                        memory[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] = memory[j] - x[i*num_columns + j];
                    }
                }
                memory = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void multiply_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++) {
                        memory[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] = memory[i]*x[i*num_columns + j];
                    }
                }
                memory = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void multiply_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++) {
                        memory[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] = memory[j]*x[i*num_columns + j];
                    }
                }
                memory = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void divide_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++) {
                        memory[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] = memory[i]/x[i*num_columns + j];
                    }
                }
                memory = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        void divide_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++) {
                        memory[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] = memory[j]/x[i*num_columns + j];
                    }
                }
                memory = m;
            }
        }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
        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>) {
                            memory[i*num_columns + j] = std::atan(x[i]/memory[i*num_columns + j]);
                        } else {
                            memory[i*num_columns + j] = std::atan2(x[i], memory[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]/memory[i]);
                        } else {
                            m[i*num_columns + j] = std::atan2(x[i*num_columns + j], memory[i]);
                        }
                    }
                }
                memory = 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>) {
                            memory[i*num_columns + j] = std::atan(x[j]/memory[i*num_columns + j]);
                        } else {
                            memory[i*num_columns + j] = std::atan2(x[j], memory[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]/memory[j]);
                        } else {
                            m[i*num_columns + j] = std::atan2(x[i*num_columns + j], memory[j]);
                        }
                    }
                }
                memory = 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++) {
                        memory[i*num_columns + j] = std::pow(memory[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(memory[i], x[i*num_columns + j]);
                    }
                }
                memory = 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++) {
                        memory[i*num_columns + j] = std::pow(memory[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(memory[j], x[i*num_columns + j]);
                    }
                }
                memory = m;
            }
        }
 
        typedef T base;
    };
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator+(buffer<T> &a,
                               buffer<T> &b) {
        if (b.size() == 1) {
            const T right = b.at(0);
            for (size_t i = 0, ie = a.size(); i < ie; i++) {
                a[i] += right;
            }
            return a;
        } else if (a.size() == 1) {
            const T left = a.at(0);
            for (size_t i = 0, ie = b.size(); i < ie; i++) {
                b[i] += 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] += b.at(i);
        }
        return a;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    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.at(i) != b.at(i)) {
                return false;
            }
        }
        return true;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator-(buffer<T> &a,
                               buffer<T> &b) {
        if (b.size() == 1) {
            const T right = b.at(0);
            for (size_t i = 0, ie = a.size(); i < ie; i++) {
                a[i] -= right;
            }
            return a;
        } else if (a.size() == 1) {
            const T left = a.at(0);
            for (size_t i = 0, ie = b.size(); i < ie; i++) {
                b[i] = left - b.at(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] -= b.at(i);
        }
        return a;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator*(buffer<T> &a,
                               buffer<T> &b) {
        if (b.size() == 1) {
            const T right = b.at(0);
            for (size_t i = 0, ie = a.size(); i < ie; i++) {
                a[i] *= right;
            }
            return a;
        } else if (a.size() == 1) {
            const T left = a.at(0);
            for (size_t i = 0, ie = b.size(); i < ie; i++) {
                b[i] *= 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] *= b.at(i);
        }
        return a;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> operator/(buffer<T> &a,
                               buffer<T> &b) {
        if (b.size() == 1) {
            const T right = b.at(0);
            for (size_t i = 0, ie = a.size(); i < ie; i++) {
                a[i] /= right;
            }
            return a;
        } else if (a.size() == 1) {
            const T left = a.at(0);
            for (size_t i = 0, ie = b.size(); i < ie; i++) {
                b[i] = left/b.at(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] /= b.at(i);
        }
        return a;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    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.at(0);
 
            if (b.size() == 1) {
                const T middle = b.at(0);
                for (size_t i = 0, ie = c.size(); i < ie; i++) {
                    if constexpr (use_fma) {
                        c[i] = std::fma(left, middle, c.at(i));
                    } else {
                        c[i] = left*middle + c.at(i);
                    }
                }
                return c;
            } else if (c.size() == 1) {
                const T right = c.at(0);
                for (size_t i = 0, ie = b.size(); i < ie; i++) {
                    if constexpr (use_fma) {
                        b[i] = std::fma(left, b.at(i), right);
                    } else {
                        b[i] = left*b.at(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.at(i), c.at(i));
                } else {
                    b[i] = left*b.at(i) + c.at(i);
                }
            }
            return b;
        } else if (b.size() == 1) {
            const T middle = b.at(0);
            if (c.size() == 1) {
                const T right = c.at(0);
                for (size_t i = 0, ie = a.size(); i < ie; i++) {
                    if constexpr (use_fma) {
                        a[i] = std::fma(a.at(i), middle, right);
                    } else {
                        a[i] = a.at(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.at(i), middle, c.at(i));
                } else {
                    a[i] = a.at(i)*middle + c.at(i);
                }
            }
            return a;
        } else if (c.size() == 1) {
            assert(a.size() == b.size() &&
                   "Size mismatch between left and middle.");
            const T right = c.at(0);
            for (size_t i = 0, ie = a.size(); i < ie; i++) {
                if constexpr (use_fma) {
                    a[i] = std::fma(a.at(i), b.at(i), right);
                } else {
                    a[i] = a.at(i)*b.at(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.at(i), b.at(i), c.at(i));
            } else {
                a[i] = a.at(i)*b.at(i) + c.at(i);
            }
        }
        return a;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    template<jit::float_scalar T>
    inline buffer<T> pow(buffer<T> &base,
                         buffer<T> &exponent) {
        if (exponent.size() == 1) {
            const T right = exponent.at(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.at(i), right);
                    }
                    return base;
                }
 
                if (right_int > 0) {
                    for (size_t i = 0, ie = base.size(); i < ie; i++) {
                        const T left = base.at(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.at(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.at(i), right);
                }
                return base;
            }
        } else if (base.size() == 1) {
            const T left = base.at(0);
            for (size_t i = 0, ie = exponent.size(); i < ie; i++) {
                exponent[i] = std::pow(left, exponent.at(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.at(i), exponent.at(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.at(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.at(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;
    }
}
 
#endif /* backend_h */