mmathlib/matrix.hpp

#pragma once
#include "vector.hpp"
#include <iostream>
#include <sstream>
#include <string>
class matrix {
private:
  long long num_entries;
  long long entry_dimension;
  vector *entries;
  bool err;

public:
  matrix(const long long num_entries, const long long entry_dimension) {
    this->num_entries = entry_dimension;
    this->entry_dimension = num_entries;
    this->err = false;
    this->entries = (vector *)malloc(sizeof(vector) * this->num_entries);
    for (long long i = 0; i < entry_dimension; i++) {
      this->entries[i] = vector(this->entry_dimension);
    }
  }
  matrix(const matrix &m) {
    this->num_entries = m.num_entries;
    this->entry_dimension = m.entry_dimension;
    this->err = m.err;
    this->entries = (vector *)malloc(sizeof(vector) * m.num_entries);
    for (long long i = 0; i < m.num_entries; i++) {
      this->entries[i] = m[i];
    }
  }
  ~matrix() {
    this->num_entries = 0;
    this->entry_dimension = 0;
    free(this->entries);
    this->entries = NULL;
    this->err = true;
  }
  const long long get_num_entries() const { return this->num_entries; }
  const long long get_entry_dimension() const { return this->entry_dimension; }
  const cnumber determinant() const {
    cnumber dsum(1, 0);
    cnumber osum(1, 0);
    vector diag = this->get_diagonal();
    vector odiag = this->get_off_diagonal();
    for (long long i = 0; i < diag.get_dimention(); i++) {
      dsum = dsum * diag[i];
      osum = osum * odiag[i];
    }
    return dsum - osum;
  }
  const vector get_diagonal() const {
    long long diag_len = this->entry_dimension;
    if (this->num_entries < diag_len) {
      diag_len = this->num_entries;
    }
    vector v(diag_len);
    for (long long i = 0; i < this->num_entries; i++) {
      for (long long j = 0; j < this->entry_dimension; j++) {
        if (i == j) {
          v[i] = this->entries[i][j];
        }
      }
    }
    return v;
  }
  const vector get_off_diagonal() const {
    return this->rotate_by_one_pi().get_diagonal();
  }
  const vector get_entry(long long index) const { return this->entries[index]; }

  const vector get_row(long long index) const {
    vector v(this->entry_dimension);
    for (long long j = 0; j < this->entry_dimension; j++) {
      v[j] = this->entries[j][index];
    }
    return v;
  }

  const matrix rotate_by_one_pi() const {
    matrix m = this->transpose();
    matrix n = matrix(m.entry_dimension, m.num_entries);
    for (long long i = 0; i < m.entry_dimension; i++) {
      int index = m.num_entries - i - 1;
      n[i] = m[index];
    }
    return n;
  }
  const matrix transpose() const {
    matrix n = matrix(this->entry_dimension, this->num_entries);
    for (long long i = 0; i < this->num_entries; i++) {
      for (long long j = 0; j < this->entry_dimension; j++) {
        n[j][i] = this->entries[i][j];
      }
    }

    return n;
  }
  const matrix conjugate() const {
    matrix n = matrix(this->num_entries, this->entry_dimension);
    for (long long i = 0; i < this->num_entries; i++) {
      for (long long j = 0; j < this->entry_dimension; j++) {
        n[i][j] = this->entries[i][j].conjugate();
      }
    }

    return n;
  }
  const matrix hermitian_conjugate() const {
    return this->transpose().conjugate();
  }
  const matrix I() const {
    matrix m = matrix(this->num_entries, this->entry_dimension);
    for (long long i = 0; i < this->num_entries; i++) {
      for (long long j = 0; j < this->entry_dimension; j++) {
        if (i == j) {
          m[i][j] = cnumber(1, 0);
        } else {
          m[i][j] = 0;
        }
      }
    }
    return m;
  }
  const bool is_diagonal() const {
    bool result = true;
    if (this->num_entries != this->entry_dimension)
      return !result;
    for (long long i = 0; i < this->num_entries; i++) {
      for (long long j = 0; j < this->num_entries; j++) {
        if (i != j && this->entries[i][j] != 0)
          return !result;
      }
    }
    return result;
  }
  const bool is_unitary() const{
	  matrix m = this->hermitian_conjugate();
	  return (*this * m == this->I());
  }
  // FIXME: This is very dumb
  const vector get_eigenvalues() const {
    if (this->is_diagonal())
      return this->get_diagonal();
    return vector(0);
  }
  const bool is_eigenvalue(cnumber z) const {
    return ((*this - (this->I() * z)).determinant() == 0);
  }
  const bool is_hermitian() const {
    if (this->entry_dimension != this->num_entries)
      return false;
    return *this == this->hermitian_conjugate();
  }
  friend std::ostream &operator<<(std::ostream &os, const matrix &m) {
    char last = '\0';
    int longest = 0;
    for (long long i = 0; i < m.num_entries; i++) {
      for (long long j = 0; j < m.entry_dimension; j++) {
        std::ostringstream oss;
        oss << m.entries[i][j];
        std::string s = oss.str();
        if (longest < s.length())
          longest = s.length();
      }
    }
    for (long long i = 0; i < m.num_entries; i++) {
      for (long long j = 0; j < m.entry_dimension; j++) {
        std::ostringstream iss;
        iss << m.entries[i][j];
        std::string s = iss.str();
        int padding = longest - s.length() +1;
        std::string symbols[3];
        symbols[0] = "|";
        std::ostringstream oss;
        oss << m.entries[i][j] ;
	oss << std::setw(padding) << "|";
        symbols[1] = oss.str();
        symbols[2] = "|";
        for (int i = 0; i < 3; i++) {
          int len = symbols[i].length() - 1;
          char cur = symbols[i][0];
          if (cur != last) {
            os << symbols[i];
          }
          last = symbols[i][len];
        }
      }
      if (i != m.num_entries - 1)
        os << std::endl << "|";
    }
    return os;
  }
  const bool operator==(const matrix &m) const {
    bool equal = true;
    for (long long i = 0; i < this->num_entries; i++) {
      for (long long j = 0; j < this->entry_dimension; j++) {
        if (this->entries[i][j] != m[i][j]) {
          equal = false;
        }
      }
    }
    return equal;
  }
  vector &operator[](const long long index) { return this->entries[index]; }
  const vector operator[](const long long index) const {
    return this->entries[index];
  }
  const matrix operator*(const cnumber z) const {
    matrix n(this->num_entries, this->entry_dimension);
    for (long long i = 0; i < this->num_entries; i++) {
      n[i] = this->entries[i] * z;
    }

    return n;
  }

  const matrix operator*(vector &v) const {
    matrix res(this->get_num_entries(), this->get_entry_dimension());
    for (long long i = 0; i < this->get_num_entries(); i++) {
      for (long long j = 0; j < this->get_entry_dimension(); j++) {
        res[i][j] = (this->get_entry(i)[j] * v.get_entry(i));
      }
    }

    return res;
  }
  const matrix operator*(const matrix m) const {
    matrix n(this->num_entries, m.entry_dimension);
    if (this->num_entries != m.entry_dimension &&
        m.num_entries != this->entry_dimension) {
      n.err = true;
      return n;
    }
    for (long long i = 0; i < this->entry_dimension; i++) {
      for (long long j = 0; j < this->entry_dimension; j++) {
        n[i][j] = this->get_entry(i) * m.get_row(j);
      }
    }

    return n;
  }

  const matrix operator+(const matrix &m) const {
    matrix n(this->num_entries, this->entry_dimension);
    for (long long i = 0; i < this->num_entries; i++) {
      n[i] = this->entries[i] + m[i];
    }

    return n;
  }
  const matrix operator-(const matrix &m) const {
    matrix n(this->num_entries, this->entry_dimension);
    for (long long i = 0; i < this->num_entries; i++) {
      n[i] = this->entries[i] - m[i];
    }

    return n;
  }
  void operator=(const matrix &m) {
    this->num_entries = m.num_entries;
    this->entry_dimension = m.entry_dimension;
    this->err = m.err;
    free(this->entries);
    this->entries = (vector *)malloc(sizeof(vector) * m.num_entries);
    for (long long i = 0; i < m.num_entries; i++)
      this->entries[i] = m[i];
  }
};