vox_geometry_rust 0.1.2

/*
 * // Copyright (c) 2021 Feng Yang
 * //
 * // I am making my contributions/submissions to this project solely in my
 * // personal capacity and am not conveying any rights to any intellectual
 * // property of any third parties.
 */

use crate::usize2::USize2;
use crate::operators::*;
use crate::vector_expression::*;
use std::ops::{Sub, Neg, Add, Mul, Div};


///
/// # Base class for matrix expression.
///
/// Matrix expression is a meta type that enables template expression pattern.
///
/// - tparam T  Real number type.
/// - tparam E  Subclass type.
///
pub trait MatrixExpression {
    /// Size of the matrix.
    fn size(&self) -> USize2;

    /// Number of rows.
    fn rows(&self) -> usize;

    /// Number of columns.
    fn cols(&self) -> usize;

    /// Returns matrix element at (i, j).
    fn eval(&self, x: usize, y: usize) -> f64;
}

//--------------------------------------------------------------------------------------------------
///
/// # Constant matrix expression.
///
/// This matrix expression represents a constant matrix.
///
/// - tparam T  Real number type.
///
pub struct MatrixConstant {
    _m: usize,
    _n: usize,
    _c: f64,
}

impl MatrixConstant {
    /// Constructs m x n constant matrix expression.
    pub fn new(m: usize, n: usize, c: f64) -> MatrixConstant {
        return MatrixConstant {
            _m: m,
            _n: n,
            _c: c,
        };
    }
}

impl MatrixExpression for MatrixConstant {
    /// Size of the matrix.
    fn size(&self) -> USize2 {
        return USize2::new(self.rows(), self.cols());
    }

    /// Number of rows.
    fn rows(&self) -> usize {
        return self._m;
    }

    /// Number of columns.
    fn cols(&self) -> usize {
        return self._n;
    }

    /// Returns matrix element at (i, j).
    fn eval(&self, _i: usize, _j: usize) -> f64 {
        return self._c;
    }
}

//--------------------------------------------------------------------------------------------------
///
/// # Identity matrix expression.
///
/// This matrix expression represents an identity matrix.
///
/// - tparam T  Real number type.
///
pub struct MatrixIdentity {
    _m: usize,
}

impl MatrixIdentity {
    /// Constructs m x m identity matrix expression.
    pub fn new(m: usize) -> MatrixIdentity {
        return MatrixIdentity {
            _m: m
        };
    }
}

impl MatrixExpression for MatrixIdentity {
    /// Size of the matrix.
    fn size(&self) -> USize2 {
        return USize2::new(self._m, self._m);
    }

    /// Number of rows.
    fn rows(&self) -> usize {
        return self._m;
    }

    /// Number of columns.
    fn cols(&self) -> usize {
        return self._m;
    }

    /// Returns matrix element at (i, j).
    fn eval(&self, i: usize, j: usize) -> f64 {
        return match i == j {
            true => 1.0,
            false => 0.0
        };
    }
}

//--------------------------------------------------------------------------------------------------
/// # MatrixUnaryOp

///
/// # Matrix expression for unary operation.
///
/// This matrix expression represents an unary matrix operation that takes
/// single input matrix expression.
///
/// - tparam T   Real number type.
/// - tparam E   Input expression type.
/// - tparam Op  Unary operation.
///
pub struct MatrixUnaryOp<E: MatrixExpression, Op: UnaryOp> {
    _u: E,
    _op: Op,
}

impl<E: MatrixExpression, Op: UnaryOp> MatrixUnaryOp<E, Op> {
    /// Constructs unary operation expression for given input expression.
    pub fn new(u: E) -> MatrixUnaryOp<E, Op> {
        return MatrixUnaryOp {
            _u: u,
            _op: Op::new(),
        };
    }
}

impl<E: MatrixExpression, Op: UnaryOp> MatrixExpression for MatrixUnaryOp<E, Op> {
    /// Size of the matrix.
    fn size(&self) -> USize2 {
        return self._u.size();
    }

    /// Number of rows.
    fn rows(&self) -> usize {
        return self._u.rows();
    }

    /// Number of columns.
    fn cols(&self) -> usize {
        return self._u.cols();
    }

    /// Returns matrix element at (i, j).
    fn eval(&self, i: usize, j: usize) -> f64 {
        return self._op.eval(self._u.eval(i, j));
    }
}

//--------------------------------------------------------------------------------------------------
///
/// # Diagonal matrix expression.
///
/// This matrix expression represents a diagonal matrix for given input matrix
/// expression.
///
/// - tparam T  Real number type.
/// - tparam E  Input expression type.
///
pub struct MatrixDiagonal<E: MatrixExpression> {
    _u: E,
    _is_diag: bool,
}

impl<E: MatrixExpression> MatrixDiagonal<E> {
    /// Constructs diagonal matrix expression for given input expression.
    /// - parameter: is_diag - True for diagonal matrix, false for off-diagonal.
    pub fn new(u: E, is_diag: bool) -> MatrixDiagonal<E> {
        return MatrixDiagonal {
            _u: u,
            _is_diag: is_diag,
        };
    }
}

impl<E: MatrixExpression> MatrixExpression for MatrixDiagonal<E> {
    /// Size of the matrix.
    fn size(&self) -> USize2 {
        return self._u.size();
    }

    /// Number of rows.
    fn rows(&self) -> usize {
        return self._u.rows();
    }

    /// Number of columns.
    fn cols(&self) -> usize {
        return self._u.cols();
    }

    /// Returns matrix element at (i, j).
    fn eval(&self, i: usize, j: usize) -> f64 {
        return if self._is_diag {
            match i == j {
                true => self._u.eval(i, j),
                false => 0.0
            }
        } else {
            match i != j {
                true => self._u.eval(i, j),
                false => 0.0
            }
        };
    }
}

//--------------------------------------------------------------------------------------------------
///
/// # Triangular matrix expression.
///
/// This matrix expression represents a triangular matrix for given input matrix
/// expression.
///
/// - tparam T  Real number type.
/// - tparam E  Input expression type.
///
pub struct MatrixTriangular<E: MatrixExpression> {
    _u: E,
    _is_upper: bool,
    _is_strict: bool,
}

impl<E: MatrixExpression> MatrixTriangular<E> {
    /// Constructs triangular matrix expression for given input expression.
    /// - parameter: is_upper - True for upper tri matrix, false for lower tri matrix.
    /// - parameter: is_strict - True for strictly upper/lower triangular matrix.
    pub fn new(u: E, is_upper: bool, is_strict: bool) -> MatrixTriangular<E> {
        return MatrixTriangular {
            _u: u,
            _is_upper: is_upper,
            _is_strict: is_strict,
        };
    }
}

impl<E: MatrixExpression> MatrixExpression for MatrixTriangular<E> {
    /// Size of the matrix.
    fn size(&self) -> USize2 {
        return self._u.size();
    }

    /// Number of rows.
    fn rows(&self) -> usize {
        return self._u.rows();
    }

    /// Number of columns.
    fn cols(&self) -> usize {
        return self._u.cols();
    }

    /// Returns matrix element at (i, j).
    fn eval(&self, i: usize, j: usize) -> f64 {
        return if i < j {
            match self._is_upper {
                true => self._u.eval(i, j),
                false => 0.0
            }
        } else if i > j {
            match !self._is_upper {
                true => self._u.eval(i, j),
                false => 0.0
            }
        } else {
            match !self._is_strict {
                true => self._u.eval(i, j),
                false => 0.0
            }
        };
    }
}

//--------------------------------------------------------------------------------------------------
/// # MatrixBinaryOp

///
/// # Matrix expression for binary operation.
///
/// This matrix expression represents a binary matrix operation that takes
/// two input matrix expressions.
///
/// - tparam T   Real number type.
/// - tparam E1  First input expression type.
/// - tparam E2  Second input expression type.
/// - tparam Op  Binary operation.
///
pub struct MatrixBinaryOp<E1: MatrixExpression, E2: MatrixExpression, Op: BinaryOp> {
    _u: E1,
    _v: E2,
    _op: Op,
}

impl<E1: MatrixExpression, E2: MatrixExpression, Op: BinaryOp> MatrixBinaryOp<E1, E2, Op> {
    /// Constructs binary operation expression for given input matrix
    /// expressions.
    pub fn new(u: E1, v: E2) -> MatrixBinaryOp<E1, E2, Op> {
        return MatrixBinaryOp {
            _u: u,
            _v: v,
            _op: Op::new(),
        };
    }
}

impl<E1: MatrixExpression, E2: MatrixExpression, Op: BinaryOp> MatrixExpression for MatrixBinaryOp<E1, E2, Op> {
    /// Size of the matrix.
    fn size(&self) -> USize2 {
        return self._v.size();
    }

    /// Number of rows.
    fn rows(&self) -> usize {
        return self._v.rows();
    }

    /// Number of columns.
    fn cols(&self) -> usize {
        return self._v.cols();
    }

    /// Returns matrix element at (i, j).
    fn eval(&self, i: usize, j: usize) -> f64 {
        return self._op.eval(self._u.eval(i, j), self._v.eval(i, j));
    }
}

//--------------------------------------------------------------------------------------------------
///
/// # Matrix expression for matrix-scalar binary operation.
///
/// This matrix expression represents a binary matrix operation that takes
/// one input matrix expression and one scalar.
///
/// - tparam T   Real number type.
/// - tparam E   Input expression type.
/// - tparam Op  Binary operation.
///
pub struct MatrixScalarBinaryOp<E: MatrixExpression, Op: BinaryOp> {
    _u: E,
    _v: f64,
    _op: Op,
}

impl<E: MatrixExpression, Op: BinaryOp> MatrixScalarBinaryOp<E, Op> {
    /// Constructs a binary expression for given matrix and scalar.
    pub fn new(u: E, v: f64) -> MatrixScalarBinaryOp<E, Op> {
        return MatrixScalarBinaryOp {
            _u: u,
            _v: v,
            _op: Op::new(),
        };
    }
}

impl<E: MatrixExpression, Op: BinaryOp> MatrixExpression for MatrixScalarBinaryOp<E, Op> {
    fn size(&self) -> USize2 {
        return self._u.size();
    }

    fn rows(&self) -> usize {
        return self._u.rows();
    }

    fn cols(&self) -> usize {
        return self._u.cols();
    }

    /// Returns matrix element at (i, j).
    fn eval(&self, i: usize, j: usize) -> f64 {
        return self._op.eval(self._u.eval(i, j), self._v);
    }
}

//--------------------------------------------------------------------------------------------------
///
/// # Vector expression for matrix-vector multiplication.
///
/// This vector expression represents a matrix-vector operation that takes
/// one input matrix expression and one vector expression.
///
/// - tparam T   Element value type.
/// - tparam ME  Matrix expression.
/// - tparam VE  Vector expression.
///
pub struct MatrixVectorMul<ME: MatrixExpression, VE: VectorExpression> {
    _m: ME,
    _v: VE,
}

impl<ME: MatrixExpression, VE: VectorExpression> MatrixVectorMul<ME, VE> {
    pub fn new(m: ME, v: VE) -> MatrixVectorMul<ME, VE> {
        return MatrixVectorMul {
            _m: m,
            _v: v,
        };
    }

    /// Returns vector element at i.
    pub fn eval(&self, i: usize) -> f64 {
        let mut sum: f64 = 0.0;
        let n = self._m.cols();
        for j in 0..n {
            sum += self._m.eval(i, j) * self._v.eval(j);
        }
        return sum;
    }

    /// Size of the vector.
    pub fn size(&self) -> usize {
        return self._v.size();
    }
}

impl<ME: MatrixExpression, VE: VectorExpression> MatrixExpression for MatrixVectorMul<ME, VE> {
    fn size(&self) -> USize2 {
        unimplemented!();
    }

    fn rows(&self) -> usize {
        unimplemented!();
    }

    fn cols(&self) -> usize {
        unimplemented!();
    }

    fn eval(&self, _x: usize, _y: usize) -> f64 {
        unimplemented!();
    }
}

//--------------------------------------------------------------------------------------------------
///
/// # Matrix expression for matrix-matrix multiplication.
///
/// This matrix expression represents a matrix-matrix operation that takes
/// two input matrices.
///
/// - tparam T   Element value type.
/// - tparam ME  Matrix expression.
/// - tparam VE  Vector expression.
///
pub struct MatrixMul<E1: MatrixExpression, E2: MatrixExpression> {
    _u: E1,
    _v: E2,
}

impl<E1: MatrixExpression, E2: MatrixExpression> MatrixMul<E1, E2> {
    /// Constructs matrix-matrix multiplication expression for given two input
    /// matrices.
    pub fn new(u: E1, v: E2) -> MatrixMul<E1, E2> {
        return MatrixMul {
            _u: u,
            _v: v,
        };
    }
}

impl<E1: MatrixExpression, E2: MatrixExpression> MatrixExpression for MatrixMul<E1, E2> {
    /// Size of the matrix.
    fn size(&self) -> USize2 {
        return USize2::new(self._u.rows(), self._v.cols());
    }

    /// Number of rows.
    fn rows(&self) -> usize {
        return self._u.rows();
    }

    /// Number of columns.
    fn cols(&self) -> usize {
        return self._v.cols();
    }

    /// Returns matrix element at (i, j).
    fn eval(&self, i: usize, j: usize) -> f64 {
        // todo: Unoptimized mat-mat-mul
        let mut sum: f64 = 0.0;
        let n = self._u.cols();
        for k in 0..n {
            sum += self._u.eval(i, k) * self._v.eval(k, j);
        }
        return sum;
    }
}

//--------------------------------------------------------------------------------------------------
/// # MatrixBinaryOp Aliases

/// Matrix-matrix addition expression.
pub type MatrixAdd<E1, E2> = MatrixBinaryOp<E1, E2, Plus>;
/// Matrix-matrix addition expression.
pub type MatrixScalarAdd<E> = MatrixScalarBinaryOp<E, Plus>;


/// Matrix-matrix subtraction expression.
pub type MatrixSub<E1, E2> = MatrixBinaryOp<E1, E2, Minus>;

/// Matrix-scalar subtraction expression.
pub type MatrixScalarSub<E> = MatrixScalarBinaryOp<E, Minus>;


/// Matrix-matrix subtraction expression with inverse order.
pub type MatrixScalarRSub<E> = MatrixScalarBinaryOp<E, RMinus>;


/// Matrix-scalar multiplication expression.
pub type MatrixScalarMul<E> = MatrixScalarBinaryOp<E, Multiplies>;

/// Matrix-scalar division expression.
pub type MatrixScalarDiv<E> = MatrixScalarBinaryOp<E, Divides>;


/// Matrix-scalar division expression with inverse order.
pub type MatrixScalarRDiv<E> = MatrixScalarBinaryOp<E, RDivides>;
//--------------------------------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------------------------------
/// Returns a matrix with opposite sign.
macro_rules! impl_neg_mat0 {
    ($Matrix:ty) => {
        impl Neg for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn neg(self) -> Self::Output {
                return MatrixScalarMul::new(self, -1.0);
            }
        }
    }
}
impl_neg_mat0!(MatrixConstant);
impl_neg_mat0!(MatrixIdentity);

macro_rules! impl_neg_mat1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression> Neg for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn neg(self) -> Self::Output {
                return MatrixScalarMul::new(self, -1.0);
            }
        }
    }
}
impl_neg_mat1!(MatrixDiagonal<E>);
impl_neg_mat1!(MatrixTriangular<E>);
impl_neg_mat1!(MatrixScalarAdd<E>);
impl_neg_mat1!(MatrixScalarSub<E>);
impl_neg_mat1!(MatrixScalarRSub<E>);
impl_neg_mat1!(MatrixScalarMul<E>);
impl_neg_mat1!(MatrixScalarDiv<E>);
impl_neg_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_neg_mat1_vec1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression> Neg for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn neg(self) -> Self::Output {
                return MatrixScalarMul::new(self, -1.0);
            }
        }
    }
}
impl_neg_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<E1>>);
impl_neg_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<E1>>);
impl_neg_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<E1>>);
impl_neg_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<E1>>);
impl_neg_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<E1>>);
impl_neg_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<E1>>);

macro_rules! impl_neg_mat1_vec2 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression, E2: VectorExpression> Neg for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn neg(self) -> Self::Output {
                return MatrixScalarMul::new(self, -1.0);
            }
        }
    }
}
impl_neg_mat1_vec2!(MatrixVectorMul<E, VectorAdd<E1, E2>>);
impl_neg_mat1_vec2!(MatrixVectorMul<E, VectorSub<E1, E2>>);
impl_neg_mat1_vec2!(MatrixVectorMul<E, VectorMul<E1, E2>>);
impl_neg_mat1_vec2!(MatrixVectorMul<E, VectorDiv<E1, E2>>);

macro_rules! impl_neg_mat2 {
    ($Matrix:ty) => {
        impl<E1: MatrixExpression, E2: MatrixExpression> Neg for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn neg(self) -> Self::Output {
                return MatrixScalarMul::new(self, -1.0);
            }
        }
    }
}
impl_neg_mat2!(MatrixMul<E1, E2>);
impl_neg_mat2!(MatrixAdd<E1, E2>);
impl_neg_mat2!(MatrixSub<E1, E2>);

//--------------------------------------------------------------------------------------------------
/// Returns a + b (element-size).
macro_rules! impl_add_mat0 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression> Add<Return> for $Matrix {
            type Output = MatrixAdd<$Matrix, Return>;

            fn add(self, rhs: Return) -> Self::Output {
                return MatrixAdd::new(self, rhs);
            }
        }
    }
}
impl_add_mat0!(MatrixConstant);
impl_add_mat0!(MatrixIdentity);

macro_rules! impl_add_mat1 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression, E: MatrixExpression> Add<Return> for $Matrix {
            type Output = MatrixAdd<$Matrix, Return>;

            fn add(self, rhs: Return) -> Self::Output {
                return MatrixAdd::new(self, rhs);
            }
        }
    }
}
impl_add_mat1!(MatrixDiagonal<E>);
impl_add_mat1!(MatrixTriangular<E>);
impl_add_mat1!(MatrixScalarAdd<E>);
impl_add_mat1!(MatrixScalarSub<E>);
impl_add_mat1!(MatrixScalarRSub<E>);
impl_add_mat1!(MatrixScalarMul<E>);
impl_add_mat1!(MatrixScalarDiv<E>);
impl_add_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_add_mat1_vec1 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression, E: MatrixExpression, E1: VectorExpression> Add<Return> for $Matrix {
            type Output = MatrixAdd<$Matrix, Return>;

            fn add(self, rhs: Return) -> Self::Output {
                return MatrixAdd::new(self, rhs);
            }
        }
    }
}
impl_add_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<E1>>);
impl_add_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<E1>>);
impl_add_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<E1>>);
impl_add_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<E1>>);
impl_add_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<E1>>);
impl_add_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<E1>>);

macro_rules! impl_add_mat1_vec2 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression, E: MatrixExpression, E1: VectorExpression, E2: VectorExpression> Add<Return> for $Matrix {
            type Output = MatrixAdd<$Matrix, Return>;

            fn add(self, rhs: Return) -> Self::Output {
                return MatrixAdd::new(self, rhs);
            }
        }
    }
}
impl_add_mat1_vec2!(MatrixVectorMul<E, VectorAdd<E1, E2>>);
impl_add_mat1_vec2!(MatrixVectorMul<E, VectorSub<E1, E2>>);
impl_add_mat1_vec2!(MatrixVectorMul<E, VectorMul<E1, E2>>);
impl_add_mat1_vec2!(MatrixVectorMul<E, VectorDiv<E1, E2>>);

macro_rules! impl_add_mat1_vec2 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression, E1: MatrixExpression, E2: MatrixExpression> Add<Return> for $Matrix {
            type Output = MatrixAdd<$Matrix, Return>;

            fn add(self, rhs: Return) -> Self::Output {
                return MatrixAdd::new(self, rhs);
            }
        }
    }
}
impl_add_mat1_vec2!(MatrixMul<E1, E2>);
impl_add_mat1_vec2!(MatrixAdd<E1, E2>);
impl_add_mat1_vec2!(MatrixSub<E1, E2>);

//--------------------------------------------------------------------------------------------------
/// Returns a + b', where every element of matrix b' is b.
macro_rules! impl_add_scalar_mat0 {
    ($Matrix:ty) => {
        impl Add<f64> for $Matrix {
            type Output = MatrixScalarAdd<$Matrix>;

            fn add(self, rhs: f64) -> Self::Output {
                return MatrixScalarAdd::new(self, rhs);
            }
        }
    }
}
impl_add_scalar_mat0!(MatrixConstant);
impl_add_scalar_mat0!(MatrixIdentity);

macro_rules! impl_add_scalar_mat1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression> Add<f64> for $Matrix {
            type Output = MatrixScalarAdd<$Matrix>;

            fn add(self, rhs: f64) -> Self::Output {
                return MatrixScalarAdd::new(self, rhs);
            }
        }
    }
}
impl_add_scalar_mat1!(MatrixDiagonal<E>);
impl_add_scalar_mat1!(MatrixTriangular<E>);
impl_add_scalar_mat1!(MatrixScalarAdd<E>);
impl_add_scalar_mat1!(MatrixScalarSub<E>);
impl_add_scalar_mat1!(MatrixScalarRSub<E>);
impl_add_scalar_mat1!(MatrixScalarMul<E>);
impl_add_scalar_mat1!(MatrixScalarDiv<E>);
impl_add_scalar_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_add_scalar_mat1_vec1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression> Add<f64> for $Matrix {
            type Output = MatrixScalarAdd<$Matrix>;

            fn add(self, rhs: f64) -> Self::Output {
                return MatrixScalarAdd::new(self, rhs);
            }
        }
    }
}
impl_add_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<E1>>);
impl_add_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<E1>>);
impl_add_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<E1>>);
impl_add_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<E1>>);
impl_add_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<E1>>);
impl_add_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<E1>>);

macro_rules! impl_add_scalar_mat1_vec2 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression, E2: VectorExpression> Add<f64> for $Matrix {
            type Output = MatrixScalarAdd<$Matrix>;

            fn add(self, rhs: f64) -> Self::Output {
                return MatrixScalarAdd::new(self, rhs);
            }
        }
    }
}
impl_add_scalar_mat1_vec2!(MatrixVectorMul<E, VectorAdd<E1, E2>>);
impl_add_scalar_mat1_vec2!(MatrixVectorMul<E, VectorSub<E1, E2>>);
impl_add_scalar_mat1_vec2!(MatrixVectorMul<E, VectorMul<E1, E2>>);
impl_add_scalar_mat1_vec2!(MatrixVectorMul<E, VectorDiv<E1, E2>>);

macro_rules! impl_add_scalar_mat2 {
    ($Matrix:ty) => {
        impl<E1: MatrixExpression, E2: MatrixExpression> Add<f64> for $Matrix {
            type Output = MatrixScalarAdd<$Matrix>;

            fn add(self, rhs: f64) -> Self::Output {
                return MatrixScalarAdd::new(self, rhs);
            }
        }
    }
}
impl_add_scalar_mat2!(MatrixMul<E1, E2>);
impl_add_scalar_mat2!(MatrixAdd<E1, E2>);
impl_add_scalar_mat2!(MatrixSub<E1, E2>);

//--------------------------------------------------------------------------------------------------
/// Returns a - b (element-size).
macro_rules! impl_sub_mat0 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression> Sub<Return> for $Matrix {
            type Output = MatrixSub<$Matrix, Return>;

            fn sub(self, rhs: Return) -> Self::Output {
                return MatrixSub::new(self, rhs);
            }
        }
    }
}
impl_sub_mat0!(MatrixConstant);
impl_sub_mat0!(MatrixIdentity);

macro_rules! impl_sub_mat1 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression, E: MatrixExpression> Sub<Return> for $Matrix {
            type Output = MatrixSub<$Matrix, Return>;

            fn sub(self, rhs: Return) -> Self::Output {
                return MatrixSub::new(self, rhs);
            }
        }
    }
}
impl_sub_mat1!(MatrixDiagonal<E>);
impl_sub_mat1!(MatrixTriangular<E>);
impl_sub_mat1!(MatrixScalarAdd<E>);
impl_sub_mat1!(MatrixScalarSub<E>);
impl_sub_mat1!(MatrixScalarRSub<E>);
impl_sub_mat1!(MatrixScalarMul<E>);
impl_sub_mat1!(MatrixScalarDiv<E>);
impl_sub_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_sub_mat1_vec1 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression, E: MatrixExpression, E1: VectorExpression> Sub<Return> for $Matrix {
            type Output = MatrixSub<$Matrix, Return>;

            fn sub(self, rhs: Return) -> Self::Output {
                return MatrixSub::new(self, rhs);
            }
        }
    }
}
impl_sub_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<E1>>);
impl_sub_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<E1>>);
impl_sub_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<E1>>);
impl_sub_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<E1>>);
impl_sub_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<E1>>);
impl_sub_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<E1>>);

macro_rules! impl_sub_mat1_vec2 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression, E: MatrixExpression, E1: VectorExpression, E2: VectorExpression> Sub<Return> for $Matrix {
            type Output = MatrixSub<$Matrix, Return>;

            fn sub(self, rhs: Return) -> Self::Output {
                return MatrixSub::new(self, rhs);
            }
        }
    }
}
impl_sub_mat1_vec2!(MatrixVectorMul<E, VectorAdd<E1, E2>>);
impl_sub_mat1_vec2!(MatrixVectorMul<E, VectorSub<E1, E2>>);
impl_sub_mat1_vec2!(MatrixVectorMul<E, VectorMul<E1, E2>>);
impl_sub_mat1_vec2!(MatrixVectorMul<E, VectorDiv<E1, E2>>);

macro_rules! impl_sub_mat2 {
    ($Matrix:ty) => {
        impl<Return: MatrixExpression, E1: MatrixExpression, E2: MatrixExpression> Sub<Return> for $Matrix {
            type Output = MatrixSub<$Matrix, Return>;

            fn sub(self, rhs: Return) -> Self::Output {
                return MatrixSub::new(self, rhs);
            }
        }
    }
}
impl_sub_mat2!(MatrixMul<E1, E2>);
impl_sub_mat2!(MatrixAdd<E1, E2>);
impl_sub_mat2!(MatrixSub<E1, E2>);

//--------------------------------------------------------------------------------------------------
/// Returns a - b', where every element of matrix b' is b.
macro_rules! impl_sub_scalar_mat0 {
    ($Matrix:ty) => {
        impl Sub<f64> for $Matrix {
            type Output = MatrixScalarSub<$Matrix>;

            fn sub(self, rhs: f64) -> Self::Output {
                return MatrixScalarSub::new(self, rhs);
            }
        }
    }
}
impl_sub_scalar_mat0!(MatrixConstant);
impl_sub_scalar_mat0!(MatrixIdentity);

macro_rules! impl_sub_scalar_mat1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression> Sub<f64> for $Matrix {
            type Output = MatrixScalarSub<$Matrix>;

            fn sub(self, rhs: f64) -> Self::Output {
                return MatrixScalarSub::new(self, rhs);
            }
        }
    }
}
impl_sub_scalar_mat1!(MatrixDiagonal<E>);
impl_sub_scalar_mat1!(MatrixTriangular<E>);
impl_sub_scalar_mat1!(MatrixScalarAdd<E>);
impl_sub_scalar_mat1!(MatrixScalarSub<E>);
impl_sub_scalar_mat1!(MatrixScalarRSub<E>);
impl_sub_scalar_mat1!(MatrixScalarMul<E>);
impl_sub_scalar_mat1!(MatrixScalarDiv<E>);
impl_sub_scalar_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_sub_scalar_mat1_vec1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression> Sub<f64> for $Matrix {
            type Output = MatrixScalarSub<$Matrix>;

            fn sub(self, rhs: f64) -> Self::Output {
                return MatrixScalarSub::new(self, rhs);
            }
        }
    }
}
impl_sub_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<E1>>);
impl_sub_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<E1>>);
impl_sub_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<E1>>);
impl_sub_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<E1>>);
impl_sub_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<E1>>);
impl_sub_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<E1>>);

macro_rules! impl_sub_scalar_mat1_vec2 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression, E2: VectorExpression> Sub<f64> for $Matrix {
            type Output = MatrixScalarSub<$Matrix>;

            fn sub(self, rhs: f64) -> Self::Output {
                return MatrixScalarSub::new(self, rhs);
            }
        }
    }
}
impl_sub_scalar_mat1_vec2!(MatrixVectorMul<E, VectorAdd<E1, E2>>);
impl_sub_scalar_mat1_vec2!(MatrixVectorMul<E, VectorSub<E1, E2>>);
impl_sub_scalar_mat1_vec2!(MatrixVectorMul<E, VectorMul<E1, E2>>);
impl_sub_scalar_mat1_vec2!(MatrixVectorMul<E, VectorDiv<E1, E2>>);

macro_rules! impl_sub_scalar_mat2 {
    ($Matrix:ty) => {
        impl<E1: MatrixExpression, E2: MatrixExpression> Sub<f64> for $Matrix {
            type Output = MatrixScalarSub<$Matrix>;

            fn sub(self, rhs: f64) -> Self::Output {
                return MatrixScalarSub::new(self, rhs);
            }
        }
    }
}
impl_sub_scalar_mat2!(MatrixMul<E1, E2>);
impl_sub_scalar_mat2!(MatrixAdd<E1, E2>);
impl_sub_scalar_mat2!(MatrixSub<E1, E2>);

//--------------------------------------------------------------------------------------------------
/// Returns a * b', where every element of matrix b' is b.
macro_rules! impl_mul_scalar_mat0 {
    ($Matrix:ty) => {
        impl Mul<f64> for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn mul(self, rhs: f64) -> Self::Output {
                return MatrixScalarMul::new(self, rhs);
            }
        }
    }
}
impl_mul_scalar_mat0!(MatrixConstant);
impl_mul_scalar_mat0!(MatrixIdentity);

macro_rules! impl_mul_scalar_mat1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression> Mul<f64> for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn mul(self, rhs: f64) -> Self::Output {
                return MatrixScalarMul::new(self, rhs);
            }
        }
    }
}
impl_mul_scalar_mat1!(MatrixDiagonal<E>);
impl_mul_scalar_mat1!(MatrixTriangular<E>);
impl_mul_scalar_mat1!(MatrixScalarAdd<E>);
impl_mul_scalar_mat1!(MatrixScalarSub<E>);
impl_mul_scalar_mat1!(MatrixScalarRSub<E>);
impl_mul_scalar_mat1!(MatrixScalarMul<E>);
impl_mul_scalar_mat1!(MatrixScalarDiv<E>);
impl_mul_scalar_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_mul_scalar_mat1_vec1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression> Mul<f64> for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn mul(self, rhs: f64) -> Self::Output {
                return MatrixScalarMul::new(self, rhs);
            }
        }
    }
}
impl_mul_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<E1>>);
impl_mul_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<E1>>);
impl_mul_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<E1>>);
impl_mul_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<E1>>);
impl_mul_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<E1>>);
impl_mul_scalar_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<E1>>);

macro_rules! impl_mul_scalar_mat1_vec2 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression, E2: VectorExpression> Mul<f64> for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn mul(self, rhs: f64) -> Self::Output {
                return MatrixScalarMul::new(self, rhs);
            }
        }
    }
}
impl_mul_scalar_mat1_vec2!(MatrixVectorMul<E, VectorAdd<E1, E2>>);
impl_mul_scalar_mat1_vec2!(MatrixVectorMul<E, VectorSub<E1, E2>>);
impl_mul_scalar_mat1_vec2!(MatrixVectorMul<E, VectorMul<E1, E2>>);
impl_mul_scalar_mat1_vec2!(MatrixVectorMul<E, VectorDiv<E1, E2>>);

macro_rules! impl_mul_scalar_mat2 {
    ($Matrix:ty) => {
        impl<E1: MatrixExpression, E2: MatrixExpression> Mul<f64> for $Matrix {
            type Output = MatrixScalarMul<$Matrix>;

            fn mul(self, rhs: f64) -> Self::Output {
                return MatrixScalarMul::new(self, rhs);
            }
        }
    }
}
impl_mul_scalar_mat2!(MatrixMul<E1, E2>);
impl_mul_scalar_mat2!(MatrixAdd<E1, E2>);
impl_mul_scalar_mat2!(MatrixSub<E1, E2>);

//--------------------------------------------------------------------------------------------------
/// Returns a * b.
macro_rules! impl_mul_vec_mat0 {
    ($Matrix:ty) => {
        macro_rules! impl_mul_vec1_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec1_mat0!(VectorScalarAdd<E1>);
        impl_mul_vec1_mat0!(VectorScalarSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarRSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarMul<E1>);
        impl_mul_vec1_mat0!(VectorScalarDiv<E1>);
        impl_mul_vec1_mat0!(VectorScalarRDiv<E1>);

        macro_rules! impl_mul_vec2_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, E2:VectorExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec2_mat0!(VectorAdd<E1, E2>);
        impl_mul_vec2_mat0!(VectorSub<E1, E2>);
        impl_mul_vec2_mat0!(VectorMul<E1, E2>);
        impl_mul_vec2_mat0!(VectorDiv<E1, E2>);
    }
}
impl_mul_vec_mat0!(MatrixConstant);
impl_mul_vec_mat0!(MatrixIdentity);

macro_rules! impl_mul_vec_mat1 {
    ($Matrix:ty) => {
        macro_rules! impl_mul_vec1_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, E: MatrixExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec1_mat0!(VectorScalarAdd<E1>);
        impl_mul_vec1_mat0!(VectorScalarSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarRSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarMul<E1>);
        impl_mul_vec1_mat0!(VectorScalarDiv<E1>);
        impl_mul_vec1_mat0!(VectorScalarRDiv<E1>);

        macro_rules! impl_mul_vec2_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, E2:VectorExpression, E: MatrixExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec2_mat0!(VectorAdd<E1, E2>);
        impl_mul_vec2_mat0!(VectorSub<E1, E2>);
        impl_mul_vec2_mat0!(VectorMul<E1, E2>);
        impl_mul_vec2_mat0!(VectorDiv<E1, E2>);
    }
}
impl_mul_vec_mat1!(MatrixDiagonal<E>);
impl_mul_vec_mat1!(MatrixTriangular<E>);
impl_mul_vec_mat1!(MatrixScalarAdd<E>);
impl_mul_vec_mat1!(MatrixScalarSub<E>);
impl_mul_vec_mat1!(MatrixScalarRSub<E>);
impl_mul_vec_mat1!(MatrixScalarMul<E>);
impl_mul_vec_mat1!(MatrixScalarDiv<E>);
impl_mul_vec_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_mul_vec_mat1_vec1 {
    ($Matrix:ty) => {
        macro_rules! impl_mul_vec1_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, E: MatrixExpression, V1: VectorExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec1_mat0!(VectorScalarAdd<E1>);
        impl_mul_vec1_mat0!(VectorScalarSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarRSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarMul<E1>);
        impl_mul_vec1_mat0!(VectorScalarDiv<E1>);
        impl_mul_vec1_mat0!(VectorScalarRDiv<E1>);

        macro_rules! impl_mul_vec2_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, E2:VectorExpression, E: MatrixExpression, V1: VectorExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec2_mat0!(VectorAdd<E1, E2>);
        impl_mul_vec2_mat0!(VectorSub<E1, E2>);
        impl_mul_vec2_mat0!(VectorMul<E1, E2>);
        impl_mul_vec2_mat0!(VectorDiv<E1, E2>);
    }
}
impl_mul_vec_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<V1>>);
impl_mul_vec_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<V1>>);
impl_mul_vec_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<V1>>);
impl_mul_vec_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<V1>>);
impl_mul_vec_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<V1>>);
impl_mul_vec_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<V1>>);

macro_rules! impl_mul_vec_mat1_vec2 {
    ($Matrix:ty) => {
        macro_rules! impl_mul_vec1_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, E: MatrixExpression, V1: VectorExpression, V2: VectorExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec1_mat0!(VectorScalarAdd<E1>);
        impl_mul_vec1_mat0!(VectorScalarSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarRSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarMul<E1>);
        impl_mul_vec1_mat0!(VectorScalarDiv<E1>);
        impl_mul_vec1_mat0!(VectorScalarRDiv<E1>);

        macro_rules! impl_mul_vec2_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, E2:VectorExpression, E: MatrixExpression, V1: VectorExpression, V2: VectorExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec2_mat0!(VectorAdd<E1, E2>);
        impl_mul_vec2_mat0!(VectorSub<E1, E2>);
        impl_mul_vec2_mat0!(VectorMul<E1, E2>);
        impl_mul_vec2_mat0!(VectorDiv<E1, E2>);
    }
}
impl_mul_vec_mat1_vec2!(MatrixVectorMul<E, VectorAdd<V1, V2>>);
impl_mul_vec_mat1_vec2!(MatrixVectorMul<E, VectorSub<V1, V2>>);
impl_mul_vec_mat1_vec2!(MatrixVectorMul<E, VectorMul<V1, V2>>);
impl_mul_vec_mat1_vec2!(MatrixVectorMul<E, VectorDiv<V1, V2>>);

macro_rules! impl_mul_vec_mat2 {
    ($Matrix:ty) => {
        macro_rules! impl_mul_vec1_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, V1: MatrixExpression, V2: MatrixExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec1_mat0!(VectorScalarAdd<E1>);
        impl_mul_vec1_mat0!(VectorScalarSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarRSub<E1>);
        impl_mul_vec1_mat0!(VectorScalarMul<E1>);
        impl_mul_vec1_mat0!(VectorScalarDiv<E1>);
        impl_mul_vec1_mat0!(VectorScalarRDiv<E1>);

        macro_rules! impl_mul_vec2_mat0 {
            ($Vector:ty) => {
                impl<E1:VectorExpression, E2:VectorExpression, V1: MatrixExpression, V2: MatrixExpression> Mul<$Vector> for $Matrix {
                    type Output = MatrixVectorMul<$Matrix, $Vector>;

                    fn mul(self, rhs: $Vector) -> Self::Output {
                        return MatrixVectorMul::new(self, rhs);
                    }
                }
            };
        }
        impl_mul_vec2_mat0!(VectorAdd<E1, E2>);
        impl_mul_vec2_mat0!(VectorSub<E1, E2>);
        impl_mul_vec2_mat0!(VectorMul<E1, E2>);
        impl_mul_vec2_mat0!(VectorDiv<E1, E2>);
    }
}
impl_mul_vec_mat2!(MatrixMul<V1, V2>);
impl_mul_vec_mat2!(MatrixAdd<V1, V2>);
impl_mul_vec_mat2!(MatrixSub<V1, V2>);

//--------------------------------------------------------------------------------------------------
/// Returns a * b.
macro_rules! impl_mul_mat0 {
    ($Matrix:ty) => {
        impl<Matrix: MatrixExpression> Mul<Matrix> for $Matrix {
            type Output = MatrixMul<$Matrix, Matrix>;

            fn mul(self, rhs: Matrix) -> Self::Output {
                return MatrixMul::new(self, rhs);
            }
        }
    }
}
impl_mul_mat0!(MatrixConstant);
impl_mul_mat0!(MatrixIdentity);

macro_rules! impl_mul_mat1 {
    ($Matrix:ty) => {
        impl<Matrix: MatrixExpression, E: MatrixExpression> Mul<Matrix> for $Matrix {
            type Output = MatrixMul<$Matrix, Matrix>;

            fn mul(self, rhs: Matrix) -> Self::Output {
                return MatrixMul::new(self, rhs);
            }
        }
    }
}
impl_mul_mat1!(MatrixDiagonal<E>);
impl_mul_mat1!(MatrixTriangular<E>);
impl_mul_mat1!(MatrixScalarAdd<E>);
impl_mul_mat1!(MatrixScalarSub<E>);
impl_mul_mat1!(MatrixScalarRSub<E>);
impl_mul_mat1!(MatrixScalarMul<E>);
impl_mul_mat1!(MatrixScalarDiv<E>);
impl_mul_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_mul_mat1_vec1 {
    ($Matrix:ty) => {
        impl<Matrix: MatrixExpression, E: MatrixExpression, E1: VectorExpression> Mul<Matrix> for $Matrix {
            type Output = MatrixMul<$Matrix, Matrix>;

            fn mul(self, rhs: Matrix) -> Self::Output {
                return MatrixMul::new(self, rhs);
            }
        }
    }
}
impl_mul_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<E1>>);
impl_mul_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<E1>>);
impl_mul_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<E1>>);
impl_mul_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<E1>>);
impl_mul_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<E1>>);
impl_mul_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<E1>>);

macro_rules! impl_mul_mat1_vec2 {
    ($Matrix:ty) => {
        impl<Matrix: MatrixExpression, E: MatrixExpression, E1: VectorExpression, E2: VectorExpression> Mul<Matrix> for $Matrix {
            type Output = MatrixMul<$Matrix, Matrix>;

            fn mul(self, rhs: Matrix) -> Self::Output {
                return MatrixMul::new(self, rhs);
            }
        }
    }
}
impl_mul_mat1_vec2!(MatrixVectorMul<E, VectorAdd<E1, E2>>);
impl_mul_mat1_vec2!(MatrixVectorMul<E, VectorSub<E1, E2>>);
impl_mul_mat1_vec2!(MatrixVectorMul<E, VectorMul<E1, E2>>);
impl_mul_mat1_vec2!(MatrixVectorMul<E, VectorDiv<E1, E2>>);

macro_rules! impl_mul_mat2 {
    ($Matrix:ty) => {
        impl<Matrix: MatrixExpression, E1: MatrixExpression, E2: MatrixExpression> Mul<Matrix> for $Matrix {
            type Output = MatrixMul<$Matrix, Matrix>;

            fn mul(self, rhs: Matrix) -> Self::Output {
                return MatrixMul::new(self, rhs);
            }
        }
    }
}
impl_mul_mat2!(MatrixMul<E1, E2>);
impl_mul_mat2!(MatrixAdd<E1, E2>);
impl_mul_mat2!(MatrixSub<E1, E2>);

//--------------------------------------------------------------------------------------------------
/// Returns a' / b, where every element of matrix a' is a.
macro_rules! impl_div_mat0 {
    ($Matrix:ty) => {
        impl Div<f64> for $Matrix {
            type Output = MatrixScalarDiv<$Matrix>;

            fn div(self, rhs: f64) -> Self::Output {
                return MatrixScalarDiv::new(self, rhs);
            }
        }
    }
}

impl_div_mat0!(MatrixConstant);
impl_div_mat0!(MatrixIdentity);

macro_rules! impl_div_mat1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression> Div<f64> for $Matrix {
            type Output = MatrixScalarDiv<$Matrix>;

            fn div(self, rhs: f64) -> Self::Output {
                return MatrixScalarDiv::new(self, rhs);
            }
        }
    }
}
impl_div_mat1!(MatrixDiagonal<E>);
impl_div_mat1!(MatrixTriangular<E>);
impl_div_mat1!(MatrixScalarAdd<E>);
impl_div_mat1!(MatrixScalarSub<E>);
impl_div_mat1!(MatrixScalarRSub<E>);
impl_div_mat1!(MatrixScalarMul<E>);
impl_div_mat1!(MatrixScalarDiv<E>);
impl_div_mat1!(MatrixScalarRDiv<E>);

macro_rules! impl_div_mat1_vec1 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression> Div<f64> for $Matrix {
            type Output = MatrixScalarDiv<$Matrix>;

            fn div(self, rhs: f64) -> Self::Output {
                return MatrixScalarDiv::new(self, rhs);
            }
        }
    }
}
impl_div_mat1_vec1!(MatrixVectorMul<E, VectorScalarAdd<E1>>);
impl_div_mat1_vec1!(MatrixVectorMul<E, VectorScalarSub<E1>>);
impl_div_mat1_vec1!(MatrixVectorMul<E, VectorScalarRSub<E1>>);
impl_div_mat1_vec1!(MatrixVectorMul<E, VectorScalarMul<E1>>);
impl_div_mat1_vec1!(MatrixVectorMul<E, VectorScalarDiv<E1>>);
impl_div_mat1_vec1!(MatrixVectorMul<E, VectorScalarRDiv<E1>>);

macro_rules! impl_div_mat1_vec2 {
    ($Matrix:ty) => {
        impl<E: MatrixExpression, E1: VectorExpression, E2: VectorExpression> Div<f64> for $Matrix {
            type Output = MatrixScalarDiv<$Matrix>;

            fn div(self, rhs: f64) -> Self::Output {
                return MatrixScalarDiv::new(self, rhs);
            }
        }
    }
}
impl_div_mat1_vec2!(MatrixVectorMul<E, VectorAdd<E1, E2>>);
impl_div_mat1_vec2!(MatrixVectorMul<E, VectorSub<E1, E2>>);
impl_div_mat1_vec2!(MatrixVectorMul<E, VectorMul<E1, E2>>);
impl_div_mat1_vec2!(MatrixVectorMul<E, VectorDiv<E1, E2>>);

macro_rules! impl_div_mat2 {
    ($Matrix:ty) => {
        impl<E1: MatrixExpression, E2: MatrixExpression> Div<f64> for $Matrix {
            type Output = MatrixScalarDiv<$Matrix>;

            fn div(self, rhs: f64) -> Self::Output {
                return MatrixScalarDiv::new(self, rhs);
            }
        }
    }
}
impl_div_mat2!(MatrixMul<E1, E2>);
impl_div_mat2!(MatrixAdd<E1, E2>);
impl_div_mat2!(MatrixSub<E1, E2>);