Numerical Explorations

// Algorithms from "A Course in Computational
// Algebraic Number Theory" by Henri Cohen
// Implemented by James Pate Williams, Jr.
// Copyright (c) 2023 All Rights Reserved

#pragma once
#include "pch.h"

template<class T> class Matrix
{
public:
	size_t m, n;
	T** data;

	Matrix() { m = 0; n = 0; data = NULL; };
	Matrix(size_t m, size_t n)
	{
		this->m = m;
		this->n = n;
		data = new T*[m];

		if (data == NULL)
			exit(-300);

		for (size_t i = 0; i < m; i++)
		{
			data[i] = new T[n];

			if (data[i] == NULL)
				exit(-301);
		}
	};
	void OutputMatrix(
		fstream& outs, char fill, int precision, int width)
	{
		for (size_t i = 0; i < m; i++)
		{
			for (size_t j = 0; j < n; j++)
			{
				outs << setfill(fill) << setprecision(precision);
				outs << setw(width) << data[i][j] << '\t';
			}

			outs << endl;
		}
	};
};

template<class T> class Vector
{
public:
	size_t n;
	T* data;

	Vector() { n = 0; data = NULL; };
	Vector(size_t n)
	{
		this->n = n;
		data = new T[n];
	};
	void OutputVector(
		fstream& outs, char fill, int precision, int width)
	{
		for (size_t i = 0; i < n; i++)
		{
			outs << setfill(fill) << setprecision(precision);
			outs << setw(width) << data[i] << '\t';
		}

		outs << endl;
	};
};

class LinearAlgebra
{
public:
	bool initialized;
	size_t m, n;
	Matrix<double> M;
	Vector<double> B;

	LinearAlgebra()
	{ 
		initialized = false;
		m = 0; n = 0;
		M.data = NULL;
		B.data = NULL;
	};
	LinearAlgebra(size_t m, size_t n)
	{
		initialized = false;
		this->m = m;
		this->n = n;
		this->M.m = m;
		this->M.n = n;
		this->B.n = n;
		this->M.data = new double*[m];
		this->B.data = new double[n];

		if (M.data != NULL)
		{
			for (size_t i = 0; i < m; i++)
			{
				this->M.data[i] = new double[n];

				for (size_t j = 0; j < n; j++)
					this->M.data[i][j] = 0;
			}
		}

		if (B.data != NULL)
		{
			this->B.data = new double[n];

			for (size_t i = 0; i < n; i++)
				this->B.data[i] = 0;
		}

		initialized = this->B.data != NULL && this->M.data != NULL;
	};
	LinearAlgebra(
		size_t m, size_t n,
		double** M,
		double* B)
	{
		this->m = m;
		this->n = n;
		this->M.m = m;
		this->M.n = n;
		this->M.data = new double*[m];
		this->B.data = new double[n];

		if (M != NULL)
		{
			for (size_t i = 0; i < m; i++)
			{
				this->M.data[i] = new double[n];

				for (size_t j = 0; j < n; j++)
					this->M.data[i][j] = M[i][j];
			}
		}

		if (B != NULL)
		{
			this->B.data = new double[n];

			for (size_t i = 0; i < m; i++)
				this->B.data[i] = B[i];
		}

		initialized = this->B.data != NULL && this->M.data != NULL;
	}
	~LinearAlgebra()
	{
		M.m = m;
		M.n = n;
		B.n = n;

		if (B.data != NULL)
			delete[] B.data;

		for (size_t i = 0; i < m; i++)
			if (M.data[i] != NULL)
				delete[] M.data[i];

		if (M.data != NULL)
			delete[] M.data;
	}
	double DblDeterminant(bool failure);
	Vector<double> DblGaussianElimination(
		bool& failure);
	// The following four methods are from the
	// textbook "Elementary Numerical Analysis
	// An Algorithmic Approach" by S. D. Conte
	// and Carl de Boor Translated from the
	// original FORTRAN by James Pate Williams, Jr.
	// Copyright (c) 2023 All Rights Reserved
	bool DblGaussianFactor(
		Vector<int>& pivot);
	bool DblGaussianSolution(
		Vector<double>& x,
		Vector<int>& pivot);
	bool DblSubstitution(
		Vector<double>& x,
		Vector<int>& pivot);
	bool DblInverse(
		Matrix<double>& A,
		Vector<int>& pivot);
	// Henri Cohen Algorithm 2.2.7
	void DblCharPolyAndAdjoint(
		Matrix<double>& adjoint,
		Vector<double>& a);
	// Henri Cohen Algorithm 2.3.1
	void DblMatrixKernel(
		Matrix<double>& X, size_t& r);
	// Henri Cohen Algorithm 2.3.1
	void DblMatrixImage(
		Matrix<double>& N, size_t& rank);
};

// pch.h: This is a precompiled header file.
// Files listed below are compiled only once, improving build performance for future builds.
// This also affects IntelliSense performance, including code completion and many code browsing features.
// However, files listed here are ALL re-compiled if any one of them is updated between builds.
// Do not add files here that you will be updating frequently as this negates the performance advantage.

#ifndef PCH_H
#define PCH_H
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
using namespace std;
#endif //PCH_H

#include "pch.h"
#include "LinearAlgebra.h"

double LinearAlgebra::DblDeterminant(
    bool failure)
{
    double deter = 1;
    Vector<int> pivot(n);

    if (!initialized || m != n)
    {
        failure = true;
        return 0.0;
    }

    if (!DblGaussianFactor(pivot))
    {
        failure = true;
        return 0.0;
    }

    for (size_t i = 0; i < n; i++)
        deter *= M.data[i][i];

    return deter;
}

Vector<double> LinearAlgebra::DblGaussianElimination(
    bool& failure)
{
    double* C = new double[m];
    Vector<double> X(n);

    X.data = new double[n];

    if (X.data == NULL)
        exit(-200);

    if (!initialized)
    {
        failure = true;
        delete[] C;
        return X;
    }
    
    for (size_t i = 0; i < m; i++)
        C[i] = -1;

    failure = false;

    for (size_t j = 0; j < n; j++)
    {
        bool found = false;
        size_t i = j;

        while (i < n && !found)
        {
            if (M.data[i][j] != 0)
                found = true;
            else
                i++;
        }

        if (!found)
        {
            failure = true;
            break;
        }

        if (i > j)
        {
            for (size_t l = j; l < n; l++)
            {
                double t = M.data[i][l];
                M.data[i][l] = M.data[j][l];
                M.data[j][l] = t;
                t = B.data[i];
                B.data[i] = B.data[j];
                B.data[j] = t;
            }
        }

        double d = 1.0 / M.data[j][j];

        for (size_t k = j + 1; k < n; k++)
            C[k] = d * M.data[k][j];

        for (size_t k = j + 1; k < n; k++)
        {
            for (size_t l = j + 1; l < n; l++)
                M.data[k][l] = M.data[k][l] - C[k] * M.data[j][l];
            
            B.data[k] = B.data[k] - C[k] * B.data[j];
        }
    }

    for (int i = (int)n - 1; i >= 0; i--)
    {
        double sum = 0;

        for (size_t j = i + 1; j < n; j++)
            sum += M.data[i][j] * X.data[j];

        X.data[i] = (B.data[i] - sum) / M.data[i][i];
    }

    delete[] C;
    return X;
}

bool LinearAlgebra::DblGaussianFactor(
    Vector<int>& pivot)
    // returns false if matrix is singular
{
    Vector<double> d(n);
    double awikod, col_max, ratio, row_max, temp;
    int flag = 1;
    size_t i_star, itemp;

    for (size_t i = 0; i < n; i++)
    {
        pivot.data[i] = i;
        row_max = 0;
        for (size_t j = 0; j < n; j++)
            row_max = max(row_max, abs(M.data[i][j]));
        if (row_max == 0)
        {
            flag = 0;
            row_max = 1;
        }
        d.data[i] = row_max;
    }
    if (n <= 1) return flag != 0;
    // factorization
    for (size_t k = 0; k < n - 1; k++)
    {
        // determine pivot row the row i_star
        col_max = abs(M.data[k][k]) / d.data[k];
        i_star = k;
        for (size_t i = k + 1; i < n; i++)
        {
            awikod = abs(M.data[i][k]) / d.data[i];
            if (awikod > col_max)
            {
                col_max = awikod;
                i_star = i;
            }
        }
        if (col_max == 0)
            flag = 0;
        else
        {
            if (i_star > k)
            {
                // make k the pivot row by
                // interchanging with i_star
                flag *= -1;
                itemp = pivot.data[i_star];
                pivot.data[i_star] = pivot.data[k];
                pivot.data[k] = itemp;
                temp = d.data[i_star];
                d.data[i_star] = d.data[k];
                d.data[k] = temp;
                for (size_t j = 0; j < n; j++)
                {
                    temp = M.data[i_star][j];
                    M.data[i_star][j] = M.data[k][j];
                    M.data[k][j] = temp;
                }
            }
            // eliminate x[k]
            for (size_t i = k + 1; i < n; i++)
            {
                M.data[i][k] /= M.data[k][k];
                ratio = M.data[i][k];
                for (size_t j = k + 1; j < n; j++)
                    M.data[i][j] -= ratio * M.data[k][j];
            }
        }

        if (M.data[n - 1][n - 1] == 0) flag = 0;
    }

    if (flag == 0)
        return false;

    return true;
}

bool LinearAlgebra::DblGaussianSolution(
    Vector<double>& x,
    Vector<int>& pivot)
{
    if (!DblGaussianFactor(pivot))
        return false;

    return DblSubstitution(x, pivot);
}

bool LinearAlgebra::DblSubstitution(
    Vector<double>& x,
    Vector<int>& pivot)
{
    double sum;
    size_t j, n1 = n - 1;

    if (n == 1)
    {
        x.data[0] = B.data[0] / M.data[0][0];
        return true;
    }

    // forward substitution

    x.data[0] = B.data[pivot.data[0]];

    for (int i = 1; i < (int)n; i++)
    {
        for (j = 0, sum = 0; j < (size_t)i; j++)
            sum += M.data[i][j] * x.data[j];

        x.data[i] = B.data[pivot.data[i]] - sum;
    }

    // backward substitution

    x.data[n1] /= M.data[n1][n1];

    for (int i = n - 2; i >= 0; i--)
    {
        double sum = 0.0;

        for (j = i + 1; j < n; j++)
            sum += M.data[i][j] * x.data[j];

        x.data[i] = (x.data[i] - sum) / M.data[i][i];
    }

    return true;
}

bool LinearAlgebra::DblInverse(
    Matrix<double>& A,
    Vector<int>& pivot)
{
    Vector<double> x(n);

    if (!DblGaussianFactor(pivot))
        return false;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            B.data[j] = 0;
    }
    
    for (size_t i = 0; i < n; i++)
    {
        B.data[i] = 1;

        if (!DblSubstitution(x, pivot))
            return false;

        B.data[i] = 0;

        for (size_t j = 0; j < n; j++)
           A.data[i][j] = x.data[pivot.data[j]];
    }

    return true;
}

void LinearAlgebra::DblCharPolyAndAdjoint(
    Matrix<double>& adjoint,
    Vector<double>& a)
{
    Matrix<double> C(n, n);
    Matrix<double> I(n, n);

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            C.data[i][j] = I.data[i][j] = 0;
    }

    for (size_t i = 0; i < n; i++)
        C.data[i][i] = I.data[i][i] = 1;

    a.data[0] = 1;

    for (size_t i = 1; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
            {
                double sum = 0.0;

                for (size_t l = 0; l < n; l++)
                    sum += M.data[j][l] * C.data[l][k];

                C.data[j][k] = sum;
            }
        }

        double tr = 0.0;

        for (size_t j = 0; j < n; j++)
            tr += C.data[j][j];

        a.data[i] = -tr / i;

        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
                C.data[j][k] += a.data[i] * I.data[j][k];
        }
    }

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            double sum = 0.0;

            for (size_t k = 0; k < n; k++)
                sum += M.data[i][k] * C.data[k][j];

            C.data[i][j] = sum;
        }
    }

    double trace = 0.0;

    for (size_t i = 0; i < n; i++)
        trace += C.data[i][i];

    trace /= n;
    a.data[n - 1] = -trace;

    double factor = 1.0;

    if ((n - 1) % 2 != 0)
        factor = -1.0;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            adjoint.data[i][j] = factor * C.data[i][j];
    }
}

void LinearAlgebra::DblMatrixKernel(Matrix<double>& X, size_t& r)
{
    double D = 0.0;
    Vector<int> c(m);
    Vector<int> d(n);

    r = 0;

    for (size_t i = 0; i < m; i++)
        c.data[i] = -1;

    size_t j, k = 1;
Step2:
    for (j = 0; j < m; j++)
    {
        if (M.data[j][k] != 0 && c.data[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d.data[k] = 0;
        goto Step4;
    }

    D = -1.0 / M.data[j][k];

    M.data[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        M.data[j][s] = D * M.data[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = M.data[i][k];
                M.data[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            M.data[i][s] += D * M.data[j][s];
        }
    }

    c.data[j] = k;
    d.data[k] = j;
Step4:
    if (k < n - 1)
    {
        k++;
        goto Step2;
    }

    X.n = n;

    if (r != 0)
    {
        for (k = 0; k < n; k++)
        {
            if (d.data[k] == 0)
            {
                for (size_t i = 0; i < n; i++)
                {
                    if (d.data[i] > 0)
                        X.data[k][i] = M.data[d.data[i]][k];
                    else if (i == k)
                        X.data[k][i] = 1;
                    else
                        X.data[k][i] = 0;
                }
            }
        }
    }
}

void LinearAlgebra::DblMatrixImage(
    Matrix<double>& N, size_t& rank)
{
    double D = 0.0;
    size_t r = 0;
    Matrix<double> copyM(m, n);
    Vector<int> c(m);
    Vector<int> d(n);

    for (size_t i = 0; i < m; i++)
        c.data[i] = -1;

    size_t j, k = 1;
    N = copyM = M;
Step2:
    for (j = 0; j < m; j++)
    {
        if (copyM.data[j][k] != 0 && c.data[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d.data[k] = 0;
        goto Step4;
    }

    D = -1.0 / copyM.data[j][k];

    copyM.data[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        copyM.data[j][s] = D * copyM.data[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = copyM.data[i][k];
                copyM.data[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            copyM.data[i][s] += D * copyM.data[j][s];
        }
    }

    c.data[j] = k;
    d.data[k] = j;
Step4:
    if (k < n - 1)
    {
        k++;
        goto Step2;
    }

    rank = n - r;

    for (j = 0; j < m; j++)
    {
        if (c.data[j] != 0)
        {
            for (size_t i = 0; i < m; i++)
            {
                N.data[i][c.data[j]] = M.data[i][c.data[j]];
            }
        }
    }
}

/*
** Cohen's linear algebra test program
** Implemented by James Pate Williams, Jr.
** Copyright (c) 2023 All Rights Reserved
*/

#include "pch.h"
#include "LinearAlgebra.h"

double GetDblNumber(fstream& inps)
{
    char ch = inps.get();
    string numberStr;

    while (ch == ' ' || ch == '\t' || ch == '\r' || ch == '\n')
        ch = inps.get();

    while (ch == '+' || ch == '-' || ch == '.' ||
        ch >= '0' && ch <= '9')
    {
        numberStr += ch;
        ch = inps.get();
    }

    double x = atof(numberStr.c_str());
    return x;
}

int GetIntNumber(fstream& inps)
{
    char ch = inps.get();
    string numberStr;

    while (ch == ' ' || ch == '\t' || ch == '\r' || ch == '\n')
        ch = inps.get();

    while (ch == '+' || ch == '-' || ch >= '0' && ch <= '9')
    {
        numberStr += ch;
        ch = inps.get();
    }

    int x = atoi(numberStr.c_str());
    return x;
}

int main()
{
    fstream inps;

    inps.open("CLATestFile.txt", fstream::in);
    
    if (inps.fail())
    {
        cout << "Input file opening error!" << endl;
        return -1;
    }

    fstream outs;

    outs.open("CLAResuFile.txt", fstream::out | fstream::ate);

    if (outs.fail())
    {
        cout << "Output file opening error!" << endl;
        return -2;
    }

    size_t m, n;
    
    while (!inps.eof())
    {
        char buffer[256] = { '\0' };

        inps.getline(buffer, 256);
        m = GetIntNumber(inps);

        if (inps.eof())
            return 0;

        if (m < 1)
        {
            cout << "The number of rows must be >= 1" << endl;
            return -100;
        }

        n = GetIntNumber(inps);

        if (n < 1)
        {
            cout << "The number of colums must be >= 1" << endl;
            return -101;
        }

        LinearAlgebra la(m, n);
        Matrix<double> copyM(m, n);
        Vector<double> copyB(n);

        inps.getline(buffer, 256);

        for (size_t i = 0; i < m; i++)
        {
            for (size_t j = 0; j < n; j++)
            {
                double x = GetDblNumber(inps);

                la.M.data[i][j] = x;
                copyM.data[i][j] = x;
            }
        }

        inps.getline(buffer, 256);

        for (size_t i = 0; i < n; i++)
        {
            la.B.data[i] = GetDblNumber(inps);
            copyB.data[i] = la.B.data[i];
        }

        bool failure = false;
        Vector<double> X = la.DblGaussianElimination(failure);

        if (!failure)
        {
            outs << "The 1st solution of the linear system of equations:" << endl;
            X.OutputVector(outs, ' ', 5, 8);
        }
        else
        {
            cout << "Cohen Gaussian elimination failure!" << endl;
            exit(-102);
        }

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        Matrix<double> A(n, n);
        Vector<int> pivot(n);

        if (!la.DblGaussianSolution(X, pivot))
            exit(-103);

        outs << "The 2nd solution of the linear system of equations:" << endl;

        X.OutputVector(outs, ' ', 5, 8);

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        double deter = la.DblDeterminant(failure);

        outs << "The determinant of the linear system of equations:" << endl;
        outs << deter << endl;

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        outs << "The inverse of the linear system of equations:" << endl;

        if (!la.DblInverse(A, pivot))
        {
            cout << "Conte Gaussian inverse matrix failure!" << endl;
            exit(-104);
        }

        else
            A.OutputMatrix(outs, ' ', 5, 8);

        Matrix<double> adjoint(n, n);
        Vector<double> a(n);

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        la.DblCharPolyAndAdjoint(adjoint, a);
        outs << "The adjoint of the linear system of equations:" << endl;
        adjoint.OutputMatrix(outs, ' ', 5, 8);
        outs << "The characteristic polynomial of the linear system of equations:" << endl;
        a.OutputVector(outs, ' ', 5, 8);

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        Matrix<double> kernel(m, n);
        size_t r = 0;

        la.DblMatrixKernel(kernel, r);

        if (r > 0)
        {
            outs << "The kernel of the matrix: " << endl;
            kernel.OutputMatrix(outs, ' ', 5, 8);
            outs << "Dimension of the kernel: " << r << endl;
        }

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        Matrix<double> N(m, n);
        size_t rank;

        la.DblMatrixImage(N, rank);

        if (rank > 0)
        {
            outs << "The image of the matrix: " << endl;
            N.OutputMatrix(outs, ' ', 5, 8);
            outs << "Rank = " << rank << endl;
        }
    }

    inps.close();
    outs.close();
}