|
|
@ -14,6 +14,7 @@ private:
|
|
|
|
bool augmented = false;
|
|
|
|
bool augmented = false;
|
|
|
|
|
|
|
|
|
|
|
|
public:
|
|
|
|
public:
|
|
|
|
|
|
|
|
// Constructors/destructor
|
|
|
|
matrix(const uint64_t num_entries, const uint64_t entry_dimension) {
|
|
|
|
matrix(const uint64_t num_entries, const uint64_t entry_dimension) {
|
|
|
|
this->num_entries = entry_dimension;
|
|
|
|
this->num_entries = entry_dimension;
|
|
|
|
this->entry_dimension = num_entries;
|
|
|
|
this->entry_dimension = num_entries;
|
|
|
@ -39,6 +40,13 @@ public:
|
|
|
|
this->entries = NULL;
|
|
|
|
this->entries = NULL;
|
|
|
|
this->err = true;
|
|
|
|
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 {
|
|
|
|
const matrix augment(vector v) const {
|
|
|
|
matrix m = *this;
|
|
|
|
matrix m = *this;
|
|
|
|
if (v.get_dimention() != this->get_num_entries()) {
|
|
|
|
if (v.get_dimention() != this->get_num_entries()) {
|
|
|
@ -57,8 +65,16 @@ public:
|
|
|
|
m.augmented = true;
|
|
|
|
m.augmented = true;
|
|
|
|
return m;
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
const uint64_t get_num_entries() const { return this->num_entries; }
|
|
|
|
const matrix conjugate() const {
|
|
|
|
const uint64_t get_entry_dimension() const { return this->entry_dimension; }
|
|
|
|
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 {
|
|
|
|
const cnumber determinant() const {
|
|
|
|
cnumber dsum(1, 0);
|
|
|
|
cnumber dsum(1, 0);
|
|
|
|
cnumber osum(1, 0);
|
|
|
|
cnumber osum(1, 0);
|
|
|
@ -70,6 +86,22 @@ public:
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return dsum - osum;
|
|
|
|
return dsum - osum;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// 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 {
|
|
|
|
const vector get_diagonal() const {
|
|
|
|
uint64_t diag_len = this->entry_dimension;
|
|
|
|
uint64_t diag_len = this->entry_dimension;
|
|
|
|
if (this->num_entries < diag_len) {
|
|
|
|
if (this->num_entries < diag_len) {
|
|
|
@ -85,10 +117,6 @@ public:
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return v;
|
|
|
|
return v;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
const vector get_off_diagonal() const {
|
|
|
|
|
|
|
|
return this->rotate().get_diagonal();
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
const matrix get_echelon() const {
|
|
|
|
const matrix get_echelon() const {
|
|
|
|
matrix m = *this;
|
|
|
|
matrix m = *this;
|
|
|
|
if (!m.is_invertible()) {
|
|
|
|
if (!m.is_invertible()) {
|
|
|
@ -97,72 +125,17 @@ public:
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return m;
|
|
|
|
return m;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// FIXME: This is very dumb
|
|
|
|
// switch row i with row j
|
|
|
|
const vector get_eigenvalues() const {
|
|
|
|
const matrix exchange_row(uint64_t i, uint64_t j) {
|
|
|
|
if (this->is_diagonal())
|
|
|
|
matrix m = *this;
|
|
|
|
return this->get_diagonal();
|
|
|
|
vector v = m.get_entry(i);
|
|
|
|
return vector(0);
|
|
|
|
m[i] = m[j];
|
|
|
|
|
|
|
|
m[j] = v;
|
|
|
|
|
|
|
|
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;
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
// add row i to row j
|
|
|
|
|
|
|
|
const matrix add_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;
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
// Multiply row by z
|
|
|
|
|
|
|
|
const matrix multiply_row(uint64_t i, cnumber z) {
|
|
|
|
|
|
|
|
matrix m = *this;
|
|
|
|
|
|
|
|
m[i] = m[i] * z;
|
|
|
|
|
|
|
|
return m;
|
|
|
|
|
|
|
|
}
|
|
|
|
}
|
|
|
|
const vector get_entry(uint64_t index) const { return this->entries[index]; }
|
|
|
|
const vector get_entry(uint64_t index) const { return this->entries[index]; }
|
|
|
|
|
|
|
|
const uint64_t get_entry_dimension() const { return this->entry_dimension; }
|
|
|
|
const vector get_column(uint64_t index) const {
|
|
|
|
const uint64_t get_num_entries() const { return this->num_entries; }
|
|
|
|
vector v(this->entry_dimension);
|
|
|
|
const vector get_off_diagonal() const {
|
|
|
|
for (uint64_t j = 0; j < this->entry_dimension; j++) {
|
|
|
|
return this->rotate().get_diagonal();
|
|
|
|
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 {
|
|
|
|
const matrix hermitian_conjugate() const {
|
|
|
|
return this->transpose().conjugate();
|
|
|
|
return this->transpose().conjugate();
|
|
|
@ -192,12 +165,6 @@ public:
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return result;
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
// 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 {
|
|
|
|
const bool is_eigenvalue(cnumber z) const {
|
|
|
|
return ((*this - (this->I() * z)).determinant() == 0);
|
|
|
|
return ((*this - (this->I() * z)).determinant() == 0);
|
|
|
|
}
|
|
|
|
}
|
|
|
@ -211,6 +178,40 @@ public:
|
|
|
|
matrix m = this->hermitian_conjugate();
|
|
|
|
matrix m = this->hermitian_conjugate();
|
|
|
|
return (*this * m == this->I());
|
|
|
|
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 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) {
|
|
|
|
friend std::ostream &operator<<(std::ostream &os, const matrix &m) {
|
|
|
|
char last = '\0';
|
|
|
|
char last = '\0';
|
|
|
|
int longest = 0;
|
|
|
|
int longest = 0;
|
|
|
@ -240,7 +241,8 @@ public:
|
|
|
|
for (int k = 0; k < 3; k++) {
|
|
|
|
for (int k = 0; k < 3; k++) {
|
|
|
|
int len = symbols[k].length() - 1;
|
|
|
|
int len = symbols[k].length() - 1;
|
|
|
|
char cur = symbols[k][0];
|
|
|
|
char cur = symbols[k][0];
|
|
|
|
if (cur != last || (m.augmented && j == m[i].get_dimention() - 1 && print)) {
|
|
|
|
if (cur != last ||
|
|
|
|
|
|
|
|
(m.augmented && j == m[i].get_dimention() - 1 && print)) {
|
|
|
|
print = false;
|
|
|
|
print = false;
|
|
|
|
os << symbols[k];
|
|
|
|
os << symbols[k];
|
|
|
|
}
|
|
|
|
}
|
|
|
|