|
|
|
#pragma once
|
|
|
|
#include "cnumber.hpp"
|
|
|
|
#include "vector.hpp"
|
|
|
|
#include <cstdint>
|
|
|
|
#include <iostream>
|
|
|
|
#include <ostream>
|
|
|
|
#include <sstream>
|
|
|
|
#include <string>
|
|
|
|
class matrix {
|
|
|
|
private:
|
|
|
|
uint64_t num_entries;
|
|
|
|
uint64_t entry_dimension;
|
|
|
|
vector *entries;
|
|
|
|
bool err;
|
|
|
|
bool augmented = false;
|
|
|
|
// recursive determinant
|
|
|
|
const cnumber rdeterminant(matrix m) const {
|
|
|
|
if (m.get_num_entries() == 2) {
|
|
|
|
cnumber dsum(1, 0);
|
|
|
|
cnumber osum(1, 0);
|
|
|
|
vector diag = m.get_diagonal();
|
|
|
|
vector odiag = m.get_off_diagonal();
|
|
|
|
for (uint64_t i = 0; i < diag.get_dimention(); i++) {
|
|
|
|
dsum = dsum * diag[i];
|
|
|
|
osum = osum * odiag[i];
|
|
|
|
}
|
|
|
|
return dsum - osum;
|
|
|
|
} else {
|
|
|
|
cnumber sum = 0;
|
|
|
|
cnumber sign(1);
|
|
|
|
for (uint64_t i = 0; i < m.entry_dimension; i++) {
|
|
|
|
matrix n = m.remove_first_row_and_column(i);
|
|
|
|
if(i % 2 == 0)
|
|
|
|
sum = sum + (m[0][i] * rdeterminant(n));
|
|
|
|
else
|
|
|
|
sum = sum - (m[0][i] * rdeterminant(n));
|
|
|
|
}
|
|
|
|
return sum;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
public:
|
|
|
|
// Constructors/destructor
|
|
|
|
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;
|
|
|
|
}
|
|
|
|
// member functions
|
|
|
|
const matrix add_row(uint64_t i, uint64_t j, cnumber multiplier = 1) {
|
|
|
|
// add row i to row j
|
|
|
|
matrix m = *this;
|
|
|
|
m[j] = m[j] + (m.multiply_row(i, multiplier)[i]);
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
const matrix augment(vector v) const {
|
|
|
|
matrix m = *this;
|
|
|
|
if (v.get_dimention() != this->get_num_entries()) {
|
|
|
|
m.err = true;
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
for (uint64_t i = 0; i < m.get_num_entries(); i++) {
|
|
|
|
vector old = m[i];
|
|
|
|
vector n(old.get_dimention() + 1);
|
|
|
|
for (uint64_t j = 0; j < old.get_dimention(); j++)
|
|
|
|
n[j] = old[j];
|
|
|
|
n[old.get_dimention()] = v[i];
|
|
|
|
m[i] = n;
|
|
|
|
}
|
|
|
|
m.entry_dimension = m.entry_dimension + 1;
|
|
|
|
m.augmented = true;
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
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 cnumber determinant() const {
|
|
|
|
if (this->num_entries != this->entry_dimension)
|
|
|
|
return 0;
|
|
|
|
return rdeterminant(*this);
|
|
|
|
}
|
|
|
|
// switch row i with row j
|
|
|
|
const matrix exchange_row(uint64_t i, uint64_t j) {
|
|
|
|
matrix m = *this;
|
|
|
|
vector v = m.get_entry(i);
|
|
|
|
m[i] = m[j];
|
|
|
|
m[j] = v;
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
const vector get_column(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 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 matrix get_echelon() const {
|
|
|
|
matrix m = *this;
|
|
|
|
if (!m.is_invertible()) {
|
|
|
|
m.err = true;
|
|
|
|
} else {
|
|
|
|
uint64_t num_entries = m.get_num_entries();
|
|
|
|
for (uint64_t i = 0; i < num_entries - 1; i++) {
|
|
|
|
for (uint64_t j = i; j < num_entries; j++) {
|
|
|
|
cnumber ratio = m[j][i] / m[i][i];
|
|
|
|
for (uint64_t k = i; k < num_entries; k++) {
|
|
|
|
m[j][k] = m[j][k] - (ratio * m[i][k]);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
// FIXME: This is very dumb
|
|
|
|
const vector get_eigenvalues() const {
|
|
|
|
if (this->is_diagonal())
|
|
|
|
return this->get_diagonal();
|
|
|
|
return vector(0);
|
|
|
|
}
|
|
|
|
const vector get_entry(uint64_t index) const { return this->entries[index]; }
|
|
|
|
const uint64_t get_entry_dimension() const { return this->entry_dimension; }
|
|
|
|
const uint64_t get_num_entries() const { return this->num_entries; }
|
|
|
|
const vector get_off_diagonal() const {
|
|
|
|
return this->rotate().get_diagonal();
|
|
|
|
}
|
|
|
|
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_broken() const { return this->err; }
|
|
|
|
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_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();
|
|
|
|
}
|
|
|
|
const bool is_invertible() const { return this->determinant() != 0; }
|
|
|
|
const bool is_unitary() const {
|
|
|
|
matrix m = this->hermitian_conjugate();
|
|
|
|
return (*this * m == this->I());
|
|
|
|
}
|
|
|
|
// Multiply row by z
|
|
|
|
const matrix multiply_row(uint64_t i, cnumber z) {
|
|
|
|
matrix m = *this;
|
|
|
|
m[i] = m[i] * z;
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
|
|
|
|
// subtract row i from row j
|
|
|
|
const matrix subtract_row(uint64_t i, uint64_t j, cnumber multiplier = 1) {
|
|
|
|
matrix m = *this;
|
|
|
|
m[j] = m[j] - (m.multiply_row(i, multiplier)[i]);
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
|
|
|
|
const matrix remove_first_row_and_column(uint64_t column) const {
|
|
|
|
matrix m(this->num_entries - 1, this->entry_dimension - 1);
|
|
|
|
uint64_t add = 0;
|
|
|
|
for (uint64_t i = 1; i < this->num_entries; i++) {
|
|
|
|
for (uint64_t j = 0; j < this->entry_dimension; j++) {
|
|
|
|
if (j != column) {
|
|
|
|
m[i - 1][j - add] = this->get_entry(i)[j];
|
|
|
|
} else {
|
|
|
|
add = 1;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
add = 0;
|
|
|
|
}
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
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;
|
|
|
|
}
|
|
|
|
// Operators
|
|
|
|
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] = "|";
|
|
|
|
bool print = true;
|
|
|
|
for (int k = 0; k < 3; k++) {
|
|
|
|
int len = symbols[k].length() - 1;
|
|
|
|
char cur = symbols[k][0];
|
|
|
|
if (cur != last ||
|
|
|
|
(m.augmented && j == m[i].get_dimention() - 1 && print)) {
|
|
|
|
print = false;
|
|
|
|
os << symbols[k];
|
|
|
|
}
|
|
|
|
last = symbols[k][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_column(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];
|
|
|
|
}
|
|
|
|
};
|