Category: Matrix Algebra

Latest Factoring Results (c) February 4, 2024, by James Pate Williams, Jr.

I am testing two factoring algorithms: Pollard-Shor-Williams’s method, a home-grown version of the venerable Pollard rho algorithm and Pollard’s factoring with cubic integers. The second recipe is from “The Development of the Number Field Sieve” edited by Arjen K. Lenstra and Hendrik W. Lenstra, Jr. I use the 20-digit test number, 2 ^ 66 + 2 = 73786976294838206466. My method is very fast with this number as shown below:

2^66+2
73786976294838206466 20
Pseudo-Random Number Generator Seed = 1
Number of Tasks = 1
2 p 1
3 p 1
11 p 2
131 p 3
2731 p 4
409891 p 6
7623851 p 7
Elapsed hrs:min:sec.MS = 00:00:00.652
Function Evaluations = 1995

The Pollard factoring with cubic integers takes a long time but is capable of factoring much larger numbers. The results of the full factorization of my test number are:

== Data Menu ==
1 Simple Number
2 Fibonacci Sequence Number
3 Lucas Sequence Number
4 Exit
Enter option (1 – 4): 1
Enter a number to be factored: 2^66+2
Enter a random number generator seed: 1
== Factoring Menu ==
1 Lenstra’s ECM
2 Lenstra’s Pollard-Rho
3 Pollard’s Factoring with Cubic Integers
Option (1 – 3): 3
73786976294838206466

Enter a lower bound : -50000
Enter a upper bound : +50000
Enter b lower bound : +1
Enter b upper bound : +50000
Enter maximum kernels : 1024
Enter algebraic prime count: 300
Enter rational prime count: 300
Enter lo large prime bound: 10000
Enter hi large prime bound: 11000
Numbers sieved = 452239293
Successes 0 gcd(a, b) is 1 = 274929004
Successes 1 rational smooth a + b * r = 181838258
Successes 2 long long smooth = 68959
Successes 3 kernels tested = 535
2 p # digits 1
3 p # digits 1
11 p # digits 2
131 p # digits 3
2731 p # digits 4
409891 p # digits 6
7623851 p # digits 7
Runtime (s) = 36383.696000

It took over ten hours to fully factor, 2 ^ 66 + 2. I am currently attempting to factor the Thirteenth Fermat number which is 2 ^ (2 ^ 13) + 1 = 2 ^ 8192 + 1. The number has 2,467 decimal digits. I am using 399 algebraic prime numbers, 600 rational prime numbers, and 316 “large prime numbers” (primes between 12,000 and 15,000). I have to find the kernels of a 1316 by 1315 matrix. I am trying the factorization using a maximum of 8192 kernels. I suspect this computation will take about a week on my desktop workstation. There is no guarantee that I will find a non-trivial factor of 2 * (2731 ^ 3) + 2.

Fast and Slow Matrix Multiplication by James Pate Williams, Jr. © January 13, 2024, All Applicable Rights Reserved

Note: For comparison with a paper on matrix operations, we use column major order also known as column major echelon order. Column major order is used in Fortran and row major order is utilized in C computations.

Click to access Lecture4.pdf

Here is the supposedly slow C source code:

void ColMajorSlowMatMul(double** a, double** b, double** c, int n)
{
    for (int j = 1; j <= n; j++)
        for (int k = 1; k <= n; k++)
            c[j][k] = 0.0;

    for (int j = 1; j <= n; j++)
    {
        for (int i = 1; i <= n; i++)
        {
            double t = 0;

            for (int k = 1; k <= n; k++)
                t += a[i][k] * b[k][j];

            c[i][j] = t;
        }
    }
}

And now we present the allegedly fast multiplication code:

void ColMajorFastMatMul(double** a, double** b, double** c, int n)
{
    for (int j = 1; j <= n; j++)
        for (int k = 1; k <= n; k++)
            c[j][k] = 0.0;
    
    for (int j = 1; j <= n; j++)
    {
        for (int k = 1; k <= n; k++)
        {
            double t = b[k][j];

            for (int i = 1; i <= n; i++)
                c[i][j] += a[i][k] * t;
        }
    }
}

Finally, we have another algorithm for matrix multiplication:

void ColMajorSemiMatMul(double** a, double** b, double** c, int n)
{
    for (int j = 1; j <= n; j++)
        for (int k = 1; k <= n; k++)
            c[j][k] = 0.0;

    for (int j = 1; j <= n; j++)
    {
        for (int i = 1; i <= n; i++)
        {
            for (int k = 1; k <= n; k++)
                c[i][j] += a[i][k] * b[k][j];
        }
    }
}

We performed the following experiments.



slow runtime in microseconds: 2736

semi runtime in microseconds: 4001

fast runtime in microseconds: 5293

n * n = 10000

equal = 10000

n * n = 10000

equal = 10000

slow runtime in microseconds: 21962

semi runtime in microseconds: 29585

fast runtime in microseconds: 28701

n * n = 40000

equal = 40000

n * n = 40000

equal = 40000

slow runtime in microseconds: 67256

semi runtime in microseconds: 101554

fast runtime in microseconds: 107969

n * n = 90000

equal = 90000

n * n = 90000

equal = 90000

slow runtime in microseconds: 183839

semi runtime in microseconds: 268616

fast runtime in microseconds: 244832

n * n = 160000

equal = 160000

n * n = 160000

equal = 160000

slow runtime in microseconds: 428263

semi runtime in microseconds: 638721

fast runtime in microseconds: 650535

n * n = 250000

equal = 250000

n * n = 250000

equal = 250000



C:\Users\james\source\repos\FastMatrixMultiplication\Debug\FastMatrixMultiplication.exe (process 29552) exited with code 0.

Press any key to close this window . . .































// FastMatrixMultiplication.cpp : This file contains the 'main' function. Program execution begins and ends there.
// https://www-users.york.ac.uk/~mijp1/teaching/grad_HPC_for_MatSci/Lecture4.pdf
// https://stackoverflow.com/questions/25483620/how-to-measure-running-time-of-specific-function-in-c-very-accurate

#include <stdlib.h>
#include <chrono>
#include <iostream>
using namespace std;

typedef chrono::high_resolution_clock Clock;

void ColMajorFastMatMul(double** a, double** b, double** c, int n)
{
    for (int j = 1; j <= n; j++)
        for (int k = 1; k <= n; k++)
            c[j][k] = 0.0;
    
    for (int j = 1; j <= n; j++)
    {
        for (int k = 1; k <= n; k++)
        {
            double t = b[k][j];

            for (int i = 1; i <= n; i++)
                c[i][j] += a[i][k] * t;
        }
    }
}

void ColMajorSemiMatMul(double** a, double** b, double** c, int n)
{
    for (int j = 1; j <= n; j++)
        for (int k = 1; k <= n; k++)
            c[j][k] = 0.0;

    for (int j = 1; j <= n; j++)
    {
        for (int i = 1; i <= n; i++)
        {
            for (int k = 1; k <= n; k++)
                c[i][j] += a[i][k] * b[k][j];
        }
    }
}

void ColMajorSlowMatMul(double** a, double** b, double** c, int n)
{
    for (int j = 1; j <= n; j++)
        for (int k = 1; k <= n; k++)
            c[j][k] = 0.0;

    for (int j = 1; j <= n; j++)
    {
        for (int i = 1; i <= n; i++)
        {
            double t = 0;

            for (int k = 1; k <= n; k++)
                t += a[i][k] * b[k][j];

            c[i][j] = t;
        }
    }
}

void GenerateMatrix(double** a, int n, int seed)
{
    srand(seed);

    for (int j = 1; j <= n; j++)
    {
        for (int k = 1; k <= n; k++)
        {
            a[j][k] = (double)rand() / RAND_MAX;
        }
    }
}

int main()
{
    for (int n = 100; n <= 500; n += 100)
    {
        double** a = new double* [n + 1];
        double** b = new double* [n + 1];
        double** c = new double* [n + 1];
        double** d = new double* [n + 1];
        double** e = new double* [n + 1];

        for (int j = 0; j < n + 1; j++)
        {
            a[j] = new double[n + 1];
            b[j] = new double[n + 1];
            c[j] = new double[n + 1];
            d[j] = new double[n + 1];
            e[j] = new double[n + 1];
        }

        GenerateMatrix(a, n, 1);
        GenerateMatrix(b, n, 2);

        auto clock1 = Clock::now();
        ColMajorSlowMatMul(a, b, c, n);
        auto clock2 = Clock::now();
        long long microseconds1 = chrono::duration_cast<chrono::microseconds>
            (clock2 - clock1).count();
    
        auto clock3 = Clock::now();
        ColMajorSemiMatMul(a, b, d, n);
        auto clock4 = Clock::now();
        long long microseconds2 = chrono::duration_cast<chrono::microseconds>
            (clock4 - clock3).count();

        auto clock5 = Clock::now();
        ColMajorFastMatMul(a, b, e, n);
        auto clock6 = Clock::now();
        long long microseconds3 = chrono::duration_cast<chrono::microseconds>
            (clock6 - clock5).count();

        cout << "slow runtime in microseconds: " << microseconds1 << endl;
        cout << "semi runtime in microseconds: " << microseconds2 << endl;
        cout << "fast runtime in microseconds: " << microseconds3 << endl;

        long long equal = 0;

        for (int j = 1; j <= n; j++)
            for (int k = 1; k <= n; k++)
                if (c[j][k] == d[j][k])
                    equal++;

        cout << "n * n = " << n * n << endl;
        cout << "equal = " << equal << endl;

        equal = 0;

        for (int j = 1; j <= n; j++)
            for (int k = 1; k <= n; k++)
                if (c[j][k] == e[j][k])
                    equal++;

        cout << "n * n = " << n * n << endl;
        cout << "equal = " << equal << endl;

        for (int j = 0; j < n + 1; j++)
        {
            delete[] a[j];
            delete[] b[j];
            delete[] c[j];
            delete[] d[j];
            delete[] e[j];
        }

        delete[] a;
        delete[] b;
        delete[] c;
        delete[] d;
        delete[] e;
    }
}

Comparison of Linear Systems Applications in C# and C++ by James Pate Williams, Jr.

Back in 2017 I created a C# application that implemented the direct methods: Cholesky decomposition, Gaussian elimination with partial pivoting, LU decomposition, and simple Gaussian elimination. The classical iterative methods Gauss-Seidel, Jacobi, Successive Overrelaxation, and gradient descent were also implemented along with the modern iterative methods: conjugate gradient descent and Modified Richardson’s method. Recently I translated the C# code to C++. I used the following test matrices: Cauchy, Lehmer, Pascal, and other. Below are some results. As is apparent the C++ runtimes are faster than the C# execution times.

linearsystems-resultsDownload

Richardson Method Translated from C Source Code to C# by James Pate Williams, Jr.

The Richardson Method is called before eliminating the system of linear equations.

Added elimination results from a vanilla C program and a C# application.

These results are not in total agreement with H. T. Lau’s results.

Finite Difference Method for Solving Second Order Ordinary Differential Equations by James Pate Williams, Jr.

My reference is Numerical Analysis: An Algorithmic Approach 3rd Edition by S. D. Conte and Carl de Boor Chapter 9.1.

Determinant Calculations: Bareiss Algorithm and Gaussian Elimination by James Pate Williams, Jr.

Below are three screenshots of two methods of calculating the determinant of a matrix, namely the Bareiss Algorithm and Gaussian Elimination:

using System;
using System.Diagnostics;
using System.Windows.Forms;

namespace MatrixInverseComparison
{
    public partial class MainForm : Form
    {
        public MainForm()
        {
            InitializeComponent();
        }
        static private string FormatNumber(double x)
        {
            string result = string.Empty;

            if (x > 0)
                result += x.ToString("F5").PadLeft(10);
            else
                result += x.ToString("F5").PadLeft(10);

            return result;
        }

        private void MultiplyPrintMatricies(
            double[,] A, double[,] B, int n)
        {
            double[,] I = new double[n, n];

            textBox1.Text += "Matrix Product:\r\n";

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

                    for (int k = 0; k < n; k++)
                        sum += A[i, k] * B[k, j];

                    I[i, j] = sum;
                }
            }

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    textBox1.Text += FormatNumber(I[i, j]) + " ";
                }

                textBox1.Text += "\r\n";
            }
        }

        private void PrintMatrix(string title, double[,] A, int n)
        {
            textBox1.Text += title + "\r\n";

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    textBox1.Text += FormatNumber(A[i, j]) + " ";
                }

                textBox1.Text += "\r\n";
            }
        }

        private void Compute(double[,] MI, int n)
        {
            double determinantGE = 1;
            double[] b = new double[n];
            double[] x = new double[n];
            double[,] MB1 = new double[n, n];
            double[,] MB2 = new double[n, n];
            double[,] MG = new double[n, n];
            double[,] MS = new double[n, n];
            double[,] IG = new double[n, n];
            double[,] IS = new double[n, n];
            int[] pivot = new int[n];
            Stopwatch sw = new();
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    MB1[i, j] = MB2[i, j] = MG[i, j] = MS[i, j] = MI[i, j];
            }
            PrintMatrix("Initial Matrix: ", MI, n);
            sw.Start();
            int flag = DirectMethods.Factor(MG, n, pivot);
            if (flag != 0)
            {
                for (int i = 0; i < n; i++)
                    determinantGE *= MG[i, i];
                determinantGE *= flag;
            }
            sw.Stop();
            PrintMatrix("Gaussian Elimination Final:", MG, n);
            for (int i = 0; i < n; i++)
                b[i] = 0;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    b[j] = 1.0;
                    DirectMethods.Substitute(MG, b, x, n, pivot);
                    for (int k = 0; k < n; k++)
                        IG[k, j] = x[k];
                    b[j] = 0.0;
                }
            }
            PrintMatrix("Gaussian Elimination Inverse:", IG, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinantGE) + "\r\n";
            MultiplyPrintMatricies(IG, MI, n);
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
            sw.Start();
            double determinant1 = BareissAlgorithm.Determinant1(MB1, n);
            sw.Stop();
            PrintMatrix("Bareiss Algorithm Final 1:", MB1, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinant1) + "\r\n";
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
            sw.Start();
            double determinant2 = BareissAlgorithm.Determinant2(MB2, n);
            sw.Stop();
            PrintMatrix("Bareiss Algorithm Final 2:", MB2, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinant2) + "\r\n";
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
        }
        private void button1_Click(object sender, EventArgs e)
        {
            int n = (int)numericUpDown1.Value;
            int seed = (int)numericUpDown2.Value;
            double[,] A = new double[n, n];
            Random random = new Random(seed);

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    double q = n * random.NextDouble();

                    while (q == 0.0)
                        q = n * random.NextDouble();

                    A[i, j] = q;
                }
            }
            
            Compute(A, n);
        }

        private void button2_Click(object sender, EventArgs e)
        {
            textBox1.Text = string.Empty;
        }
    }
}

using System;

namespace MatrixInverseComparison
{
    public class DirectMethods
    {
        // Substitute and Factor translated from FORTRAN 77
        // source code found in "Elementary Numerical Analysis:
        // An Algorithmic Approach" by S. D. Conte and Carl de
        // Boor. Translator: James Pate Williams, Jr. (c)
        // August 14 - 17, 2023
        static public void Substitute(
            double[,] w,
            double[] b,
            double[] x,
            int n,
            int[] pivot)
        {
            double sum;
            int i, j, n1 = n - 1;

            if (n == 1)
            {
                x[0] = b[0] / w[0, 0];
                return;
            }

            // forward substitution
            x[0] = b[pivot[0]];

            for (i = 0; i < n; i++)
            {
                for (j = 0, sum = 0.0; j < i; j++)
                    sum += w[i, j] * x[j];
                x[i] = b[pivot[i]] - sum;
            }

            // backward substitution
            x[n1] /= w[n1, n1];

            for (i = n - 2; i >= 0; i--)
            {
                for (j = i + 1, sum = 0.0; j < n; j++)
                    sum += w[i, j] * x[j];
                x[i] = (x[i] - sum) / w[i, i];
            }
        }

        // Factor returns +1 if an even number of exchanges
        // Factor returns -1 if an odd  number of exchanges
        // Factor retrurn 0  if matrix is singular
        static public int Factor(
            double[,] w, int n, int[] pivot)
        // returns 0 if matrix is singular
        {
            double awikod, col_max, ratio, row_max, temp;
            double[] d = new double[n];
            int flag = 1, i, i_star, j, k;

            for (i = 0; i < n; i++)
            {
                pivot[i] = i;
                row_max = 0;
                for (j = 0; j < n; j++)
                    row_max = Math.Max(row_max, Math.Abs(w[i, j]));
                if (row_max == 0)
                {
                    flag = 0;
                    row_max = 1;
                }
                d[i] = row_max;
            }
            if (n <= 1)
                return flag;
            // factorization
            for (k = 0; k < n - 1; k++)
            {
                // determine pivot row the row i_star
                col_max = Math.Abs(w[k, k]) / d[k];
                i_star = k;
                for (i = k + 1; i < n; i++)
                {
                    awikod = Math.Abs(w[i, k]) / d[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;
                        i = pivot[i_star];
                        pivot[i_star] = pivot[k];
                        pivot[k] = i;
                        temp = d[i_star];
                        d[i_star] = d[k];
                        d[k] = temp;
                        for (j = 0; j < n; j++)
                        {
                            temp = w[i_star, j];
                            w[i_star, j] = w[k, j];
                            w[k, j] = temp;
                        }
                    }
                    // eliminate x[k]
                    for (i = k + 1; i < n; i++)
                    {
                        w[i, k] /= w[k, k];
                        ratio = w[i, k];
                        for (j = k + 1; j < n; j++)
                            w[i, j] -= ratio * w[k, j];
                    }
                }
            }

            if (w[n - 1, n - 1] == 0)
                flag = 0;

            return flag;
        }
    }
}

namespace MatrixInverseComparison
{
    // One implementation is based on https://en.wikipedia.org/wiki/Bareiss_algorithm
    // Another perhaps better implementation is found on the following webpage
    // https://cs.stackexchange.com/questions/124759/determinant-calculation-bareiss-vs-gauss-algorithm
    class BareissAlgorithm
    {
        static public double Determinant1(double[,] M, int n)
        {
            double M00;

            for (int k = 0; k < n; k++)
            {
                if (k - 1 == -1)
                    M00 = 1;
                else
                    M00 = M[k - 1, k - 1];

                for (int i = k + 1; i < n; i++)
                {
                    for (int j = k + 1; j < n; j++)
                    {
                        M[i, j] = (
                            M[i, j] * M[k, k] -
                            M[i, k] * M[k, j]) / M00;
                    }
                }
            }

            return M[n - 1, n - 1];
        }

        static public double Determinant2(double[,] A, int dim)
        {
            if (dim <= 0)
            {
                return 0;
            }

            double sign = 1;

            for (int k = 0; k < dim - 1; k++)
            {
                //Pivot - row swap needed
                if (A[k, k] == 0)
                {
                    int m;
                    for (m = k + 1; m < dim; m++)
                    {
                        if (A[m, k] != 0)
                        {
                            double[] tempRow = new double[dim];

                            for (int i = 0; i < dim; i++)
                                tempRow[i] = A[m, i];
                            for (int i = 0; i < dim; i++)
                                A[m, i] = A[k, i];
                            for (int i = 0; i < dim; i++)
                                A[k, i] = tempRow[i];
                            sign = -sign;
                            break;
                        }
                    }
                    //No entries != 0 found in column k -> det = 0
                    if (m == dim)
                    {
                        return 0;
                    }
                }
                //Apply formula
                for (int i = k + 1; i < dim; i++)
                {
                    for (int j = k + 1; j < dim; j++)
                    {
                        A[i, j] = A[k, k] * A[i, j] - A[i, k] * A[k, j];

                        if (k >= 1)
                        {
                            A[i, j] /= A[k - 1, k - 1];
                        }
                    }
                }
            }

            return sign * A[dim - 1, dim - 1];
        }
    }
}

Blast from the Past 1997 Image of a Matrix by James Pate Williams, Jr.

/*
  Author:  Pate Williams (c) 1997

  "Algorithm 2.3.2 (Image of a Matrix). Given an
  m by n matrix M = (m[i][i]) with 1 <= i <= m and
  1 <= j <= n having coefficients in the field K,
  this algorithm outputs a basis of the image of
  M, i. e. vector space spanned by the columns of
  M." -Henri Cohen- See "A Course in Computational
  Algebraic Number Theory" by Henri Cohen pages
  58-59. We specialize the code to the field Zp.
*/

#include <stdio.h>
#include <stdlib.h>

long** create_matrix(long m, long n)
{
    long i, ** matrix = (long**)calloc(m, sizeof(long*));

    if (!matrix) {
        fprintf(stderr, "fatal error\ninsufficient memory\n");
        fprintf(stderr, "from create_matrix\n");
        exit(1);
    }
    for (i = 0; i < m; i++) {
        matrix[i] = (long*)calloc(n, sizeof(long));
        if (!matrix[i]) {
            fprintf(stderr, "fatal error\ninsufficient memory\n");
            fprintf(stderr, "from create_matrix\n");
            exit(1);
        }
    }
    return matrix;
}

void delete_matrix(long m, long** matrix)
{
    long i;

    for (i = 0; i < m; i++) free(matrix[i]);
    free(matrix);
}

void Euclid_extended(long a, long b, long* u,
    long* v, long* d)
{
    long q, t1, t3, v1, v3;

    *u = 1, * d = a;
    if (b == 0) {
        *v = 0;
        return;
    }
    v1 = 0, v3 = b;
#ifdef DEBUG
    printf("----------------------------------\n");
    printf(" q    t3   *u   *d   t1   v1   v3\n");
    printf("----------------------------------\n");
#endif
    while (v3 != 0) {
        q = *d / v3;
        t3 = *d - q * v3;
        t1 = *u - q * v1;
        *u = v1, * d = v3;
#ifdef DEBUG
        printf("%4ld %4ld %4ld ", q, t3, *u);
        printf("%4ld %4ld %4ld %4ld\n", *d, t1, v1, v3);
#endif
        v1 = t1, v3 = t3;
    }
    *v = (*d - a * *u) / b;
#ifdef DEBUG
    printf("----------------------------------\n");
#endif
}

long inv(long number, long modulus)
{
    long d, u, v;

    Euclid_extended(number, modulus, &u, &v, &d);
    if (d == 1) return u;
    return 0;
}

void image(long m, long n, long p,
    long** M, long** X, long* r)
{
    int found;
    long D, i, j, k, s;
    long* c = (long*)calloc(m, sizeof(long));
    long* d = (long*)calloc(n, sizeof(long));
    long** N = create_matrix(m, n);

    if (!c || !d) {
        fprintf(stderr, "fatal error\ninsufficient memory\n");
        fprintf(stderr, "from kernel\n");
        exit(1);
    }
    for (i = 0; i < m; i++) {
        c[i] = -1;
        for (j = 0; j < n; j++) N[i][j] = M[i][j];
    }
    *r = 0;
    for (k = 0; k < n; k++) {
        found = 0, j = 0;
        while (!found && j < m) {
            found = M[j][k] != 0 && c[j] == -1;
            if (!found) j++;
        }
        if (found) {
            D = p - inv(M[j][k], p);
            M[j][k] = p - 1;
            for (s = k + 1; s < n; s++)
                M[j][s] = (D * M[j][s]) % p;
            for (i = 0; i < m; i++) {
                if (i != j) {
                    D = M[i][k];
                    M[i][k] = 0;
                    for (s = k + 1; s < n; s++) {
                        M[i][s] = (M[i][s] + D * M[j][s]) % p;
                        if (M[i][s] < 0) M[i][s] += p;
                    }
                }
            }
            c[j] = k;
            d[k] = j;
        }
        else {
            *r = *r + 1;
            d[k] = -1;
        }
    }
    for (j = 0; j < m; j++) {
        if (c[j] != -1) {
            for (i = 0; i < n; i++) {
                if (i < m) X[i][j] = N[i][c[j]];
                else X[i][j] = 0;
            }
        }
    }
    delete_matrix(m, N);
    free(c);
    free(d);
}

void print_matrix(long m, long n, long** a)
{
    long i, j;

    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++)
            printf("%2ld ", a[i][j]);
        printf("\n");
    }
}

int main(void)
{
    long i, j, m = 8, n = 8, p = 13, r;
    long a[8][8] = { {0,  0,  0,  0,  0,  0,  0,  0},
                    {2,  0,  7, 11, 10, 12,  5, 11},
                    {3,  6,  3,  3,  0,  4,  7,  2},
                    {4,  3,  6,  4,  1,  6,  2,  3},
                    {2, 11,  8,  8,  2,  1,  3, 11},
                    {6, 11,  8,  6,  2,  6, 10,  9},
                    {5, 11,  7, 10,  0, 11,  6, 12},
                    {3,  3, 12,  5,  0, 11,  9, 11} };
    long** M = create_matrix(m, n);
    long** X = create_matrix(n, n);

    for (i = 0; i < m; i++)
        for (j = 0; j < n; j++)
            M[i][j] = a[j][i];
    printf("the original matrix is as follows:\n");
    print_matrix(m, n, M);
    image(m, n, p, M, X, &r);
    printf("the image of the matrix is as follows:\n");
    print_matrix(n, n - r, X);
    printf("the rank of the matrix is: %ld\n", n - r);
    delete_matrix(m, M);
    delete_matrix(n, X);
    return 0;
}

First Blast from the Past (1997) Computing the Inverse Image of a Matrix by James Pate Williams, Jr.

/*
  Author:  Pate Williams (c) 1997

  "Algorithm 2.3.5 (Inverse Image Matrix). Let M be
  an m by n matrix and V be an m by r matrix, where
  n <= m. This algorithm either outputs a message
  saying that some column vector of V is not in the
  image of M, or outputs an n by r matrix X such
  that V = MX." -Henri Cohen- See "A Course in Com-
  putational Algebraic Number Theory" by Henri
  Cohen pages 60-61. We specialize to the field Q.
*/

#include <stdio.h>
#include <stdlib.h>

double** create_matrix(long m, long n)
{
    double** matrix = (double**)calloc(m, sizeof(double*));
    long i;

    if (!matrix) {
        fprintf(stderr, "fatal error\ninsufficient memory\n");
        fprintf(stderr, "from create_matrix\n");
        exit(1);
    }
    for (i = 0; i < m; i++) {
        matrix[i] = (double*)calloc(n, sizeof(double));
        if (!matrix[i]) {
            fprintf(stderr, "fatal error\ninsufficient memory\n");
            fprintf(stderr, "from create_matrix\n");
            exit(1);
        }
    }
    return matrix;
}

void delete_matrix(long m, double** matrix)
{
    long i;

    for (i = 0; i < m; i++) free(matrix[i]);
    free(matrix);
}

void inverse_image_matrix(long m, long n, long r,
    double** M, double** V,
    double** X)
{
    double ck, d, sum, t;
    double** B = create_matrix(m, r);
    int found;
    long i, j, k, l;

    for (i = 0; i < m; i++)
        for (j = 0; j < r; j++)
            B[i][j] = V[i][j];
    for (j = 0; j < n; j++) {
        found = 0, i = j;
        while (!found && i < m) {
            found = M[i][j] != 0;
            if (!found) i++;
        }
        if (!found) {
            fprintf(stderr, "fatal error\nnot linearly independent\n");
            fprintf(stderr, "from inverse_image_matrix\n");
            exit(1);
        }
        if (i > j) {
            for (l = 0; l < n; l++)
                t = M[i][l], M[i][l] = M[j][l], M[j][l] = t;
            for (l = 0; l < r; l++)
                t = B[i][l], B[i][l] = B[j][l], B[j][l] = t;
        }
        d = 1.0 / M[j][j];
        for (k = j + 1; k < m; k++) {
            ck = d * M[k][j];
            for (l = j + 1; l < n; l++)
                M[k][l] -= ck * M[j][l];
            for (l = 0; l < r; l++)
                B[k][l] -= ck * B[j][l];
        }
    }
    for (i = n - 1; i >= 0; i--) {
        for (k = 0; k < r; k++) {
            sum = 0;
            for (j = i + 1; j < n; j++)
                sum += M[i][j] * X[j][k];
            X[i][k] = (B[i][k] - sum) / M[i][i];
        }
    }
    for (k = n + 1; k < m; k++) {
        for (j = 0; j < r; j++) {
            sum = 0;
            for (i = 0; i < n; i++)
                sum += M[k][i] * X[i][j];
            if (sum != B[k][j]) {
                fprintf(stderr, "fatal error\ncolumn not in image\n");
                fprintf(stderr, "from inverse_image_matrix\n");
                exit(1);
            }
        }
    }
    delete_matrix(m, B);
}

void matrix_multiply(long m, long n, long r,
    double** a, double** b,
    double** c)
    /* c = a * b */
{
    double sum;
    long i, j, k;

    for (i = 0; i < m; i++) {
        for (j = 0; j < r; j++) {
            sum = 0.0;
            for (k = 0; k < n; k++)
                sum += a[i][k] * b[k][j];
            c[i][j] = sum;
        }
    }
}

void print_matrix(long m, long n, double** a)
{
    long i, j;

    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++)
            printf("%+10.6lf ", a[i][j]);
        printf("\n");
    }
}

int main(void)
{
    long i, j, m = 4, n = 4, r = 4;
    double** c = create_matrix(m, n);
    double** M = create_matrix(m, n);
    double** V = create_matrix(m, r);
    double** X = create_matrix(n, r);

    for (i = 0; i < m; i++) {
        c[i][i] = M[i][i] = 2.0;
        if (i > 0)
            c[i][i - 1] = M[i][i - 1] = -1.0;
        if (i < m - 1)
            c[i][i + 1] = M[i][i + 1] = -1.0;
    }
    for (i = 0; i < m; i++)
        for (j = 0; j < r; j++)
            V[i][j] = i + j + 1;
    printf("M is as follows:\n");
    print_matrix(m, n, M);
    printf("V is as follows:\n");
    print_matrix(m, r, V);
    inverse_image_matrix(m, n, r, M, V, X);
    printf("X is as follows:\n");
    print_matrix(n, r, X);
    matrix_multiply(m, n, r, c, X, M);
    printf("MX is as follows:\n");
    print_matrix(m, r, M);
    delete_matrix(m, c);
    delete_matrix(m, M);
    delete_matrix(m, V);
    delete_matrix(n, X);
    return 0;
}

C++ Linear Algebra Package Extension Implemented by James Pate Williams, Jr.

Input file:

The dimensions of the linear system of equations (m and n, m = n):
2
2
The matrix of the linear system of equations (n by n):
1	1
1	2
The right-hand side of the linear system of equations (n by 1):
7	11
The dimensions of the linear system of equations (m and n, m = 2):
2
2
The matrix of the linear system of equations (n by n):
1	1
1	3
The right-hand side of the linear system of equations (n by 1):
7	11
The dimensions of the linear system of equations (m and n, m = n):
2
2
The matrix of the linear system of equations (n by n):
6	3
4	8
The right-hand side of the linear system of equations (n by 1):
5	6
The dimensions of the linear system of equations (m and n, m = n):
2
2
The matrix of the linear system of equations (n by n):
5	3
10	4
The right-hand side of the linear system of equations (n by 1):
8	6
The dimensions of the linear system of equations (m and n, m = n):
3
3
The matrix of the linear system of equations (n by n):
2	1	-1
-3	-1	2
-2	1	2
The right-hand side of the linear system of equations (n by 1):
8	-11	-3

Output file:

The 1st solution of the linear system of equations:
       3	       4	
The 2nd solution of the linear system of equations:
       3	       4	
The determinant of the linear system of equations:
1
The inverse of the linear system of equations:
       2	      -1	
      -1	       1	
The adjoint of the linear system of equations:
       2	      -0	
      -0	       2	
The characteristic polynomial of the linear system of equations:
       1	       2	
The image of the matrix: 
       1	      -1	
       1	       2	
Rank = 2
The 1st solution of the linear system of equations:
       5	       2	
The 2nd solution of the linear system of equations:
       5	       2	
The determinant of the linear system of equations:
2
The inverse of the linear system of equations:
     1.5	    -0.5	
    -0.5	     0.5	
The adjoint of the linear system of equations:
       3	      -0	
      -0	       3	
The characteristic polynomial of the linear system of equations:
       1	       3	
The image of the matrix: 
       1	      -1	
       1	       3	
Rank = 2
The 1st solution of the linear system of equations:
 0.61111	 0.44444	
The 2nd solution of the linear system of equations:
 0.61111	 0.44444	
The determinant of the linear system of equations:
36
The inverse of the linear system of equations:
 0.22222	-0.11111	
-0.083333	 0.16667	
The adjoint of the linear system of equations:
      48	      -0	
      -0	      48	
The characteristic polynomial of the linear system of equations:
       1	      48	
The image of the matrix: 
       6	      -1	
       4	       8	
Rank = 2
The 1st solution of the linear system of equations:
    -1.4	       5	
The 2nd solution of the linear system of equations:
    -1.4	       5	
The determinant of the linear system of equations:
-10
The inverse of the linear system of equations:
    -0.4	       1	
     0.3	    -0.5	
The adjoint of the linear system of equations:
      20	      -0	
      -0	      20	
The characteristic polynomial of the linear system of equations:
       1	      20	
The image of the matrix: 
       5	      -1	
      10	       4	
Rank = 2
The 1st solution of the linear system of equations:
       2	       3	      -1	
The 2nd solution of the linear system of equations:
       2	       3	      -1	
The determinant of the linear system of equations:
1
The inverse of the linear system of equations:
       4	       5	      -2	
       3	       4	      -2	
      -1	      -1	       1	
The adjoint of the linear system of equations:
      -4	       0	       0	
       0	      -4	       0	
       0	       0	      -4	
The characteristic polynomial of the linear system of equations:
       1	      -7	       4	
The image of the matrix: 
       2	      -1	       2	
      -3	       0	      -1	
      -2	       0	       4	
Rank = 3

// 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();
}

Part of a New Linear Algebra Package in C++ Implemented by James Pate Williams, Jr.

// 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];

		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(size_t n, bool failure);
	Vector<double> DblGaussianElimination(
		bool& failure);
	// The following three 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(
		size_t n,
		Vector<int>& pivot);
	bool DblGaussianSolution(
		int n,
		Vector<double>& x,
		Vector<int>& pivot);
	bool DblSubstitution(
		size_t n,
		Vector<double>& x,
		Vector<int>& pivot);
	bool DblInverse(
		size_t n,
		Matrix<double>& A,
		Vector<int>& pivot);
};

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

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

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

    if (!DblGaussianFactor(n, 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(
    size_t n,
    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(
    int n,
    Vector<double>& x,
    Vector<int>& pivot)
{
    if (!DblGaussianFactor(n, pivot))
        return false;

    return DblSubstitution(n, x, pivot);
}

bool LinearAlgebra::DblSubstitution(
    size_t n, 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(
    size_t n,
    Matrix<double>& A,
    Vector<int>& pivot)
{
    Vector<double> x(n);

    if (!DblGaussianFactor(n, 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(n, 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;
}

/*
** 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())
    {
        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 rows must be >= 1" << endl;
            return -101;
        }

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

        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;
            }
        }

        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)
            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(n, X, pivot))
            exit(-103);

        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(n, failure);

        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];
            }
        }

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

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

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

2
2
1	1
1	2
7	11
2
2
1	1
1	3
7	11
2
2
6	3
4	8
5	6
2
2
5	3
10	4
8	6
3
3
2	1	-1
-3	-1	2
-2	1	2
8	-11	-3

       3	       4	
       3	       4	
1
       2	      -1	
      -1	       1	
       5	       2	
       5	       2	
2
     1.5	    -0.5	
    -0.5	     0.5	
 0.61111	 0.44444	
 0.61111	 0.44444	
36
 0.22222	-0.11111	
-0.083333	 0.16667	
    -1.4	       5	
    -1.4	       5	
-10
    -0.4	       1	
     0.3	    -0.5	
       2	       3	      -1	
       2	       3	      -1	
1
       4	       5	      -2	
       3	       4	      -2	
      -1	      -1	       1