Author: jamespatewilliamsjr

My whole legal name is James Pate Williams, Jr. I was born in LaGrange, Georgia approximately 70 years ago. I barely graduated from LaGrange High School with low marks in June 1971. Later in June 1979, I graduated from LaGrange College with a Bachelor of Arts in Chemistry with a little over a 3 out 4 Grade Point Average (GPA). In the Spring Quarter of 1978, I taught myself how to program a Texas Instruments desktop programmable calculator and in the Summer Quarter of 1978 I taught myself Dayton BASIC (Beginner's All-purpose Symbolic Instruction Code) on LaGrange College's Data General Eclipse minicomputer. I took courses in BASIC in the Fall Quarter of 1978 and FORTRAN IV (Formula Translator IV) in the Winter Quarter of 1979. Professor Kenneth Cooper, a genius poly-scientist taught me a course in the Intel 8085 microprocessor architecture and assembly and machine language. We would hand assemble our programs and insert the resulting machine code into our crude wooden box computer which was designed and built by Professor Cooper. From 1990 to 1994 I earned a Bachelor of Science in Computer Science from LaGrange College. I had a 4 out of 4 GPA in the period 1990 to 1994. I took courses in C, COBOL, and Pascal during my BS work. After graduating from LaGrange College a second time in May 1994, I taught myself C++. In December 1995, I started using the Internet and taught myself client-server programming. I created a website in 1997 which had C and C# implementations of algorithms from the "Handbook of Applied Cryptography" by Alfred J. Menezes, et. al., and some other cryptography and number theory textbooks and treatises.

Take the Long Road Home by James Pate Williams, Jr. Excerpts from a by Gone Journal

03-04-2009 Math Grades at Georgia Tech (1980 – 1983)

The applied numerical analysis that I have been doing lately reminds of my not-so-great experience at Georgia Tech in 1980 to 1983. In that long bygone era, I was an untreated schizophrenic later diagnosed as bipolar undergoing some pretty difficult times. I was having some strangely seriously disabling delusions about the CIA and other branches of the federal government. As I have probably stated in other slide presentations on this site, perhaps even too many times, I knew I was delusional, but I could not stop dwelling on these fictitious thoughts. Anyway, I miraculously was able to take and complete a few senior level engineering or math major math courses and three graduate level math courses. Here are the courses and my grades:

MATH 4582 – Advanced Engineering Mathematics – B, 3 Quarter Hours

MATH 4347 – Partial Differential Equations – B, 3 Hours

MATH 4583 – Vector Analysis – B, 3 Hours

MATH 4320 – Complex Analysis – B, 3 Hours

MATH 4338 – Partial Differential Equations – A, 3 Hours

MATH 6581 – Calculus of Variations – A, 3 Hours

MATH 4640 – Scientific Computing I – B, 3 Hours

MATH 4311 – Introduction to Analysis I – C, 4 Hours

MATH 6341 – Partial Differential Equations I – B, 3 Hours

MATH 4312 – Introduction to Analysis II – B, 4 Hours

MATH 6342 – Partial Differential Equations – D, 3 Hours

I have perhaps legitimate excuses for the C and D. I failed to show up for one of my analysis tests and I received a 0. In the PDE course in which I made a D, I failed to do a term computer project, because I temporarily lost access to LaGrange College’s DG Eclipse minicomputer, and I was scared to death of the CDC Cyber mainframe on the GT campus. I had a CDC FORTRAN manual, but I had no training with the CDC Compass OS. Also, I don’t remember if the DG Eclipse in the chemistry x-ray crystallography lab had either a BASIC interpreter or Fortran or Pascal compilers.

A-grades at GT in those days were very hard to get. I can’t remember whether my introduction to Bessel functions and Legendre polynomials was in MATH 4582 or MATH 4338. I had not taken a formal ordinary differential equations course at LaGrange College, but I was an autodidact (self-taught person) in that area of advanced enginneering mathematics.

I think I could take an applied numerical analysis or partial differential equation course at any institution of higher learning today and do fairly well.

03-04-2009 CHEM 6151 1982 Georgia Tech

I took CHEM 6151 X-Ray Crystallography in the winter of 1982 under the Head of the Chemistry Department. The professor was a top-notched x-ray crystallographer and researcher. Anyway, I wound up with a B in the course since I did not do the optional actual x-ray crystallography analysis of a previously unknown crystal structure. A true chemical, physical, and mathematical genius who had 4.0 averages at Tech both as an undergraduate and graduate student did very well on his crystal and, of course, made an A in the course. This guy was studying nuclear chemistry, and he could easily unravel and decipher any type of spectral data whether it was done by alpha, beta, or gamma spectroscopy, IR, mass, NMR, UV, or visible spectroscopy. I could with great difficulty analyze some spectral data, but I was not even marginally as proficient as the genius.

Again, I think I was afraid of the x-ray crystallography equipment since I thought the government could use spurious x-rays to read my thoughts. I could have used an aluminum hat to deflect the x-rays, but I was even too paranoid for that insane gesture.

03-06-2009 Bailout Insanity

I voted for President Obama, and I really don’t think I had any other viable option with the exception of a vote for the unelectable Ralph Nader. I do wonder about the sanity of the financial sector bailout initiated by former President Bush (good riddance) to the tune of $350 billion out of $700 billion available and also the most recent economic stimulus package signed by President Obama. All these “bailouts or stimuli” are very reminiscent of the desperate measures taken by FDR during the onset of the great depression. Don’t get me wrong, I like FDR a lot, but the massive public works programs he started were nothing more than a somewhat ineffective economic band-aid. It took the Imperialistic Japanese and the mentally ill Nazis with WWII to really recover to a certain extent our economy. Then there was the great economic bust immediately following WWII when the factories geared down production from the peak WWII levels, released all the highly skilled female employees, and the large number of armed services personnel returned to no jobs. Only the GI bill, Marshall Plan, and other ingenious economic downturn countermeasures saved us for the baby booming times of the late 1940s and early 1950s. You can’t just wantonly throw federal money at a set of serious economic problems without some very long term strategic as well as tactical planning. Just creating a lot of construction jobs is a 1930s approach and inherently destabilizing today. As the old actor (“Water World” villain whose name evades me) states in a television commercial about retirement: “You have to have a plan”. Well, I don’t think either the former President, current President, or Congress really has long range plan and what we are looking at now is a lot of panic, put your finger in the dike thinking.

P. S. Fortunately, for me, President Bush enacted the second tier of unemployment benefits that I am collecting for about the next 10 or so weeks.

P. S. S. I thought the pundits said that the bottom of the stock market was at the INDU level of 7,000 now it is about 6,500. AFLAC stock has lost almost 5/6’s of its peak value about last March. Google has lost over half its stock value. We have similar situations with Amazon and Apple stocks. The technology stocks are still better off than the financial related stocks with Bank of America stock worth less than a $1.00 a share. In LaGrange, Georgia, U. S., the January unemployment rate was over 14%, a hefty increase from 6.7% a year ago. My county, Troup County’s unemployment rate was 12.2% up from 10.3% in December 2008 and 6.7% in December 2007. These are the scariest economic times that I can remember. I am ready to be a greeter at WAL-MART, if and only if, I can get the job.

Multiple Precision Arithmetic in C++ Implemented by James Pate Williams, Jr.

I have developed a long integer package in C++ using H. T. Lau’s “A Numerical Library in C for Scientists and Engineers”. That library is based on the NUMAL Numerical Algorithms in Algol library. I use 32-bit integers (long) as the basis of the LongInteger type. The base (radix) is 10000 which is the largest base using 32-bit integers. As a test of the library, I use Pollard’s original rho factorization method. I utilized the Alfred J. Menezes et al “Handbook of Applied Cryptography” Miller-Rabin algorithm and ANSI X9.17 pseudo random number generator with Triple-DES as the encryption algorithm. I translated Hans Riesel Pascal code for Euclid’s algorithm and the power modulus technique. I don’t use dynamic long integers a la Hans Riesel’s Pascal multiple precision library. The single precision is 32 32-bit longs and multiple precision 64 32-bit longs.

Here is a typical factorization:

Number to be factored, N = 3 1234 5678 9012
Factor = 3
Is Factor prime? 1
Factor = 4
Is Factor prime? 0
Factor = 102 8806 5751
Is Factor prime? 1
Function Evaluations = 6
Number to be factored, N = 0

The first number of N is the number of base 10000 digits. I verified that 10288065751 was prime using Miller-Rabin and the table found online below:

http://compoasso.free.fr/primelistweb/page/prime/liste_online_en.php

Solving the Low-Density Subset Sum Problem with the LLL Lattice Reduction Algorithm C# Implementation by James Pate Williams, Jr.

A few years ago, I implemented the LLL lattice reduction algorithm found in two references: “Handbook of Applied Cryptography” Edited by Alfred J. Menezes et al and “A Course in Computational Algebraic Number Theory” by Henri Cohen. My new implementation in C# uses 64-bit integers and BigIntegers. Here are some results.

Comparison of Brent’s Modification of the Pollard rho Factoring Method and the Original Pollard rho Method C# Implementations by James Pate Williams, Jr.

27-Decimal Digit Number 123456789012345678901234567

31-Digit Decimal Number 1234567890123456789012345678901 Pollard Failed

The same factorizations by the Classical Shor-Pollard- James Pate Williams, Jr. algorithm.

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

Powering Algorithms from Heri Cohen’s Textbook 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

using System.Collections.Generic;
using System.Numerics;

namespace PoweringAlgorithms
{
    public interface IAlgorithms
    {
        // Extended Euclidean Algorithm
        BigInteger Algorithm_1_2_1(
            BigInteger g, BigInteger n);
        BigInteger Algorithm_1_2_2(
            int e, BigInteger g, BigInteger n);
        public BigInteger Algorithm_1_2_3(
            int e, BigInteger g, BigInteger n);
        void Algorithm_1_3_6(
            BigInteger a, BigInteger b,
            out BigInteger u, out BigInteger v,
            out BigInteger d);
        List<int> GetBits(BigInteger N);
        BigInteger Inverse(
            BigInteger a, BigInteger n);
    }
}

using System.Collections.Generic;
using System.Numerics;

namespace PoweringAlgorithms
{
    public class Algorithms : IAlgorithms
    { 
        public BigInteger Algorithm_1_2_1(
            BigInteger g, BigInteger n)
        {
            BigInteger y = 1;

            if (n == 0)
                return y;

            BigInteger N, z;

            if (n < 0)
            {
                N = -n;
                z = Inverse(g, n);
            }

            else
            {
                N = n;
                z = g;
            }

            while (N > 0)
            {
                if (N % 2 == 1)
                    y *= z;

                N /= 2;
                z *= z;
            }

            return y;
        }

        public BigInteger Algorithm_1_2_2(
            int e, BigInteger g, BigInteger n)
        {
            BigInteger E, N, y, z;

            if (n == 0)
                return BigInteger.One;

            if (n < 0)
            {
                N = -n;
                z = Inverse(g, n);
            }

            else
            {
                N = n;
                z = g;
            }
            
            y = z;
            E = BigInteger.Pow(2, e);
            N -= E;

            while (E > 1)
            {
                E /= 2;
                y *= y;

                if (N >= E)
                {
                    N -= E;
                    y *= z;
                }
            }

            return y;
        }

        public BigInteger Algorithm_1_2_3(
            int e, BigInteger g, BigInteger n)
        {
            BigInteger y = 1;

            if (n == 0)
                return y;

            BigInteger N, z;

            if (n < 0)
            {
                N = -n;
                z = Inverse(g, n);
            }

            else
            {
                N = n;
                z = g;
            }

            y = z;

            int f = e;
            List<int> bits = GetBits(N);

            while (f > 0)
            {
                f--;
                y *= y;
                if (bits[f] == 1)
                    y *= z;
            }

            return y;
        }

        public void Algorithm_1_3_6(
                    BigInteger a, BigInteger b,
                    out BigInteger u, out BigInteger v,
                    out BigInteger d)
        {
            u = 1;
            d = a;

            if (b == 0)
            {
                v = 0;
                return;
            }

            BigInteger q, t1, t3, v1 = 0, v3 = b;

            while (v3 > 0)
            {
                q = d / v3;
                t3 = d % v3;
                t1 = u - q * v1;
                u = v1;
                d = v3;
                v1 = t1;
                v3 = t3;
            }

            v = (d - a * u) / b;
        }

        public List<int> GetBits(BigInteger N)
        {
            int b;
            List<int> temp = new();

            while (N > 0)
            {
                b = (int)(N % 2);
                N /= 2;
                temp.Add(b);
            }

            List<int> bits = new();

            for (int i = temp.Count; i >= 0; i--)
                bits.Add(temp[i]);

            return bits;
        }

        public BigInteger Inverse(BigInteger a, BigInteger n)
        {
            BigInteger d, u, v;

            Algorithm_1_3_6(
                a, n, out u, out v, out d);
            
            return u;
        }
    }
}

using System;
using System.Numerics;
using System.Windows.Forms;

namespace PoweringAlgorithms
{
    public partial class TestForm : Form
    {
        public TestForm()
        {
            InitializeComponent();
        }

        private static BigInteger Horner(string text)
        {
            BigInteger a = text[0] - '0';

            for (int i = 1; i < text.Length; i++)
                a = 10 * a + text[i] - '0';

            return a;
        }

        private void button1_Click(object sender, EventArgs e)
        {
            int expon = 0;
            Algorithms algos = new Algorithms();
            BigInteger g = Horner(textBox1.Text);
            BigInteger n = Horner(textBox2.Text);
            BigInteger a = algos.Algorithm_1_2_1(g, n);
            BigInteger m = BigInteger.Abs(n);

            while (m > 0)
            {
                expon++;
                m /= 2;
            }
            
            BigInteger b = algos.Algorithm_1_2_2(expon - 1, g, n);
            BigInteger c = algos.Algorithm_1_2_2(expon - 1, g, n);

            textBox3.Text = a.ToString();
            textBox4.Text = b.ToString();
            textBox5.Text = c.ToString();
        }

        private void button2_Click(object sender, EventArgs e)
        {
            Close();
        }
    }
}

Huffman Compression in C# Implemented by James Pate Williams, Jr.

Algorithms are found in the textbook “Introduction to Algorithms” by Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest p. 340.

using System;
using System.Collections.Generic;
using System.Windows.Forms;

namespace Huffman
{
    public partial class MainForm : Form
    {
        private int leafNodes;

        public MainForm()
        {
            InitializeComponent();
        }

        private void InorderTraversal(BinaryTreeNode<CharFreq> node)
        {
            if (node != null)
            {
                InorderTraversal(node.Left);

                CharFreq cf = node.Value;
                int ord = (int)cf.ch;

                if (node.Left == null && node.Right == null)
                {
                    textBox2.Text += leafNodes.ToString("F0").PadLeft(3) + '\t';
                    textBox2.Text += "'" + new string(cf.ch, 1) + "' " + '\t';
                    textBox2.Text += node.Value.freq.ToString() + "\r\n";
                    leafNodes++;
                }

                InorderTraversal(node.Right);
            }
        }

        private void button1_Click(object sender, EventArgs e)
        {
            string s = textBox1.Text;
            int n = s.Length;
            List<CharFreq> list = new List<CharFreq>();

            textBox2.Text = string.Empty;

            for (int i = 0; i < n; i++)
            {
                bool found = false;
                char c = s[i];
                CharFreq cf = new CharFreq();

                for (int j = 0; !found && j < list.Count; j++)
                {
                    if (c == list[j].ch)
                    {
                        found = true;
                        cf.ch = c;
                        cf.freq = 1 + list[j].freq;
                        list.RemoveAt(j);
                        list.Add(cf);
                    }
                }

                if (!found)
                {
                    cf.ch = c;
                    cf.freq = 1;
                    list.Add(cf);
                }
            }
                        
            HuffmanTree ht = new HuffmanTree();
            BinaryTreeNode<CharFreq> root = ht.Build(list, list.Count);

            InorderTraversal(root);
            textBox2.Text += "\r\n# characters = " + n.ToString() + "\r\n";
            textBox2.Text += "# leaf nodes = " + leafNodes.ToString() + "\r\n";
            textBox2.Text += "% compressed = " + 
                (100.0 - 100.0 * ((double)leafNodes) / n).ToString("F2") + "\r\n"; 
        }
    }
}

namespace Huffman
{
    public class BinaryTreeNode<T> : Node<T>
    {
        public BinaryTreeNode() : base() { }
        public BinaryTreeNode(T data) : base(data, null) { }
        public BinaryTreeNode(T data, BinaryTreeNode<T> left, BinaryTreeNode<T> right)
        {
            base.Value = data;
            NodeList<T> children = new NodeList<T>(2)
            {
                [0] = left,
                [1] = right
            };

            base.Neighbors = children;
        }

        public BinaryTreeNode<T> Left
        {
            get
            {
                if (base.Neighbors == null)
                    return null;
                else
                    return (BinaryTreeNode<T>)base.Neighbors[0];
            }
            set
            {
                if (base.Neighbors == null)
                    base.Neighbors = new NodeList<T>(2);

                base.Neighbors[0] = value;
            }
        }

        public BinaryTreeNode<T> Right
        {
            get
            {
                if (base.Neighbors == null)
                    return null;
                else
                    return (BinaryTreeNode<T>)base.Neighbors[1];
            }
            set
            {
                if (base.Neighbors == null)
                    base.Neighbors = new NodeList<T>(2);

                base.Neighbors[1] = value;
            }
        }
    }
}

using System.Collections.Generic;

namespace Huffman
{
    public class HuffmanTree
    {
        public BinaryTreeNode<CharFreq> Build(List<CharFreq> charFreq, int n)
        {
            PriorityQueue Q = new PriorityQueue();

            for (int i = 0; i < n; i++)
            {
                BinaryTreeNode<CharFreq> z = new BinaryTreeNode<CharFreq>(charFreq[i]);

                Q.insert(z);
            }

            Q.buildHeap();

            for (int i = 0; i < n - 1; i++)
            {
                BinaryTreeNode<CharFreq> x = Q.extractMin();
                BinaryTreeNode<CharFreq> y = Q.extractMin();
                CharFreq chFreq = new CharFreq();

                chFreq.ch = (char)((int)x.Value.ch + (int)y.Value.ch);
                chFreq.freq = x.Value.freq + y.Value.freq;

                BinaryTreeNode<CharFreq> z = new BinaryTreeNode<CharFreq>(chFreq);

                z.Left = x;
                z.Right = y;
                Q.insert(z);
            }

            return Q.extractMin();
        }
    }
}

namespace Huffman
{
    public class Node<T>
    {
        // Private member-variables
        private T data;
        private NodeList<T> neighbors = null;

        public Node() { }
        public Node(T data) : this(data, null) { }
        public Node(T data, NodeList<T> neighbors)
        {
            this.data = data;
            this.neighbors = neighbors;
        }

        public T Value
        {
            get
            {
                return data;
            }
            set
            {
                data = value;
            }
        }

        protected NodeList<T> Neighbors
        {
            get
            {
                return neighbors;
            }
            set
            {
                neighbors = value;
            }
        }
    }
}

using System.Collections.ObjectModel;

namespace Huffman
{
    public class NodeList<T> : Collection<Node<T>>
    {
        public NodeList() : base() { }

        public NodeList(int initialSize)
        {
            // Add the specified number of items
            for (int i = 0; i < initialSize; i++)
                base.Items.Add(default(Node<T>));
        }

        public Node<T> FindByValue(T value)
        {
            // search the list for the value
            foreach (Node<T> node in Items)
                if (node.Value.Equals(value))
                    return node;

            // if we reached here, we didn't find a matching node
            return null;
        }
    }
}

using System.Collections.Generic;

namespace Huffman
{
    public struct CharFreq
    {
        public char ch;
        public int freq;
    }

    public class PriorityQueue
    {
        int heapSize;
        List<BinaryTreeNode<CharFreq>> nodeList;

        public List<BinaryTreeNode<CharFreq>> NodeList
        {
            get
            {
                return nodeList;
            }
        }

        public PriorityQueue()
        {
            nodeList = new List<BinaryTreeNode<CharFreq>>();
        }

        public PriorityQueue(List<BinaryTreeNode<CharFreq>> nl)
        {
            heapSize = nl.Count;
            nodeList = new List<BinaryTreeNode<CharFreq>>();

            for (int i = 0; i < nl.Count; i++)
                nodeList.Add(nl[i]);
        }

        public void exchange(int i, int j)
        {
            BinaryTreeNode<CharFreq> temp = nodeList[i];

            nodeList[i] = nodeList[j];
            nodeList[j] = temp;
        }

        public void heapify(int i)
        {
            int l = 2 * i + 1;
            int r = 2 * i + 2;
            int largest = -1;

            if (l < heapSize && nodeList[l].Value.ch > nodeList[i].Value.ch)
                largest = l;
            else
                largest = i;
            if (r < heapSize && nodeList[r].Value.ch > nodeList[largest].Value.ch)
                largest = r;
            if (largest != i)
            {
                exchange(i, largest);
                heapify(largest);
            }
        }

        public void buildHeap()
        {
            for (int i = heapSize / 2; i >= 0; i--)
                heapify(i);
        }

        public int size()
        {
            return heapSize;
        }

        public BinaryTreeNode<CharFreq> elementAt(int i)
        {
            return nodeList[i];
        }

        public void heapSort()
        {
            int temp = heapSize;

            buildHeap();

            for (int i = heapSize - 1; i >= 1; i--)
            {
                exchange(0, i);
                heapSize--;
                heapify(0);
            }

            heapSize = temp;
        }

        public BinaryTreeNode<CharFreq> extractMin()
        {
            if (heapSize < 1)
                return null;

            heapSort();

            exchange(0, heapSize - 1);
            heapSize--;

            BinaryTreeNode<CharFreq> node = nodeList[heapSize];

            nodeList.RemoveAt(heapSize);
            heapSize = nodeList.Count;
            return node;
        }

        public void insert(BinaryTreeNode<CharFreq> node)
        {
            nodeList.Add(node);
            heapSize = nodeList.Count;
            buildHeap();
        }
    }
}

Huffman Compression in C++ Implemented by James Pate Williams, Jr.

The original string is:
abbbccddeefffghhhhijkllmm

# characters = 25
The compressed codes and frequencies are:
  0     e         2
  1     l         2
  2     f         3
  3     b         3
  4     i         1
  5     c         2
  6     g         1
  7     a         1
  8     d         2
  9     m         2
 10     k         1
 11     j         1
 12     h         4
# leaf nodes = 13
% compressed = 48

C:\Users\james\source\repos\HuffmanCodes\Debug\HuffmanCodes.exe (process 36772) exited with code 0.
Press any key to close this window . . .
// Algorithm is found in the textbook
// "Introduction to Algorithms"
// by Thomas H. Cormen, Charles E.
// Leiserson, Ronald L. Rivest p. 340

#include "pch.h"

int leafNodes = 0;

void InorderTraversal(BinaryTreeNode<CharFreq>* node)
{
    if (node != NULL)
    {
        InorderTraversal(node->lt);

        if (node->lt == NULL && node->rt == NULL)
        {
            CharFreq cf = node->data;

            std::cout << setw(3) << leafNodes << '\t';
            std::cout << cf.ch << '\t';
            std::cout << setw(3) << cf.freq << '\n';
            leafNodes++;
        }
        
        InorderTraversal(node->rt);
    }
}

int main()
{
    int f[128] = { 0 };
    string str = "abbbccddeefffghhhhijkllmm";
    BinaryTreeNode<CharFreq> charFreqTree;
    vector<BinaryTreeNode<CharFreq>> v;

    std::cout << "The original string is: " << endl;
    std::cout << str << endl << endl;

    for (size_t i = 0; i < strlen(str.c_str()); i++)
    {
        bool found = false;
        char ch = str.c_str()[i];

        for (auto iter = v.begin(); !found &&
            iter != v.end(); iter++)
        {
            BinaryTreeNode<CharFreq> node = *iter;

            if (node.data.ch == ch)
            {
                node.data.freq++;
                *iter = node;
                found = true;
            }
        }

        if (!found)
        {
            BinaryTreeNode<CharFreq> node;

            node.data.ch = ch;
            node.data.freq = 1;
            node.lt = node.rt = NULL;
            v.push_back(node);
        }
    }

    priority_queue<BinaryTreeNode<CharFreq>, vector<BinaryTreeNode<CharFreq>>,
        greater<BinaryTreeNode<CharFreq>>> Q(v.begin(), v.end());

    size_t n = Q.size();
   
    for (size_t i = 0; i < n - 1; i++)
    {
        BinaryTreeNode<CharFreq>* x = new
            BinaryTreeNode<CharFreq>();
        BinaryTreeNode<CharFreq>* y = new
            BinaryTreeNode<CharFreq>();
        *x = Q.top();
        Q.pop();
        *y = Q.top();
        Q.pop();

        CharFreq charFreq;
        charFreq.ch = (char)(x->data.ch + y->data.ch);
        charFreq.freq = x->data.freq + y->data.freq;

        BinaryTreeNode<CharFreq>* z = new
            BinaryTreeNode<CharFreq>(charFreq, x, y);

        Q.push(*z);
    }

    BinaryTreeNode<CharFreq> root = Q.top();
    std::cout << "# characters = " << strlen(str.c_str()) << endl;
    std::cout << "The compressed codes and frequencies are:" << endl;
    InorderTraversal(&root);
    std::cout << "# leaf nodes = " << leafNodes << endl;
    std::cout << "% compressed = " <<
        (100.0 - 100.0 * ((double)leafNodes) / strlen(str.c_str())) << endl;
    return 0;
}

#pragma once
#include "pch.h"
using namespace std;

template <class T>
	class BinaryTreeNode
	{
	public:
		T data;
		BinaryTreeNode* lt;
		BinaryTreeNode* rt;

		BinaryTreeNode() { 
			lt = rt = nullptr;
		};

		BinaryTreeNode(T data)
		{
			this->data = data;
			lt = rt = nullptr;
		};

		BinaryTreeNode(T data, BinaryTreeNode* lt, BinaryTreeNode* rt)
		{
			this->data = data;
			this->lt = lt;
			this->rt = rt;
		};

		friend bool operator > (BinaryTreeNode lhs, BinaryTreeNode rhs)
		{
			return lhs.data > rhs.data;
		};

		friend bool operator < (BinaryTreeNode lhs, BinaryTreeNode rhs)
		{
			return lhs.data < rhs.data;
		};

		friend bool operator == (BinaryTreeNode lhs, BinaryTreeNode rhs)
		{
			return lhs.data == rhs.data;
		};
	};

#pragma once
class CharFreq
{
public:
	char ch;
	int freq;
	
	CharFreq()
	{
		ch = '\0';
		freq = 0;
	};
	CharFreq(char c)
	{
		ch = c;
		freq = 0;
	};
	CharFreq(char c, int f)
	{
		ch = c;
		freq = f;
	};

	friend int operator - (CharFreq lhs, CharFreq rhs)
	{
		return lhs.freq - rhs.freq;
	}

	friend bool operator > (CharFreq lhs, CharFreq rhs)
	{
		return lhs.freq > rhs.freq;
	};

	friend bool operator < (CharFreq lhs, CharFreq rhs)
	{
		return lhs.freq < rhs.freq;
	};

	friend bool operator == (CharFreq lhs, CharFreq rhs)
	{
		return lhs.freq == rhs.freq && lhs.ch == rhs.ch;
	};
};

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

// add headers that you want to pre-compile here
#include "BinaryTreeNode.h"
#include "CharFreq.h"
#include <iomanip>
#include <iostream>
#include <list>
#include <queue>
#include <string>
using namespace std;
#endif //PCH_H

Infix Notation to Postfix Notation and Postfix Evaluator Implemented by James Pate Williams, Jr.

I translated a Pascal program that is found in “Applied Data Structures Using Pascal” by Guy J. Hale and Richard J. Easton. The original Pascal program only used single digits and four arithmetic operators: ‘*’, ‘/’, ‘+’, and ‘-‘. I extended the code to multiple digit positive integers and added an exponentiation operator ‘^’. The priority of the operators is ‘^’, ‘*’, and ‘/’ highest value and ‘+’ and ‘-‘ lowest priority. I could easily add a modulus operator ‘%’ and Boolean bit operators. Extension to negative integers should be a facile operation. Below is the output from my C++ application.

3 + 7 * 8 - 5
3 7 8 * + 5 - 
Positive integer value = 54
3 * 7 - 4 / 2
3 7 * 4 2 / - 
Positive integer value = 19
(3 + 7) * 8 - 5
3 7 + 8 * 5 - 
Positive integer value = 75
(3 + 4) * 8 - (7 * 3 - 4)
3 4 + 8 * 7 3 * 4 - - 
Positive integer value = 39
(100 + 50) * 20 - 100 / 2
100 50 + 20 * 100 2 / - 
Positive integer value = 2950
2 ^ 16 - 5 * 100
2 16 ^ 5 100 * - 
Positive integer value = 65036


#pragma once
#include <list>
#include <stack>
#include <string>
#include <vector>
using namespace std;

const char Exp = '^';
const char Mul = '*';
const char Div = '/';
const char Add = '+';
const char Sub = '-';
const char LtParen = '(';
const char RtParen = ')';

class InfixToPostFix
{
public:
	stack<char> numberStack;
	
	int Convert(
		string infixStr, string& postfixStr);
	int EvaluatePostFix(string postfixStr);
	int Priority(char opcode);
};

#include "pch.h"
#include "InfixToPostFix.h"
#include <vector>
using namespace std;

int InfixToPostFix::Convert(
	string infixStr, string& postfixStr)
{
	char ch = 0;
	int number = 0, opcode = 0, opcode1 = 0, parenLevel = 0;
	string numberStr;

	for (size_t i = 0; i < infixStr.size();)
	{
		//while (i < infixStr.size() && infixStr[i] == ' ')
			//i++;

		while (i < infixStr.size() && infixStr[i] >= '0' &&
			infixStr[i] <= '9')
			numberStr.push_back(infixStr[i++]);

		if (numberStr.size() != 0)
		{
			for (size_t j = 0; j < numberStr.size(); j++)
				postfixStr.push_back(numberStr[j]);

			postfixStr.push_back(' ');
			numberStr.clear();
		}

		//while (i < infixStr.size() && infixStr[i] == ' ')
			//i++;

		if (infixStr[i] == '(')
		{
			numberStack.push(infixStr[i]);
			parenLevel++;
		}

		//while (i < infixStr.size() && infixStr[i] == ' ')
			//i++;

		if (infixStr[i] == ')')
		{
			ch = numberStack.top();
			numberStack.pop();

			parenLevel--;

			//while (i < infixStr.size() && infixStr[i] == ' ')
				//i++;

			while (ch != '(')
			{
				postfixStr.push_back(ch);
				postfixStr.push_back(' ');

				parenLevel++;

				ch = numberStack.top();
				numberStack.pop();
			}
		}

		//while (i < infixStr.size() && infixStr[i] == ' ')
			//i++;

		if (infixStr[i] == '^' || infixStr[i] == '*' ||
			infixStr[i] == '/' || infixStr[i] == '+' ||
			infixStr[i] == '-')
		{
			if (numberStack.empty())
				numberStack.push(infixStr[i]);
			else
			{
				ch = numberStack.top();
				numberStack.pop();

				//while (i < infixStr.size() && infixStr[i] == ' ')
					//i++;

				while (Priority(ch) >= Priority(infixStr[i]) &&
					!numberStack.empty() && ch != '(')
				{
					postfixStr.push_back(ch);
					postfixStr.push_back(' ');

					ch = numberStack.top();
					numberStack.pop();
				}

				//while (i < infixStr.size() && infixStr[i] == ' ')
					//i++;

				if (ch != '(')
				{
					if (Priority(infixStr[i]) <= Priority(ch))
					{
						postfixStr.push_back(ch);
						postfixStr.push_back(' ');

						numberStack.push(infixStr[i]);
					}

					else
					{
						numberStack.push(ch);
						numberStack.push(infixStr[i]);
					}
				}
				else
				{
					numberStack.push(ch);
					numberStack.push(infixStr[i]);
				}
			}
		}

		i++;
	}

	while (!numberStack.empty())
	{
		ch = numberStack.top();
		numberStack.pop();

		postfixStr.push_back(ch);
		postfixStr.push_back(' ');
	}

	return 0;
}

int InfixToPostFix::EvaluatePostFix(string postfixStr)
{
	char opcode = 0;
	int charValue = 0, result = 0, value1 = 0, value2 = 0;
	stack<int> intStack;

	for (size_t i = 0; i < postfixStr.size();)
	{
		string numberStr;

		if (postfixStr[i] != ' ')
		{
			while (postfixStr[i] >= '0' && postfixStr[i] <= '9')
				numberStr.push_back(postfixStr[i++]);

			if (!numberStr.empty())
				intStack.push(atoi(numberStr.c_str()));

			else
			{
				value2 = intStack.top();
				intStack.pop();
				value1 = intStack.top();
				intStack.pop();

				opcode = postfixStr[i++];

				if (opcode == '^')
					result = (int)pow(value1, value2);
				else if (opcode == '*')
					result = value1 * value2;
				else if (opcode == '/')
					result = value1 / value2;
				else if (opcode == '+')
					result = value1 + value2;
				else if (opcode == '-')
					result = value1 - value2;

				intStack.push(result);
			}
		}

		else
			i++;
	}

	result = intStack.top();
	intStack.pop();

	return result;
}

int InfixToPostFix::Priority(char opcode)
{
	int result = 0;

	switch (opcode)
	{
	case Exp:
		result = 2;
		break;
	case Mul:
		result = 2;
		break;
	case Div:
		result = 2;
		break;
	case Add:
		result = 1;
		break;
	case Sub:
		result = 1;
		break;
	case LtParen:
		result = 0;
	}

	return result;
}

#include "pch.h"
#include "InfixToPostFix.h"
#include <fstream>
#include <iostream>
#include <string>
using namespace std;

int main()
{
    fstream inps, outs;
    char line[256];
    
    inps.open("TestFile.txt", fstream::in);
    outs.open("ResuFile.txt", fstream::out | fstream::ate);

    while (!inps.eof())
    {
        string postfixStr;

        inps.getline(line, 256);

        if (strlen(line) == 0)
            break;

        string infixStr(line, strlen(line));
        InfixToPostFix translate;

        int con = translate.Convert(
            infixStr, postfixStr);

        if (con != 0)
        {
            cout << "Conversion error!" << endl;
            cout << "Error code = " << con << endl;
        }

        else
        {
            char newline[] = { '\n' };

            outs.write(infixStr.c_str(), infixStr.size());
            outs.write(newline, 1);
            outs.write(postfixStr.c_str(), postfixStr.size());
            outs.write(newline, 1);

            int val = translate.EvaluatePostFix(postfixStr);

            if (val < 0)
            {
                cout << "Evaluation error!" << endl;
                cout << "Error code = " << val << endl;
            }

            else
            {
                char buffer[256] = { '\0' };
                string str = "Positive integer value = ";

                _itoa_s(val, buffer, 10);
                
                outs.write(str.c_str(), strlen(str.c_str()));
                outs.write(buffer, strlen(buffer));
                outs.write(newline, 1);
            }
        }
    }

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

    return 0;
}