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.
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();
}
// 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
// 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();
}
}
}
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();
}
}
}
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
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;
}
I translated to C++ and enhanced a Pascal arithmetic expression evaluator program found in “Applied Data Structures Using Pascal” by Guy J. Hale and Richard E. Easton. The original code used single digit numbers. As an exercise I enhanced the application to utilize multiple digit numbers. Below is my test file and output from my program. I used the C++ standard library stack data structure.
20 + 3 * 4 + 50 * 4 * 2 + 6 * 2 - 8 / 2 + 2 ^ 5
#pragma once
#include <fstream>
#include <stack>
using namespace std;
class Expression
{
public:
char ch, sign, termOp;
int number;
stack<int> stk;
void GetChar(fstream& inps);
void GetExpression(fstream& inps);
void GetFactor(fstream& inps);
void GetTerm(fstream& inps);
};
#include "pch.h"
#include "Expression.h"
#include <math.h>
#include <fstream>
#include <stack>
#include <string>
using namespace std;
void Expression::GetChar(
fstream& inps)
{
while (!inps.eof())
{
ch = (char)inps.get();
if (inps.eof())
exit(1);
if (ch >= '0' && ch <= '9')
return;
if (ch == '^' || ch == '*' || ch == '/')
return;
if (ch == '+' || ch == '-')
return;
if (ch == ' ' || ch == '\t' || ch == '\n')
continue;
if (ch == ';')
return;
}
}
void Expression::GetExpression(fstream& inps)
{
int num, num1, num2;
if (ch == '+' || ch == '-')
{
sign = ch;
GetChar(inps);
GetTerm(inps);
if (sign == '-')
{
num = stk.top();
stk.pop();
num = -num;
stk.push(num);
}
}
GetTerm(inps);
while (ch == '+' || ch == '-')
{
termOp = ch;
GetChar(inps);
GetTerm(inps);
num2 = stk.top();
stk.pop();
num1 = stk.top();
stk.pop();
if (termOp == '+')
num = num1 + num2;
else
num = num1 - num2;
stk.push(num);
}
}
void Expression::GetFactor(fstream& inps)
{
if (ch >= '0' && ch <= '9')
{
string str;
while (ch >= '0' && ch <= '9' && !inps.eof())
{
str += ch;
ch = (char)inps.get();
if (ch == ' ' || ch == '\t' || ch == '\n')
break;
else if (ch == '^' || ch == '+' || ch == '-' ||
ch == '*' || ch == '/')
break;
}
if (inps.eof())
exit(-1);
while (ch == ' ' || ch == '\t' || ch == '\n')
ch = (char)inps.get();
number = atoi(str.c_str());
stk.push(number);
return;
}
else
{
GetChar(inps);
GetExpression(inps);
GetChar(inps);
}
}
void Expression::GetTerm(fstream& inps)
{
char factOp;
int num, num1, num2;
GetFactor(inps);
while (ch == '*' || ch == '/' || ch == '^')
{
factOp = ch;
GetChar(inps);
GetFactor(inps);
num2 = stk.top();
stk.pop();
num1 = stk.top();
stk.pop();
if (factOp == '*')
num = num1 * num2;
else if (factOp == '/')
num = num1 / num2;
else if (factOp == '^')
num = (int)pow(num1, num2);
stk.push(num);
}
}
#include "pch.h"
#include "Expression.h"
#include <fstream>
#include <iostream>
using namespace std;
int main()
{
bool validExp;
int expVal;
fstream inps;
Expression ex;
inps.open("TestExp.txt", fstream::in);
if (!inps.eof())
{
validExp = true;
ex.GetChar(inps);
ex.GetExpression(inps);
expVal = ex.stk.top();
ex.stk.pop();
cout << "Value = " << expVal << endl;
}
else
{
validExp = false;
expVal = 0;
}
return 0;
}
I am working my way through two compiler textbooks: “Design of Compilers Techniques of Programming Language Translation” by Karen A. Lemone and “Modern Compiler Implementation in Java” by Andrew W. Appel. My first exercise is a single line by line assignment statement parser.
Here is my source code and my translation structures:
X1:=a+bb*12; X2:=a/2+bb*12; Identifiers: X1 a bb Literals: 12 Operators: := + * Punctuation: ; Identifiers: X2 a bb Literals: 2 12 Operators: := / + * Punctuation: ;
#pragma once
#include "RegularExpAssignStm.h"
#include <string>
#include <vector>
using namespace std;
class RegularExpAssignStm
{
public:
string punctuation[3] = { ";", "(", ")" };
string upperCase =
"ABCDEFGHIJKLMNOPQRSTUVWXYZ";
string lowerCase =
"abcdefghijklmnopqrstuvwxyz";
string dig = "0123456789";
string ops[5] = { "+", "-", "*", "/", ":="};
vector<string> identifier;
vector<string> liter;
vector<string> oper;
vector<string> punc;
RegularExpAssignStm() { };
bool IdContains(char key);
size_t Search(size_t pos, string key, string match);
size_t SingleCharSearch(char key, size_t index, string match[]);
bool GetIdentifier(string assignStm);
bool GetLiteral(string assignStm);
bool Parse(string assignStm);
};
#include "pch.h"
#include "RegularExpAssignStm.h"
#include <functional>
#include <iostream>
#include <string>
using namespace std;
bool RegularExpAssignStm::IdContains(char key)
{
bool dg = false, uc = false, lc = false;
for (size_t i = 0; !uc && i < upperCase.size(); i++)
uc = key == upperCase[i];
if (uc)
return true;
for (size_t i = 0; !lc && i < lowerCase.size(); i++)
lc = key == lowerCase[i];
if (lc)
return true;
for (size_t i = 0; !dg && i < dig.size(); i++)
dg = key == dig[i];
if (dg)
return true;
return false;
}
size_t RegularExpAssignStm::Search(size_t pos, string key, string match)
{
bool found = false;
size_t i;
for (i = 0; !found && i < match.size(); i++)
found = key[pos] == match[i];
if (!found)
i = 4294967295;
else
i--;
return i;
}
size_t RegularExpAssignStm::SingleCharSearch(
char key, size_t index, string match[])
{
bool found = false;
size_t i;
for (i = 0; !found && i < match[index].size(); i++)
found = key == match[index].c_str()[i];
if (!found)
i = 4294967295;
else
i--;
return i;
}
bool RegularExpAssignStm::GetIdentifier(string assignStm)
{
string idStr;
for (size_t i = 0; i < dig.size(); i++)
if (assignStm[0] == dig[i])
return false;
for (size_t i = 0; i < assignStm.size(); i++)
{
if (IdContains(assignStm[i]))
idStr.push_back(assignStm[i]);
else
break;
}
if (idStr.size() > 0)
identifier.push_back(idStr);
return idStr.size() > 0;
}
bool RegularExpAssignStm::GetLiteral(string assignStm)
{
bool start = false;
string litStr;
for (size_t i = 0; !start && i < assignStm.size(); i++)
{
if (assignStm[0] == dig[i])
start = true;
}
if (start)
{
liter.push_back("");
for (size_t i = 0; i < assignStm.size(); i++)
{
if (assignStm[i] >= '0' && assignStm[i] <= '9')
liter[liter.size() - 1].push_back(assignStm[i]);
else
return liter.size() > 0;
}
}
return false;
}
bool RegularExpAssignStm::Parse(string assignStm)
{
if (GetIdentifier(assignStm))
assignStm.erase(0, identifier[identifier.size() - 1].size());
else
return false;
size_t assignOpPos = Search(0, assignStm, ops[4]);
if (assignOpPos != 4294967295)
{
assignStm.erase(0, ops[4].size());
oper.push_back(ops[4]);
}
else
return false;
while (true)
{
if (GetLiteral(assignStm))
{
assignStm.erase(0, liter[liter.size() - 1].size());
if (assignStm.size() <= 0)
return false;
}
else if (GetIdentifier(assignStm) &&
identifier[identifier.size() - 1].size() != 0)
{
assignStm.erase(0, identifier[identifier.size() - 1].size());
if (assignStm.size() <= 0)
return false;
}
size_t plusPos, minusPos, timesPos, divPos;
plusPos = SingleCharSearch(assignStm[0], 0, ops);
minusPos = SingleCharSearch(assignStm[0], 1, ops);
timesPos = SingleCharSearch(assignStm[0], 2, ops);
divPos = SingleCharSearch(assignStm[0], 3, ops);
if (plusPos != 4294967295)
{
oper.push_back(ops[0]);
assignStm.erase(0, 1);
}
else if (minusPos != 4294967295)
{
oper.push_back(ops[1]);
assignStm.erase(0, 1);
}
else if (timesPos != 4294967295)
{
oper.push_back(ops[2]);
assignStm.erase(0, 1);
}
else if (divPos != 4294967295)
{
oper.push_back(ops[3]);
assignStm.erase(0, 1);
}
else
return false;
if (assignStm.size() <= 0)
return false;
if (GetLiteral(assignStm))
{
assignStm.erase(0, liter[liter.size() - 1].size());
if (assignStm.size() <= 0)
return false;
}
else if (GetIdentifier(assignStm) &&
identifier[identifier.size() - 1].size() != 0)
{
assignStm.erase(0, identifier[identifier.size() - 1].size());
if (assignStm.size() <= 0)
return false;
}
size_t puns = SingleCharSearch(assignStm[0], 0, punctuation);
if (puns != 4294967295)
{
punc.push_back(punctuation[puns]);
assignStm.erase(0, punc[punc.size() - 1].size());
if (assignStm.size() <= 0)
return false;
}
}
return true;
}
Numerical Explorations 