You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

281 lines
8.0 KiB

#pragma once
#include "vector.hpp"
#include <iostream>
#include <sstream>
#include <string>
class matrix {
private:
uint64_t num_entries;
uint64_t entry_dimension;
vector *entries;
bool err;
public:
matrix(const uint64_t num_entries, const uint64_t 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 (uint64_t 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 (uint64_t 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 uint64_t get_num_entries() const { return this->num_entries; }
const uint64_t 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 (uint64_t i = 0; i < diag.get_dimention(); i++) {
dsum = dsum * diag[i];
osum = osum * odiag[i];
}
return dsum - osum;
}
const vector get_diagonal() const {
uint64_t diag_len = this->entry_dimension;
if (this->num_entries < diag_len) {
diag_len = this->num_entries;
}
vector v(diag_len);
for (uint64_t i = 0; i < this->num_entries; i++) {
for (uint64_t 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().get_diagonal();
}
const matrix get_echelon() const {
matrix m = *this;
if (!m.is_invertible()) {
m.err = true;
} else {
for (uint64_t i = 0; i < m.get_num_entries(); i++) {
}
}
return m;
}
const vector get_entry(uint64_t index) const { return this->entries[index]; }
const vector get_row(uint64_t index) const {
vector v(this->entry_dimension);
for (uint64_t j = 0; j < this->entry_dimension; j++) {
v[j] = this->entries[j][index];
}
return v;
}
const matrix rotate() const {
matrix m = this->transpose();
matrix n = matrix(m.entry_dimension, m.num_entries);
for (uint64_t 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 (uint64_t i = 0; i < this->num_entries; i++) {
for (uint64_t 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 (uint64_t i = 0; i < this->num_entries; i++) {
for (uint64_t 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 (uint64_t i = 0; i < this->num_entries; i++) {
for (uint64_t 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 (uint64_t i = 0; i < this->num_entries; i++) {
for (uint64_t 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 (uint64_t i = 0; i < m.num_entries; i++) {
for (uint64_t 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 (uint64_t i = 0; i < m.num_entries; i++) {
for (uint64_t 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 (uint64_t i = 0; i < this->num_entries; i++) {
for (uint64_t j = 0; j < this->entry_dimension; j++) {
if (this->entries[i][j] != m[i][j]) {
equal = false;
}
}
}
return equal;
}
vector &operator[](const uint64_t index) { return this->entries[index]; }
const vector operator[](const uint64_t index) const {
return this->entries[index];
}
const matrix operator*(const cnumber z) const {
matrix n(this->num_entries, this->entry_dimension);
for (uint64_t 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 (uint64_t i = 0; i < this->get_num_entries(); i++) {
for (uint64_t 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 (uint64_t i = 0; i < this->entry_dimension; i++) {
for (uint64_t 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 (uint64_t 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 (uint64_t 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 (uint64_t i = 0; i < m.num_entries; i++)
this->entries[i] = m[i];
}
};