cmatrix.cc

/*
 *	Matrix maths
 *
 *	This program is distributed under the terms of the GPL v3.0 or later
 *	Download the GNU Public License (GPL) from <https://www.gnu.org>
 *
 *	Copyright(C) 2025, Free Software Foundation, Inc.
 *	Written by Nicholas Christopoulos <mailto:netnic@proton.me>
 */

#ifndef CMATRIX_H
#define CMATRIX_H	0x10000

#include <cstddef>
#include <cassert>
#include <cmath>

class CMatrix;
CMatrix operator+(const CMatrix& a, const CMatrix& b);
CMatrix operator-(const CMatrix& a, const CMatrix& b);
CMatrix operator*(const CMatrix& a, const CMatrix& b);

class CMatrix {
protected:
	size_t _rows = 0, _cols = 0;
	real_t *data = nullptr;
	
public:
	CMatrix() { _rows = _cols = 0; data = nullptr; }
	CMatrix(size_t r, size_t c, real_t val=0.0)	{ init(r, c, val); }
	CMatrix(const CMatrix &src) {
		init(src.rows(), src.cols());
		memcpy(data, src.data, src.size() * sizeof(real_t));
		}
	virtual ~CMatrix() { clear(); }

	inline void clear() {
		_rows = _cols = 0;
		if ( data )
			delete[] data;
		data = nullptr;
		}
	
	inline size_t rows() const { return _rows; }
	inline size_t cols() const { return _cols; }
	inline size_t size() const { return _rows * _cols; }

	inline real_t get(size_t r, size_t c) const
		{ assert(size() > r*_cols+c); return data[r*_cols+c]; }
	inline void set(size_t r, size_t c, real_t z)
		{ assert(size() > r*_cols+c); data[r*_cols+c] = z; }

	inline bool empty() const
		{ return (size() == 0); }

	// returns value of given location
	inline const real_t& operator()(size_t r, size_t c) const
		{ assert(size() > r*_cols+c); return data[r*_cols+c]; }
	inline real_t& operator()(size_t r, size_t c)
		{ assert(size() > r*_cols+c); return data[r*_cols+c]; }

	// assign
	const CMatrix& operator= (const CMatrix& src) {
		assert(src.size());
		clear();
		if ( src.size() ) {
			init(src.rows(), src.cols());
			memcpy(data, src.data, src.size() * sizeof(real_t));
			}
		return *this;
		}

	// fill
	inline const CMatrix& set(real_t k) {
		for ( size_t i = 0; i < size(); i ++ )
			data[i] = k;
		return *this;
		}

	// equal
	bool operator==(const CMatrix& src) const {
		if ( rows() != src.rows() ) return false;
		if ( cols() != src.cols() ) return false;
		for ( size_t i = 0; i < src.size(); i ++ )
			if ( data[i] != src.data[i] )
				return false;
		return true;
		}
	inline bool operator!=(const CMatrix& src) const
		{ return !(*this == src); }


	// unary -M
	CMatrix operator- () const {
		CMatrix x(rows(), cols());
		for ( size_t i = 0; i < size(); i ++ )
			x.data[i] = -data[i];
		return x;
		}

	//
	CMatrix operator+ (real_t k) const {
		CMatrix x(rows(), cols());
		for ( size_t i = 0; i < size(); i ++ )
			x.data[i] = data[i] + k;
		return x;
		}

	CMatrix operator- (real_t k) const {
		CMatrix x(rows(), cols());
		for ( size_t i = 0; i < size(); i ++ )
			x.data[i] = data[i] - k;
		return x;
		}

	CMatrix operator* (real_t k) const {
		CMatrix x(rows(), cols());
		for ( size_t i = 0; i < size(); i ++ )
			x.data[i] = data[i] * k;
		return x;
		}

	CMatrix operator/ (real_t k) const {
		if ( k == 0.0 ) throw std::runtime_error("Division by zero");
		CMatrix x(rows(), cols());
		for ( size_t i = 0; i < size(); i ++ )
			x.data[i] = data[i] / k;
		return x;
		}

	CMatrix power(int k) const {
		CMatrix x(rows(), cols(), 0.0);
		if ( size() == 0 || rows() != cols() )
			throw std::runtime_error("Matrix dimensions are not square (op ^)");
		for ( int n = 1; n < k; n ++ )
			x = x + (*this * *this);
		return x;
		}

	CMatrix identity() const;
	CMatrix transpose() const;
	CMatrix inverse() const;
	real_t determinant() const;
	
	// utils - print matrix (ASCII draw)
	void print() const;
	// utils - print matrix, CExp's syntax
	void lprint() const;

protected:
	void init(size_t r, size_t c, real_t val=0.0) {
		if ( r * c > 0 ) {
			_rows = r; _cols = c;
			data = new real_t[r*c];
			for ( size_t i = 0; i < size(); i ++ )
				data[i] = val;
			}
		}
	
	CMatrix get_minor(size_t row, size_t col) const;
	CMatrix complements() const;
	};

#endif

/* -------------------------------------------------------------------------------- */
#ifdef CMATRIX_IMPL
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
using namespace std;
	
// return sub matrix, delete row and col
CMatrix CMatrix::get_minor(size_t row, size_t col) const {
	if ( cols() > 1 && rows() > 1 ) {
		CMatrix minor(rows() - 1, cols() - 1);
		for ( size_t i = 0, x = 0; i < rows(); i ++ ) {
			if ( row != i ) {
				for ( size_t j = 0, y = 0; j < cols(); j ++ ) {
					if ( col != j ) {
						minor.set(x, y, get(i, j));
						y ++;
						}
					}
				x ++;
				}
	        }
		return minor;
		}
	throw runtime_error("cannot get minor Matrix");
	}

//
real_t CMatrix::determinant() const {
    real_t	det, t, res = 0.0;

	if ( cols() > 0 && cols() == rows() ) {
		if ( rows() == 1 )
			return get(0, 0);
		else {
			res = 0;
			for ( size_t i = 0; i < cols(); i ++ ) {
				CMatrix minor(get_minor(i, 0));
				det = minor.determinant();
				t = get(i, 0) * det;
				if ( (2 + i) % 2 == 1 )
					t = -t;
				res += t;
				}
			}
		return res;
		}
	throw runtime_error("not square matrix");
	}

//
CMatrix CMatrix::complements() const {
	real_t det_sub = 0.0;
		
	if ( cols() > 0 && cols() == rows() ) {
		CMatrix res(rows(), cols());
		for ( size_t i = 0; i < res.rows(); i ++ ) {
			for ( size_t j = 0; j < res.cols(); j ++ ) {
				CMatrix tmp(get_minor(i, j));
				det_sub = tmp.determinant();
				res.set(i, j, det_sub);
				if ( i + j % 2 == 1 )
					res.set(i, j, -res.get(i, j));
				}
			}
		return res;
		}
	throw runtime_error("not square matrix");
	}

//
CMatrix CMatrix::identity() const {
	if ( rows() != cols() )
		throw runtime_error("Matrix dimensions are not square (op I)");
	CMatrix x(rows(), cols());
	for ( size_t r = 0; r < x.rows(); r ++ ) {
		for ( size_t c = 0; c < x.cols(); c ++ )
			x.set(r, c, ( r == c ) ? 1.0 : 0.0);
		}
	return x;
	}

CMatrix CMatrix::transpose() const {
	CMatrix x(cols(), rows());
	for ( size_t r = 0; r < x.rows(); r ++ )
		for ( size_t c = 0; c < x.cols(); c ++ )
			x(r, c) = get(c, r);
	return x;
	}

//
CMatrix CMatrix::inverse() const {
	real_t det;

	CMatrix dop(complements());
	det = determinant();
	CMatrix dopT(dop.transpose());
	for ( size_t i = 0; i < dopT.rows(); i ++ )
		for ( size_t j = 0; j < dopT.cols(); j ++ )
			dopT.set(i, j, dopT.get(i, j) / det);
    return dopT;
	}

CMatrix operator+ (const CMatrix& a, const CMatrix& b) {
	if ( a.rows() != b.rows() || a.cols() != b.cols() )
		throw runtime_error("Matrices dimensions are not the same (op+)");
	CMatrix x(a.rows(), a.cols());
	for ( size_t r = 0; r < x.rows(); r ++ ) {
		for ( size_t c = 0; c < x.cols(); c ++ )
			x.set(r, c, a.get(r, c) + b.get(r, c));
		}
	return x;
	}

CMatrix operator- (const CMatrix& a, const CMatrix& b) {
	if ( a.rows() != b.rows() || a.cols() != b.cols() )
		throw runtime_error("Matrices dimensions are not the same (op+)");
	CMatrix x(a.rows(), a.cols());
	for ( size_t r = 0; r < x.rows(); r ++ ) {
		for ( size_t c = 0; c < x.cols(); c ++ )
			x.set(r, c, a.get(r, c) - b.get(r, c));
		}
	return x;
	}

// multiply
CMatrix operator*(const CMatrix& a, const CMatrix& b) {
	if ( a.rows() != b.cols() || a.cols() != b.rows() )
		throw runtime_error("Matrices dimensions are not the reversed (op*)");
	CMatrix x(a.rows(), b.cols());
	for ( size_t i = 0; i < a.rows(); i ++ ) {
		for ( size_t j = 0; j < b.cols(); j ++ ) {
			x.set(i, j, 0);
			for ( size_t k = 0; k < a.cols(); k ++ )
				x.set(i, j, x.get(i, j) + a.get(i, k) * b.get(k, j) );
			}
		}
	return x;
	}

// print matrix ; multiline
void CMatrix::print() const {
	size_t r, c;
	printf("\t");
	for ( c = 0; c < cols(); c ++ )
		printf("%zd\t", c+1);
	printf("\n");
	
	for ( r = 0; r < rows(); r ++ ) {
		printf("%zd:\t", r+1);
		for ( c = 0; c < cols(); c ++ )
			printf("%g\t", get(r, c));
		printf("\n");
		}
	}

// print matrix ; single line (`;` = row separator)
void CMatrix::lprint() const {
	size_t r, c;
	printf("[ ");
	for ( r = 0; r < rows(); r ++ ) {
		for ( c = 0; c < cols(); c ++ ) {
			printf("%g", get(r, c));
			if ( c < cols() - 1 )
				printf(", ");
			}
		if ( r < rows() - 1 )
			printf(" ; ");
		}
	printf(" ]\n");
	}

#endif // implementation