main.c

// Artemii Miasoedov
// a.miasoedov@innopolis.university
#include <iostream>
#include <memory>
#include <functional>
#include <type_traits>
#include <sstream>
#include <iomanip>
#include <limits>
#include <cmath>
#include <fstream>

namespace la {
	/**
	 * Exception, thrown, when invalid matrix size is encountered
	 */
	class dimensions_error : public std::runtime_error {
	public:
		explicit dimensions_error()
			: runtime_error("Error: the dimensional problem occurred") {
		}

		explicit dimensions_error(const std::string &message)
			: runtime_error(message) {
		}
	};

	/**
	 * Steps of elimination process
	 */
	enum EliminationStep {
		ELIMINATION,
		PERMUTATION,
		BACK_ELIMINATION,
		DIAG_NORM
	};

	// epsilon values for float and double types, used to compare with 0 and other values
	static thread_local float f_epsilon = std::numeric_limits<float>::epsilon();
	static thread_local double d_epsilon = std::numeric_limits<double>::epsilon();

	template<typename T>
	requires std::is_arithmetic_v<T>
	inline T epsilon();

	// returns epsilon for float type
	template<>
	inline float epsilon<float>() {
		return f_epsilon;
	}

	// returns epsilon for double type
	template<>
	inline double epsilon<double>() {
		return d_epsilon;
	}

	template<typename T>
	requires std::is_arithmetic_v<T>
	inline bool is_zero(T t);

	// compares int value with 0
	template<>
	inline bool is_zero<int>(int t) {
		return t == 0;
	}

	// compares float value with 0
	template<>
	inline bool is_zero<float>(float t) {
		return std::fabs(t) < epsilon<float>();
	}

	// compares double value with 0
	template<>
	inline bool is_zero<double>(double t) {
		return std::fabs(t) < epsilon<double>();
	}

	template<typename T>
	requires std::is_floating_point_v<T>
	bool is_less_than_10_nn(T value, int n) {
		auto t = std::abs(value) * std::pow(10, n);
		if (t > 1) return false;
		return t * 10 < 5;
		// return std::abs(value) < std::pow(10, (T) -n);
	}

	template<typename T> requires std::is_arithmetic_v<T>
	class Matrix;

	template<typename T> requires std::is_arithmetic_v<T>
	class ColumnVector;

	/**
	 * Generic Matrix class, supporting different operations such as: multiplication, summation, and others
	 * @tparam T type of Matrix values
	 */
	template<typename T> requires std::is_arithmetic_v<T>
	class Matrix {
	protected:

		using dim_t = int;

		dim_t rows_ = 0, cols_ = 0;
		std::unique_ptr<std::unique_ptr<T[]>[]> row_vectors_;

	public:

		Matrix()
			: row_vectors_(nullptr) {
		}

		/**
		 * Creates new matrix with specified width and height filled with 0s
		 * @param height number of rows in matrix
		 * @param width number of cols in matrix
		 */
		Matrix(dim_t height, dim_t width)
			: cols_(width), rows_(height), row_vectors_(std::make_unique<std::unique_ptr<T[]>[]>(rows_)) {
			for (dim_t i = 0; i < rows_; ++i) {
				std::unique_ptr<T[]> &row = (row_vectors_[i] = std::make_unique<T[]>(cols_));
				std::fill(&row[0], &row[cols_], 0);
			}
		}

		/**
		 * Creates new square matrix filled with 0s
		 * @param size number of rows and cols in the square matrix
		 */
		explicit Matrix(dim_t size)
			: Matrix(size, size) {
		}

		/**
		 * Creates matrix from initializer list
		 * @param m initializer list
		 */
		Matrix(const std::initializer_list<std::initializer_list<T>> &m)
			: Matrix(m.size(), m.begin()->size()) {
			dim_t y = 0, x = 0;

			for (const auto &row: m) {
				if (row.size() != cols_) {
					throw dimensions_error("All rows must have the same number of columns");
				}

				for (const auto &item: row) {
					at(y, x) = item;
					x++;
				}
				x = 0;
				y++;
			}
		}

		/**
		 * Copy constructor
		 * @param m matrix to copy
		 */
		Matrix(const Matrix<T> &m)
			: Matrix(m.rows_, m.cols_) {
			for (dim_t i = 0; i < rows_; ++i) {
				std::copy(&m.row_vectors_[i][0], &m.row_vectors_[i][cols_], &row_vectors_[i][0]);
			}
		}

		/**
		 * Move constructor
		 * @param m matrix to move
		 */
		Matrix(Matrix<T> &&m) noexcept
			: rows_(m.rows_), cols_(m.cols_) {
			row_vectors_ = std::move(m.row_vectors_);
		}

		/**
		 * Creates identity matrix with specified size
		 * @param size number of rows and cols in the matrix
		 * @return identity matrix with specified size
		 */
		static Matrix<T> identity(dim_t size) {
			Matrix<T> result(size);

			for (int i = 0; i < size; ++i) {
				result[i][i] = 1;
			}

			return result;
		}

		/**
		 * Creates permutation matrix for specified matrix
		 * @param m matrix, which rows/cols should be swapped
		 * @param a first row/column to swap
		 * @param b second row/column to swap
		 * @return new permutation matrix
		 */
		static Matrix<T> permutation(const Matrix<T> &m, dim_t a, dim_t b) {
			Matrix<T> result = Matrix<T>::identity(m.rows_);

			result[a][a] = result[b][b] = 0;
			result[a][b] = result[b][a] = 1;

			return result;
		}

		/**
		 * Finds pivot in the specified column
		 * @param col where to find a pivot
		 * @return row where pivot is located and pivot itself
		 */
		std::tuple<dim_t, T> pivot(dim_t col) {
			dim_t pivot_row = col;
			T pivot = std::abs(at(col, col));

			for (dim_t row = col + 1; row < rows_; ++row) {
				// todo:
				T e = std::abs(at(row, col));
				if (pivot < e) {
					pivot = e;
					pivot_row = row;
				}
			}

			return std::make_tuple(pivot_row, at(pivot_row, col));
		}

		using elimination_callback = std::function<void(EliminationStep, dim_t, dim_t, T)>;

		/**
		 * Subtracts one row from another
		 * @param row_a subtrahend row
		 * @param row_b minuend row
		 * @param scalar factor of a subtrahend row
		 */
		void eliminate(dim_t row_a, dim_t row_b, T scalar) {
			for (int i = 0; i < cols_; ++i) {
				at(row_a, i) -= scalar * at(row_b, i);
			}
		}

		/**
		 * Swaps two rows in the matrix
		 * @param row_a first row index to swap
		 * @param row_b second row index to swap
		 */
		void swap(dim_t row_a, dim_t row_b) {
			std::swap(row_vectors_[row_a], row_vectors_[row_b]);
		}

		/**
		 * Divides row by the specified factor
		 * @param row row to divide
		 * @param scalar quotient
		 */
		void div_row(dim_t row, T scalar) {
			auto &r = row_vectors_[row];
			for (int i = 0; i < cols_; ++i) {
				r[i] /= scalar;
			}
		}

		// todo: check cols >= rows
		/**
		 * Performs row echelon form elimination
		 * @param callback callback function for each step of REF elimination
		 */
		void ref(const elimination_callback &callback = elimination_callback()) {
			for (dim_t col = 0; col < rows_ - 1; ++col) {
				auto [p_row, p] = pivot(col);
				if (la::is_zero(p)) continue;

				if (p_row != col) {
					std::swap(row_vectors_[col], row_vectors_[p_row]);
					callback(PERMUTATION, col, p_row, 0);
					p_row = col;
				}

				for (dim_t e_row = col + 1; e_row < rows_; ++e_row) {
					T scalar = at(e_row, col) / p;
					at(e_row, col) = 0;

					for (dim_t e_col = col + 1; e_col < cols_; ++e_col) {
						at(e_row, e_col) -= (scalar * at(p_row, e_col));
					}

					callback(ELIMINATION, e_row, p_row, scalar);
				}
			}
		}

		/**
		 * Performs reduced row echelon form elimination
		 * @param callback callback function for each step of RREF elimination
		 */
		void rref(const elimination_callback &callback = elimination_callback()) {
			auto max_pivots = std::min(rows_, cols_);

			for (dim_t col = 0; col < max_pivots; ++col) {
				auto [p_row, p] = pivot(col);
				if (is_zero(p)) continue;

				if (p_row != col) {
					std::swap(row_vectors_[col], row_vectors_[p_row]);
					callback(PERMUTATION, col, p_row, 0);
					p_row = col;
				}

				for (dim_t e_row = col + 1; e_row < rows_; ++e_row) {
					T scalar = at(e_row, col) / p;
					if (is_zero(scalar)) continue;

					at(e_row, col) = 0;

					for (dim_t e_col = col; e_col < cols_; ++e_col) {
						at(e_row, e_col) -= (scalar * at(p_row, e_col));
					}

					callback(ELIMINATION, e_row, p_row, scalar);
				}
			}

			callback(BACK_ELIMINATION, 0, 0, 0);

			for (dim_t col = max_pivots - 1; col >= 0; --col) {
				T p = at(col, col);
				if (is_zero(p)) continue;

				for (dim_t e_row = col - 1; e_row >= 0; --e_row) {
					T scalar = at(e_row, col) / p;
					if (is_zero(scalar)) continue;

					at(e_row, col) = 0;

					for (dim_t e_col = col; e_col < cols_; ++e_col) {
						at(e_row, e_col) -= (scalar * at(col, e_col));
					}

					callback(ELIMINATION, e_row, col, scalar);
				}
			}

			for (int i = 0; i < diag(); ++i) {
				T scalar = at(i, i);
				if (is_zero(scalar) || is_zero(1 - scalar)) continue;

				div_row(i, scalar);
				callback(DIAG_NORM, i, 0, scalar);
			}
		}

		/**
		 * Assignment operator
		 * @param m value to be copied
		 * @return copied value
		 */
		Matrix<T> &operator =(const Matrix<T> &m) {
			cols_ = m.cols_;
			rows_ = m.rows_;
			row_vectors_ = std::make_unique<std::unique_ptr<T[]>[]>(rows_);

			for (dim_t i = 0; i < rows_; ++i) {
				std::unique_ptr<T[]> &row = (row_vectors_[i] = std::make_unique<T[]>(cols_));
				std::copy(&row[0], &row[cols_], &row_vectors_[0]);
			}

			return *this;
		}

		/**
		 * Assignment operator
		 * @param m value to moved
		 * @return moved value
		 */
		Matrix<T> &operator =(Matrix<T> &&m) noexcept {
			cols_ = m.cols_;
			rows_ = m.rows_;
			row_vectors_ = std::move(m.row_vectors_);

			return *this;
		}

		/**
		 * Returns element at specified position
		 * @param row index of value row
		 * @param col index of value column
		 * @return value at specified position
		 */
		[[nodiscard]] virtual const T &at(dim_t row, dim_t col) const {
			return row_vectors_[row][col];
		}

		/**
		 * Returns element at specified position
		 * @param row index of value row
		 * @param col index of value column
		 * @return value at specified position
		 */
		virtual T &at(dim_t row, dim_t col) {
			return row_vectors_[row][col];
		}

		/**
		 * Returns row with specified index
		 * @param row index of row to return
		 * @return row with specified index
		 */
		const T *operator [](dim_t row) const {
			return &row_vectors_[row][0];
		}

		/**
		 * Returns row with specified index
		 * @param row index of row to return
		 * @return row with specified index
		 */
		virtual T *operator [](dim_t row) {
			return &row_vectors_[row][0];
		}

		/**
		 * Sums 2 matrices
		 * @param b summand matrix
		 * @return sum of 2 matrices
		 */
		virtual Matrix<T> operator +(const Matrix<T> &b) const {
			if (b.rows_ != rows_ || b.cols_ != cols_) {
				throw dimensions_error();
			}

			Matrix<T> result(rows_, cols_);

			for (dim_t i = 0; i < rows_; ++i) {
				for (dim_t j = 0; j < cols_; ++j) {
					result[i][j] = at(i, j) + b[i][j];
				}
			}

			return result;
		}

		/**
		 * Subtracts 2 matrices
		 * @param b subtrahend matrix
		 * @return difference of 2 matrices
		 */
		virtual Matrix<T> operator -(const Matrix<T> &b) const {
			if (b.rows_ != rows_ || b.cols_ != cols_) {
				throw dimensions_error();
			}

			Matrix<T> result(rows_, cols_);

			for (dim_t i = 0; i < rows_; ++i) {
				for (dim_t j = 0; j < cols_; ++j) {
					result[i][j] = at(i, j) - b[i][j];
				}
			}

			return result;
		}

		/**
		 * Multiplies 2 matrices
		 * @param b multiplier matrix
		 * @return the result of matrix multiplication
		 */
		virtual Matrix<T> operator *(const Matrix<T> &b) const {
			if (b.rows_ != cols_) {
				throw dimensions_error();
			}

			Matrix<T> result(rows_, b.cols_);

			for (dim_t i = 0; i < rows_; ++i) {
				for (dim_t j = 0; j < b.cols_; ++j) {
					T e_ij = 0;
					for (dim_t k = 0; k < cols_; ++k) {
						e_ij += at(i, k) * b[k][j];
					}
					result[i][j] = e_ij;
				}
			}

			return result;
		}

		/**
		 * Multiplies matrix by scalar
		 * @param scalar scalar to be multiplied by
		 * @return the result of multiplication of matrix by scalar
		 */
		virtual Matrix<T> operator *(const T scalar) const {
			Matrix<T> result(rows_, cols_);

			for (dim_t i = 0; i < rows_; ++i) {
				for (dim_t j = 0; j < cols_; ++j) {
					result[i][j] = scalar * at(i, j);
				}
			}

			return result;
		}

		virtual ColumnVector<T> operator *(const ColumnVector<T> &vector) const;

		/**
		 * Transposes the matrix
		 * @return transposition result
		 */
		[[nodiscard]] virtual Matrix<T> transpose() const {
			T transposed_width = rows_;
			T transposed_height = cols_;
			Matrix<T> result(transposed_height, transposed_width);

			for (dim_t i = 0; i < rows_; ++i) {
				for (dim_t j = 0; j < cols_; ++j) {
					result[j][i] = at(i, j);
				}
			}

			return result;
		}

		[[nodiscard]] virtual Matrix<T> inverse() const {
			if (cols_ != rows_) {
				throw dimensions_error();
			}

			Matrix<T> tmp = *this;
			Matrix<T> result = Matrix<T>::identity(rows_);

			tmp.rref(
				[&result](la::EliminationStep step, int row_a, int row_b, double value) {
					switch (step) {
					case la::ELIMINATION:
						result.eliminate(row_a, row_b, value);
						break;
					case la::PERMUTATION:
						result.swap(row_a, row_b);
						break;
					case la::DIAG_NORM:
						result.div_row(row_a, value);
						break;
					default:
						break;
					}
				}
			);

			return result;
		}

		[[nodiscard]] dim_t rows() const {
			return rows_;
		}

		[[nodiscard]] dim_t cols() const {
			return cols_;
		}

		[[nodiscard]] dim_t diag() const {
			return std::min(cols_, rows_);
		}

		[[nodiscard]] dim_t order() const {
			return rows_ * cols_;
		}

		/**
		 * Reads matrix from stream
		 * @param is stream to read matrix from
		 * @param m matrix, where values will be read
		 * @return passed in stream
		 */
		friend std::istream &operator >>(std::istream &is, Matrix<T> &m) {
			for (dim_t i = 0; i < m.rows_; ++i) {
				for (dim_t j = 0; j < m.cols_; ++j) {
					is >> m[i][j];
				}
			}

			return is;
		}

		/**
		 * Prints the row from matrix
		 * @param os stream to print to
		 * @param row index of row to be printed
		 * @return passed in stream
		 */
		std::ostream &print_row(std::ostream &os, dim_t row) const {
			auto precision = os.precision();

			for (dim_t j = 0; j < cols_; ++j) {
				const auto value = at(row, j);
				os << (la::is_less_than_10_nn(value, precision) ? 0 : value);
				if (j < cols_ - 1) {
					os << ' ';
				}
			}

			return os;
		}

		/**
		 * Prints the matrix to the stream
		 * @param os stream to print matrix to
		 * @param m matrix to be printed
		 * @return passed in stream
		 */
		friend std::ostream &operator <<(std::ostream &os, const Matrix<T> &m) {
			for (dim_t i = 0; i < m.rows_; ++i) {
				m.print_row(os, i);
				os << std::endl;
			}

			return os;
		}
	};

	template<>
	std::ostream &Matrix<int>::print_row(std::ostream &os, la::Matrix<int>::dim_t row) const {
		for (dim_t j = 0; j < cols_; ++j) {
			os << at(row, j);
			if (j < cols_ - 1) {
				os << ' ';
			}
		}
		return os;
	}

	/**
	 * Column vector
	 * @tparam T type of elements in the ColumnVector
	 */
	template<typename T> requires std::is_arithmetic_v<T>
	class ColumnVector {
		using dim_t = unsigned int;

		dim_t rows_ = 0;
		std::unique_ptr<T[]> column_;

	public:

		/**
		 * Creates new ColumnVector with specified number of rows
		 * @param rows number of rows in ColumnVector
		 */
		explicit ColumnVector(dim_t rows)
			: rows_(rows), column_(std::make_unique<T[]>(rows)) {
		}

		ColumnVector(const ColumnVector<T> &m)
			: ColumnVector(m.rows_) {
			for (dim_t i = 0; i < rows_; ++i) {
				this->column_[i] = m.column_[i];
			}
		}

		ColumnVector<T> &operator =(const ColumnVector<T> &m) {
			rows_ = m.rows_;
			column_ = std::make_unique<T[]>(rows_);

			for (dim_t i = 0; i < rows_; ++i) {
				this->column_[i] = m.column_[i];
			}

			return *this;
		}

		[[nodiscard]] dim_t rows() const {
			return rows_;
		}

		T &operator [](dim_t row) {
			return column_[row];
		}

		const T &operator [](dim_t row) const {
			return column_[row];
		}

		[[nodiscard]] T norm() const {
			T r = 0;
			for (int i = 0; i < rows_; ++i) {
				r += (column_[i] * column_[i]);
			}
			return std::sqrt(r);
		}

		/**
		 * Sums 2 ColumnVectors
		 * @param a first summand
		 * @param b second summand
		 * @return sum of 2 vectors
		 */
		friend ColumnVector<T> operator +(const ColumnVector<T> &a, const ColumnVector<T> &b) {
			if (a.rows_ != b.rows_) {
				throw dimensions_error("Vectors should have the same number of rows");
			}

			ColumnVector<T> r(a.rows_);
			for (int i = 0; i < a.rows_; ++i) {
				r[i] = a[i] + b[i];
			}

			return r;
		}

		friend ColumnVector<T> operator -(const ColumnVector<T> &a, const ColumnVector<T> &b) {
			if (a.rows_ != b.rows_) {
				throw dimensions_error("Vectors should have the same number of rows");
			}

			ColumnVector<T> r(a.rows_);
			for (int i = 0; i < a.rows_; ++i) {
				r[i] = a[i] - b[i];
			}

			return r;
		}

		/**
		 * Calculates the dot-product of 2 vectors
		 * @param a first vector
		 * @param b second vector
		 * @return dot product of 2 vectors
		 */
		friend T operator *(const ColumnVector<T> &a, const ColumnVector<T> &b) {
			if (a.rows_ != b.rows_) {
				throw dimensions_error("Vectors should have the same number of rows");
			}
			T r = 0;

			for (int i = 0; i < a.rows_; ++i) {
				r += a[i] * b[i];
			}

			return r;
		}

		/**
		 * Reads ColumnVector from stream
		 * @param is stream to read from
		 * @param v ColumnVector to read to
		 * @return passed in stream
		 */
		friend std::istream &operator >>(std::istream &is, ColumnVector<T> &v) {
			for (dim_t i = 0; i < v.rows_; ++i) {
				is >> v[i];
			}

			return is;
		}

		/**
		 * Prints ColumnVector to stream
		 * @param os stream to print to
		 * @param v ColumnVector to print
		 * @return passed in stream
		 */
		friend std::ostream &operator <<(std::ostream &os, const ColumnVector<T> &v) {
			auto precision = os.precision();

			for (dim_t i = 0; i < v.rows_; ++i) {
				const auto value = v[i];
				os << (la::is_less_than_10_nn(value, precision) ? 0 : value) << std::endl;
			}

			return os;
		}
	};

	using mati = Matrix<int>;
	using matd = Matrix<double>;
	using cvecd = ColumnVector<double>;
}

int main() {
	la::d_epsilon = 10e-9;
	std::istream &in = std::cin;
	std::ofstream out("points.dat");
	out << std::fixed << std::setprecision(4);

	int dataset_size, degree;
	in >> dataset_size;
	double d[dataset_size];
	la::cvecd b(dataset_size);
	for (int i = 0; i < dataset_size; ++i) {
		double x, y;
		in >> x >> y, out << x << ' ' << y << std::endl;
		d[i] = x, b[i] = y;
	}

	in >> degree;
	la::matd a(dataset_size, degree + 1);
	for (int i = 0; i < a.rows(); ++i) {
		for (int j = 0; j < a.cols(); ++j) {
			a[i][j] = std::pow(d[i], j);
		}
	}

	la::matd a_t = a.transpose();
	b = (a_t * a).inverse() * (a_t * b);

	FILE *gnuplot = popen("gnuplot -persistent", "w");
	fprintf(gnuplot, "set title 'Least squares approximation'\n");
	fprintf(gnuplot, "plot 'points.dat' with points title 'Data', \\\n");
	for (int i = 0; i < b.rows(); ++i) {
		fprintf(gnuplot, "%f*x**%d", b[i], i);
		if (i < b.rows() - 1) fprintf(gnuplot, "+");
	}
	fprintf(gnuplot, " title 'Approximation'\n");
	pclose(gnuplot);
	return 0;
}

template<typename T>
requires std::is_arithmetic_v<T>
la::ColumnVector<T> la::Matrix<T>::operator *(const ColumnVector<T> &vector) const {
	if (cols_ != vector.rows()) {
		throw la::dimensions_error();
	}

	ColumnVector<T> result(rows_);

	for (dim_t i = 0; i < rows_; ++i) {
		for (dim_t j = 0; j < cols_; ++j) {
			result[i] += vector[j] * at(i, j);
		}
	}

	return result;
}