#pragma once #include "vector.hpp" #include #include #include using namespace std; 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() { 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; } // Fixme: This is very dumb const vector get_eigenvalues() { 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; matrix m = this->hermitian_conjugate(); bool equal = true; for (long long i = 0; i < m.num_entries; i++) { for (long long j = 0; j < m.entry_dimension; j++) { if (m[i][j] != this->entries[i][j]) { equal = false; } } } return equal; } friend ostream &operator<<(ostream &os, const matrix &m) { char last = '\0'; for (long long i = 0; i < m.num_entries; i++) { for (long long j = 0; j < m.entry_dimension; j++) { string symbols[3]; symbols[0] = "|"; ostringstream oss; oss << " " << m.entries[i][j] << " "; 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 << 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]; } };