wmathrs 0.1.0

use crate::scalar::{ One, Zero, Latex, Similar, Neg, Con, Pow, DivLeft }; 
use std::fmt::{ Formatter, Display, Result }; 
use std::cmp::{ PartialEq, Eq }; 
use std::ops::{ Add, Sub, Mul, Div }; 

#[derive(Debug, Clone)]
/// 矩阵
pub struct Matrix<T> {
    kernel: Vec<T>, 
    row: usize, 
    col: usize, 
}

#[macro_export]
/// mat!(k, r, c) => Matrix::new(k, r, c)
macro_rules! mat {
    ($kernel:expr, $row:expr, $col:expr) => {
        Matrix::new($kernel, $row, $col)
    }; 
}

impl <T: Clone> Matrix<T> {
    pub fn new(kernel: &Vec<T>, row: usize, col: usize) -> Self {
        assert!(row > 0 && col > 0); 
        assert!(kernel.len() == row * col); 
        Matrix {
            kernel: kernel.clone(), 
            row: row, 
            col: col, 
        }
    }

    pub fn get_kernel(&self) -> &Vec<T> {
        &self.kernel
    }

    /// [危险函数，外部操作可能会改变kernel长度]
    pub fn get_kernel_mut(&mut self) -> &mut Vec<T> {
        &mut self.kernel
    }

    pub fn get_row(&self) -> &usize {
        &self.row
    }

    pub fn get_col(&self) -> &usize {
        &self.col
    }
}

impl <T: Display> Display for Matrix<T> {
    fn fmt(&self, f: &mut Formatter) -> Result {
        write!(f, "Matrix [[")?; 
        for (i, item) in self.kernel.iter().enumerate() {
            if i % self.col == self.col - 1 {
                if i == self.kernel.len() - 1 { write!(f, "{}]]", item)?; }
                else { write!(f, "{}], \n        [", item)?; }
            } else {
                write!(f, "{}, ", item)?; 
            }
        }
        Ok(())
    }
}

impl <T: Latex> Latex for Matrix<T> {
    fn latex(&self) -> String {
        let mut string = String::from("\\begin{pmatrix}\n"); 
        for (i, item) in self.kernel.iter().enumerate() {
            if i % self.col == self.col - 1 {
                string.push_str(&format!("{{{}}}\\\\\n", item.latex())); 
            } else {
                string.push_str(&format!("{{{}}}&", item.latex())); 
            }
        }
        string.push_str("\\end{pmatrix}"); 
        return string; 
    }
}

impl <T: Clone + One + Zero + PartialEq> One for Matrix<T> {
    fn get_one(&self) -> Self {
        let mut kernel = Vec::with_capacity(self.kernel.len()); 
        for (i, item) in self.kernel.iter().enumerate() {
            if i / self.col == i % self.col {
                kernel.push(item.get_one()); 
            } else {
                kernel.push(item.get_zero()); 
            }
        }
        Matrix::new(&kernel, self.row, self.col)
    }

    fn eq_one(&self) -> bool {
        for (i, item) in self.kernel.iter().enumerate() {
            if i / self.col == i % self.col {
                if !item.eq_one() { return false; }
            } else {
                if !item.eq_zero() { return false; }
            }
        }
        return true; 
    }

    fn similar_one(&self, precision: f64) -> bool {
        for (i, item) in self.kernel.iter().enumerate() {
            if i / self.col == i % self.col {
                if !item.similar_one(precision) { return false; }
            } else {
                if !item.similar_zero(precision) { return false; }
            }
        }
        return true; 
    }
}

impl <T: Clone + Zero + PartialEq> Zero for Matrix<T> {
    fn get_zero(&self) -> Self {
        let kernel = self.kernel.iter().map(|x| x.get_zero()).collect::<Vec<T>>(); 
        Matrix::new(&kernel, self.row, self.col)
    }

    fn eq_zero(&self) -> bool {
        !self.kernel.iter().any(|x| !x.eq_zero())
    }

    fn similar_zero(&self, precision: f64) -> bool {
        !self.kernel.iter().any(|x| !x.similar_zero(precision))
    }
}

impl <T: Neg> Neg for Matrix<T> where <T as Neg>::Output: Clone {
    type Output = Matrix<<T as Neg>::Output>; 
    fn neg(&self) -> Self::Output {
        let kernel = self.kernel.iter().map(|x| x.neg()).collect(); 
        Matrix::new(&kernel, self.row, self.col)
    }
}

impl <T: Con> Con for Matrix<T> where <T as Con>::Output: Clone {
    type Output = Matrix<<T as Con>::Output>; 
    fn con(&self) -> Self::Output {
        let kernel = self.kernel.iter().map(|x| x.con()).collect(); 
        Matrix::new(&kernel, self.row, self.col)
    }
}

impl <T: Clone + Mul<Output = T> + Add<Output = T> + One + Zero + PartialEq> Pow for Matrix<T> {
    type Output = Self; 
    fn pow(&self, power: u32) -> Self::Output {
        assert!(self.row == self.col); 
        let mut result = self.get_one(); 
        let mut now = self.clone(); 
        let mut power = power; 
        while power != 0 {
            if power & 1 == 1 {
                result = result.clone() * now.clone(); 
            }
            now = now.clone() * now.clone(); 
            power >>= 1; 
        }
        return result; 
    }
}

impl <T: PartialEq> PartialEq for Matrix<T> {
    fn eq(&self, other: &Self) -> bool {
        if (self.row != other.row) || (self.col != other.col) { return false; }
        !self.kernel.iter().zip(other.kernel.iter()).any(|(x, y)| x != y)
    }
}

impl <T: Eq> Eq for Matrix<T> {}

impl <T: Similar> Similar for Matrix<T> {
    fn similar(&self, other: &Self, precision: f64) -> bool {
        if (self.row != other.row) || (self.col != other.col) { return false; }
        !self.kernel.iter().zip(other.kernel.iter()).any(|(x, y)| !x.similar(y, precision))
    }
}

impl <T, U, V> Add<Matrix<U>> for Matrix<T> where
    T: Clone + Add<U, Output = V>, 
    U: Clone, 
    V: Clone {
    type Output = Matrix<V>; 
    fn add(self, other: Matrix<U>) -> Self::Output {
        assert!(self.row == other.row && self.col == other.col); 
        let kernel = self.kernel.iter().zip(other.kernel.iter()).map(|(x, y)| x.clone() + y.clone()).collect(); 
        Matrix::new(&kernel, self.row, self.col)
    }
}

impl <T, U, V> Sub<Matrix<U>> for Matrix<T> where
    T: Clone + Sub<U, Output = V>, 
    U: Clone, 
    V: Clone {
    type Output = Matrix<V>; 
    fn sub(self, other: Matrix<U>) -> Self::Output {
        assert!(self.row == other.row && self.col == other.col); 
        let kernel = self.kernel.iter().zip(other.kernel.iter()).map(|(x, y)| x.clone() - y.clone()).collect(); 
        Matrix::new(&kernel, self.row, self.col)
    }
}

impl <T, U, V> Mul<Matrix<U>> for Matrix<T> where
    T: Clone + Mul<U, Output = V>, 
    U: Clone, 
    V: Clone + Add<Output = V> {
    type Output = Matrix<V>; 
    fn mul(self, other: Matrix<U>) -> Self::Output {
        assert!(self.col == other.row); 
        let mut kernel = Vec::with_capacity(self.row * other.col); 
        for i in 0..self.row {
            for j in 0..other.col {
                // kernel.push(self.kernel[(i * self.col)..((i+1) * self.col)].iter().zip(other.kernel.iter().skip(j).step_by(other.row)).map(|(x, y)| x.clone() * y.clone()).sum::<V>()); 
                let mut sum = self.kernel[i * self.col].clone() * other.kernel[j].clone(); 
                for k in 1..self.col {
                    sum = sum.clone() + self.kernel[i * self.col + k].clone() * other.kernel[k * other.col + j].clone(); 
                }
                kernel.push(sum); 
            }
        }
        Matrix::new(&kernel, self.row, other.col)
    }
}

impl <T> Div for Matrix<T> where
    T: Clone + One + Zero + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + PartialEq {
    type Output = Self; 
    fn div(self, other: Self) -> Self::Output {
        if let Some(inv) = other.inverse() {
            return self * inv; 
        } else { panic!("被除数不可逆！"); }
    }
}

impl <T> Matrix<T> where
    T: Clone + One + Zero + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + PartialEq {
    pub fn try_div(self, other: Self) -> Option<Self> {
        if let Some(inv) = other.inverse() {
            return Some(self * inv); 
        } else { return None; }
    }
}

impl <T> DivLeft for Matrix<T> where
    T: Clone + One + Zero + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + PartialEq {
    type Output = Self; 
    fn div_left(&self, other: &Self) -> Self::Output {
        if let Some(inv) = other.inverse_left() {
            return inv * self.clone(); 
        } else { panic!("被除数不可逆！"); }
    }
}

impl <T> Matrix<T> where
    T: Clone + One + Zero + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + PartialEq {
    pub fn try_div_left(&self, other: &Self) -> Option<Self> {
        if let Some(inv) = other.inverse_right() {
            return Some(inv * self.clone()); 
        } else { return None; }
    }
}

impl <T: Clone> Matrix<T> {
    /// 转置
    pub fn transpose(&self) -> Self {
        let mut kernel = Vec::with_capacity(self.kernel.len()); 
        for r in 0..self.col {
            for c in 0..self.row {
                kernel.push(self.kernel[c * self.col + r].clone()); 
            }
        }
        return Matrix::new(&kernel, self.col, self.row); 
    }

    /// 行交换
    pub fn trans_rows(&self, r1: usize, r2: usize) -> Self {
        let mut kernel = self.kernel.clone(); 
        for c in 0..self.col {
            kernel[r1 * self.col + c] = self.kernel[r2 * self.col + c].clone(); 
            kernel[r2 * self.col + c] = self.kernel[r1 * self.col + c].clone(); 
        }
        return Matrix::new(&kernel, self.row, self.col); 
    }

    /// 列交换
    pub fn trans_cols(&self, c1: usize, c2: usize) -> Self {
        let mut kernel = self.kernel.clone(); 
        for r in 0..self.row {
            kernel[r * self.col + c1] = self.kernel[r * self.col + c2].clone(); 
            kernel[r * self.col + c2] = self.kernel[r * self.col + c1].clone(); 
        }
        return Matrix::new(&kernel, self.row, self.col); 
    }
}

impl <T> Matrix<T> where
    T: Clone + One + Zero + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T> + PartialEq {
    /// 秩
    pub fn rank(&self) -> usize {
        let (_, independent_cols, _) = self.step(false, false); 
        return independent_cols.len(); 
    }
    
    /// 阶梯化<br/>
    /// standardize: 是否进行标准型<br/>
    /// simplify: 是否进行简化标准型<br/>
    /// 返回值中的(neg: bool)为ture则说明step过程中改变了矩阵行列式的符号
    pub fn step(&self, standardize: bool, simplify: bool) -> (Self, Vec<usize>, bool) {
        let mut _self = self.clone(); 
        let mut independent_cols: Vec<usize> = vec![]; 
        let mut neg = false; 
        let (mut row, mut col) = (0, 0); 
        while (row < _self.row) && (col < _self.col) {
            // skip rows beginning with 'zero'
            let mut _row = row; 
            while (_row < _self.row) && _self.kernel[_row * _self.col + col].eq_zero() { _row += 1; }
            if _row == _self.row {
                col += 1; 
                continue; 
            }
            independent_cols.push(col); 
            // upper triangle method
            for r in (_row+1).._self.row {
                let times = _self.kernel[r * _self.col + col].clone() / _self.kernel[_row * _self.col + col].clone(); 
                _self.kernel[r * _self.col + col] = _self.kernel[r * _self.col + col].get_zero(); 
                for c in (col+1).._self.col {
                    _self.kernel[r * _self.col + c] = _self.kernel[r * _self.col + c].clone() - times.clone() * _self.kernel[_row * _self.col + c].clone(); 
                }
            }
            // switch position if necesary
            if _row != row {
                for c in 0.._self.col {
                    let temp = _self.kernel[row * _self.col + c].clone(); 
                    _self.kernel[row * _self.col + c] = _self.kernel[_row * _self.col + c].clone(); 
                    _self.kernel[_row * _self.col + c] = temp; 
                }
                neg = !neg; 
            }
            row += 1; 
            col += 1; 
        }
        if simplify {
            let mut dependent_cols: Vec<usize> = vec![]; 
            for c in 0.._self.col {
                if !independent_cols.contains(&c) { dependent_cols.push(c); }
            }
            for i in 0..independent_cols.len() {
                let ic = independent_cols.len() - 1 - i; 
                for r in 0..ic {
                    let times = _self.kernel[r * _self.col + independent_cols[ic]].clone() / _self.kernel[ic * _self.col + independent_cols[ic]].clone(); 
                    for dc in dependent_cols.iter() {
                        if dc > &independent_cols[ic] {
                            _self.kernel[r * _self.col + dc] = _self.kernel[r * _self.col + dc].clone() - times.clone() * _self.kernel[ic * _self.col + dc].clone(); 
                        }
                    }
                    _self.kernel[r * _self.col + independent_cols[ic]] = _self.kernel[r * _self.col + independent_cols[ic]].get_zero(); 
                }
                for dc in dependent_cols.iter() {
                    if dc > &independent_cols[ic] {
                        _self.kernel[ic * _self.col + dc] = _self.kernel[ic * _self.col + dc].clone() / _self.kernel[ic * _self.col + independent_cols[ic]].clone(); 
                    }
                }
                _self.kernel[ic * _self.col + independent_cols[ic]] = _self.kernel[ic * _self.col + independent_cols[ic]].get_one(); 
            }
        } else if standardize {
            for i in 0..independent_cols.len() {
                let ic = independent_cols.len() - 1 - i; 
                for c in (independent_cols[ic]+1).._self.col {
                    _self.kernel[ic * _self.col + c] = _self.kernel[ic * _self.col + c].clone() / _self.kernel[ic * _self.col + independent_cols[ic]].clone(); 
                }
                _self.kernel[ic * _self.col + independent_cols[ic]] = _self.kernel[ic * _self.col + independent_cols[ic]].get_one(); 
            }
        }
        return (_self, independent_cols.clone(), neg); 
    }
    
    /// 非奇异矩阵的逆
    pub fn inverse(&self) -> Option<Self> {
        if self.row != self.col { return None; }
        let mut kernel = self.kernel.clone(); 
        // create kernel as result
        let mut inverse = Vec::with_capacity(kernel.len()); 
        let one = kernel[0].get_one(); 
        let zero = kernel[0].get_zero(); 
        for i in 0..kernel.len() {
            if i / self.col == i % self.col {
                inverse.push(one.clone()); 
            } else {
                inverse.push(zero.clone()); 
            }
        }
        // upper triangle
        for col in 0..self.col {
            // skip rows beginning with a 'zero'
            let mut row = col; 
            while row < self.row && kernel[row * self.col + col].eq_zero() { row += 1; }
            if row == self.col { return None; }

            for _row in (row+1)..self.row {
                let times = kernel[_row * self.col + col].clone() / kernel[row * self.col + col].clone(); 
                for _col in (col+1)..self.col {
                    kernel[_row * self.col + _col] = kernel[_row * self.col + _col].clone() - times.clone() * kernel[row * self.col + _col].clone(); 
                }
                for _col in 0..self.col {
                    inverse[_row * self.col + _col] = inverse[_row * self.col + _col].clone() - times.clone() * inverse[row * self.col + _col].clone(); 
                }
            }
            // switch position if necessary
            if row != col {
                for i in 0..self.col {
                    let temp = kernel[row * self.col + i].clone(); 
                    kernel[row * self.col + i] = kernel[col * self.col + i].clone(); 
                    kernel[col * self.col + i] = temp; 
                    let temp = inverse[row * self.col + i].clone(); 
                    inverse[row * self.col + i] = inverse[col * self.col + i].clone(); 
                    inverse[col * self.col + i] = temp; 
                }
            }
        }
        // lower triangle
        for col_inv in 0..self.col {
            let col = self.col - 1 - col_inv; 
            for _row in 0..col {
                let times = kernel[_row * self.col + col].clone() / kernel[col * self.col + col].clone(); 
                for _col in 0..self.col {
                    inverse[_row * self.col + _col] = inverse[_row * self.col + _col].clone() - times.clone() * inverse[col * self.col + _col].clone(); 
                }
            }
            for _col in 0..self.col {
                inverse[col * self.col + _col] = inverse[col * self.col + _col].clone() / kernel[col * self.col + col].clone(); 
            }
        }

        return Some(Matrix::new(&inverse, self.row, self.col)); 
    }

    /// 矩阵的左逆
    pub fn inverse_left(&self) -> Option<Self> {
        if self.row < self.col { return None; }
        if self.row == self.col { return self.inverse(); }
        if let Some(inv) = (self.transpose() * self.clone()).inverse() {
            return Some(inv * self.transpose()); 
        } else {
            return None; 
        }
    }

    /// 矩阵的右逆
    pub fn inverse_right(&self) -> Option<Self> {
        if self.row > self.col { return None; }
        if self.row == self.col { return self.inverse(); }
        if let Some(inv) = (self.clone() * self.transpose()).inverse() {
            return Some(self.transpose() * inv); 
        } else {
            return None; 
        }
    }
}

impl <T> Matrix<T> where
    T: Clone + One + Zero + Add<Output = T> + Sub<Output = T> + Neg<Output = T> + Mul<Output = T> + Div<Output = T> + PartialEq {
    /// 求矩阵的行列式（上三角法，适用于可阶梯化的类型）
    pub fn det(&self) -> T {
        assert!(self.row == self.col); 
        let (_self, independent_cols, neg) = self.step(false, false); 
        if independent_cols.len() < _self.row { return _self.kernel[0].get_zero(); }
        let mut det = _self.kernel[0].clone(); 
        for i in 1.._self.row { det = det.clone() * _self.kernel[i * _self.col + i].clone(); }
        if neg { det = det.neg(); }
        return det; 
    }
}

impl <T> Matrix<T> where
    T: Clone + Mul<Output = T> + Add<Output = T> + Sub<Output = T> + Zero {
    /// 按照定义求矩阵的行列式<br/>
    /// 速度慢，是对于不适用上三角法的类型的备用方法。
    /// 比如在求特征多项式时可用来求多项式矩阵的行列式。
    pub fn det_def(&self) -> T {
        assert!(self.row == self.col); 
        return self._det_def((0..self.row).collect::<Vec<usize>>(), (0..self.col).collect::<Vec<usize>>()); 
    }

    fn _det_def(&self, rows: Vec<usize>, cols: Vec<usize>) -> T {
        assert!(rows.len() == cols.len()); 
        if rows.len() == 1 { return self.kernel[rows[0] * self.col + cols[0]].clone(); }
        let mut det = self.kernel[0].get_zero(); 
        for i in 0..rows.len() {
            let (mut _rows, mut _cols) = (rows.clone(), cols.clone()); 
            _rows.remove(i); 
            _cols.remove(0); 
            if i % 2 == 1 {
                det = det.clone() - self.kernel[rows[i] * self.col + cols[0]].clone() * self._det_def(_rows, _cols); 
            } else {
                det = det.clone() + self.kernel[rows[i] * self.col + cols[0]].clone() * self._det_def(_rows, _cols); 
            }
        }
        return det; 
    }
}