Knight’s Tour with Windows C++ Mutual Exclusion and Threading by James Pate Williams, Jr.

Menu
1 N-Queens Problem (Backtracking)
2 Knight's Tour Problem (Backtracking)
3 Warnsdorff Heuristic Solution
4 Print Knight Move Count Matrix
5 Exit
Option (1 - 5): 2
n = 6
seed = 1
starting row = 1
starting col = 3
ending row = 1
ending col = 3
length = 36
maximum runtime per instance in seconds = 60
Chessboard Rows and Columns = 6
Pseudo-random Number Generator Seed = 1
Row = 1
Col = 3
Sequence Length = 36
18 29 26  9 12 35
27  8 19  0 25 10
30 17 28 11 34 13
 7 20  1 32  3 24
16 31 22  5 14 33
21  6 15  2 23  4
Current Runtime in Seconds = 50.747000
Menu
1 N-Queens Problem (Backtracking)
2 Knight's Tour Problem (Backtracking)
3 Warnsdorff Heuristic Solution
4 Print Knight Move Count Matrix
5 Exit
Option (1 - 5):
Menu
1 N-Queens Problem (Backtracking)
2 Knight's Tour Problem (Backtracking)
3 Warnsdorff Heuristic Solution
4 Print Knight Move Count Matrix
5 Exit
Option (1 - 5): 3
n = 6
seed = 1
starting row = 1
starting col = 3
ending row = 1
ending col = 3
length = 36
maximum runtime per instance in seconds = 60
Chessboard Rows and Columns = 6
Pseudo-random Number Generator Seed = 1
Row = 1
Col = 3
Sequence Length = 36
18 29 26  9 12 35
27  8 19  0 25 10
30 17 28 11 34 13
 7 20  1 32  3 24
16 31 22  5 14 33
21  6 15  2 23  4
Current Runtime in Seconds = 46.983000
Menu
1 N-Queens Problem (Backtracking)
2 Knight's Tour Problem (Backtracking)
3 Warnsdorff Heuristic Solution
4 Print Knight Move Count Matrix
5 Exit
Option (1 - 5):

// KnightsTourCPP.cpp : Defines the entry point for the console application.
//
/*
Author:	James Pate Williams, Jr. (c) 2000 - 2022
Backtracking algorithm applied to the N-Queens CSP.
Algorithm from _Foundations of Constraint Satisfaction_
by E. P. K. Tsang Chapter 2 page 37.
Warnsdorff's solution by James Pate Williams, Jr.
*/

#include "stdafx.h"
#include <windows.h>
#include <limits.h>
#include <stdlib.h>
#include <synchapi.h>
#include <time.h>
#include <algorithm>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <ostream>
#include <string>
#include <vector>
#include <chrono>
using namespace std::chrono;
using namespace std;

// https://www.geeksforgeeks.org/measure-execution-time-function-cpp/

const int lBufferLength = 8192;

class Pair {
public:
	int row, col;
	static int n;
	Pair() { row = col = -1; };
	Pair(int row, int col) {
		this->row = row;
		this->col = col;
	};
	inline friend bool operator ==(
		const Pair lt,
		const Pair rt) {
		return
			lt.row == rt.row &&
			lt.col == rt.col;
	};
	inline friend bool operator <(
		const Pair lt,
		const Pair rt) {
		int lval = lt.row * n + lt.col;
		int rval = rt.row * n + rt.col;
		return lval < rval;
	};
	inline friend bool operator >(
		const Pair lt,
		const Pair rt) {
		int lval = lt.row * n + lt.col;
		int rval = rt.row * n + rt.col;
		return lval > rval;
	};
	inline static int BinarySearch(
		vector<Pair> arr, int p, int r, Pair num) {
		if (r >= p) {
			int mid = p + (r - p) / 2;
			if (arr[mid] == num)
				return mid;
			if (arr[mid] > num)
				return BinarySearch(arr, p, mid - 1, num);
			if (arr[mid] < num)
				return BinarySearch(arr, mid + 1, r, num);
		}
		return -1;
	};
	inline static bool BinarySearchSTL(vector<Pair> solution, Pair test) {
		return binary_search(solution.begin(), solution.end(), test);
	};
	inline static bool LinearSearch(vector<Pair> solution, Pair test, int& index) {
		for (int i = 0; i < (int)solution.size(); i++) {
			if (test == solution[i]) {
				index = i;
				return true;
			}
		}
		index = -1;
		return false;
	};
	inline static bool FindSTL(vector<Pair> solution, Pair test, int& index)
	{
		vector<Pair>::iterator sit = find(solution.begin(), solution.end(), test);
		index = sit - solution.begin();
		return sit != solution.end();
	};
};

int Pair::n = 6;

class Board {
public:
	ofstream file;
	bool warnsdorff, writeToFile;
	char* buffer;
	double currentRuntime, cumulativeRuntime;
	int bufferLength;
	int cumulativeCount, length, n;
	int startingRow, startingCol;
	int endingRow, endingCol;
	int maxRuntimeSeconds;
	unsigned int seed;
	bool** marked;
	int** legalMoveCount;
	vector<int> unlabeled;
	vector<Pair> compoundLabel;
	vector<Pair> Dx;
	vector<Pair> solution;
	Board(
		bool warnsdorff, 
		bool writeToFile,
		char* buffer,
		int bufferLength,
		int n,
		int startingRow, 
		int startingCol,
		int endingRow,
		int endingCol,
		int length,
		int maxRuntimeSeconds,
		unsigned int seed) {
		this->warnsdorff = warnsdorff;
		this->writeToFile = writeToFile;
		this->buffer = buffer;
		this->bufferLength = bufferLength;
		this->n = n;
		this->startingRow = startingRow;
		this->startingCol = startingCol;
		this->endingRow = endingRow;
		this->endingCol = endingCol;
		this->length = length;
		this->maxRuntimeSeconds = maxRuntimeSeconds;
		this->seed = seed;
		Initialization();
		currentRuntime = 0;
		cumulativeCount = 0;
		cumulativeRuntime = 0;
	}
	~Board() {
		if (marked != NULL) {
			for (int i = 0; i < n; i++)
				delete[] marked[i];
			delete[] marked;
		}
		if (legalMoveCount != NULL) {
			for (int i = 0; i < n; i++)
				delete[] legalMoveCount[i];
			delete[] legalMoveCount;
		}
	}
	vector<Pair> GetLegalKnightMoves(
		int i, int j);
	void Initialization();
	int KTconstraintsViolated(Pair pair);
	bool KTBacktrackingSolution(
		int count, int prevx, int prevy,
		vector<int> unlabeled,
		vector<Pair>& compounLabel,
		vector<Pair> Dx,
		vector<Pair>& solution);
	void DoBTSolution(
		int row, int col,
		HANDLE hBTMutex);
	void Board::GetPossibleMoveCountMatrix();
	void PrintPossibleMoveCountMatrix();
	bool WarnsdorffSolution(
		int count, int prevx, int prevy,
		vector<int> unlabeled,
		vector<Pair>& compoundLabel,
		vector<Pair> Dx,
		vector<Pair>& solution);
	void DoWarnsdorffSolution(
		int row, int col,
		HANDLE hWDMutex);
	void PrintPairVector();
	bool WriteBoardDataToFile(
		bool writeSolution, int row, int col);
};

vector<Pair> Board::GetLegalKnightMoves(int i, int j)
{
	int count = 0;
	Pair pair;
	vector<Pair> legal;
	// quadrant 1
	if (i + 1 < n && j + 2 < n)
	{
		if (!marked[i + 1][j + 2])
		{
			pair.row = i + 1;
			pair.col = j + 2;
			legal.push_back(pair);
		}
	}
	if (i + 2 < n && j + 1 < n)
	{
		if (!marked[i + 2][j + 1])
		{
			pair.row = i + 2;
			pair.col = j + 1;
			legal.push_back(pair);
		}
	}
	// quadrant 2
	if (i - 1 >= 0 && j + 2 < n)
	{
		if (!marked[i - 1][j + 2])
		{
			pair.row = i - 1;
			pair.col = j + 2;
			legal.push_back(pair);
		}
	}
	if (i - 2 >= 0 && j + 1 < n)
	{
		if (!marked[i - 2][j + 1])
		{
			pair.row = i - 2;
			pair.col = j + 1;
			legal.push_back(pair);
		}
	}
	// quadrant 3
	if (i - 1 >= 0 && j - 2 >= 0)
	{
		if (!marked[i - 1][j - 2])
		{
			pair.row = i - 1;
			pair.col = j - 2;
			legal.push_back(pair);
		}
	}
	if (i - 2 >= 0 && j - 1 >= 0)
	{
		if (!marked[i - 2][j - 1])
		{
			pair.row = i - 2;
			pair.col = j - 1;
			legal.push_back(pair);
		}
	}
	// quadrant 4
	if (i + 1 < n && j - 2 >= 0)
	{
		if (!marked[i + 1][j - 2])
		{
			pair.row = i + 1;
			pair.col = j - 2;
			legal.push_back(pair);
		}
	}
	if (i + 2 < n && j - 1 >= 0)
	{
		if (!marked[i + 2][j - 1])
		{
			pair.row = i + 2;
			pair.col = j - 1;
			legal.push_back(pair);
		}
	}
	return legal;
}

void Board::Initialization() {
	if (marked == NULL) {
		marked = new bool* [n];
		for (int i = 0; i < n; i++) {
			marked[i] = new bool[n];
		}
	}
	if (legalMoveCount == NULL) {
		legalMoveCount = new int* [n];
		for (int i = 0; i < n; i++) {
			legalMoveCount[i] = new int[n];
		}
	}
	if (warnsdorff)
		GetPossibleMoveCountMatrix();
	srand(seed);
	unlabeled.clear();
	compoundLabel.clear();
	Dx.clear();
	solution.clear();
	for (int i = 0; i < n * n; i++)
		unlabeled.push_back(i);
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			Pair move;
			move.row = i;
			move.col = j;
			marked[i][j] = false;
			Dx.push_back(move);
		}
	}
}

int Board::KTconstraintsViolated(Pair pair) {
	int cv = 0;
	int row = pair.row;
	int col = pair.col;
	int xi = row, yi = col;

	for (int i = 0; i < (int)compoundLabel.size(); i++)
	{
		int x = compoundLabel[i].row;
		int y = compoundLabel[i].col;
		int r1 = x + 1, r2 = x + 2, r3 = x - 1, r4 = x - 2;
		int c1 = y + 1, c2 = y + 2, c3 = y - 1, c4 = y - 2;
		if (r1 < n && xi == r1 && c1 < n && yi == c1 ||
			r1 < n && xi == r1 && c2 < n && yi == c2 ||
			r1 < n && xi == r1 && c3 >= 0 && yi == c3 ||
			r1 < n && xi == r1 && c4 >= 0 && yi == c4 ||
			r2 < n && xi == r2 && c1 < n && yi == c1 ||
			r2 < n && xi == r2 && c2 < n && yi == c2 ||
			r2 < n && xi == r2 && c3 >= 0 && yi == c3 ||
			r2 < n && xi == r2 && c4 >= 0 && yi == c4 ||
			r3 >= 0 && xi == r3 && c1 < n && yi == c1 ||
			r3 >= 0 && xi == r3 && c2 < n && yi == c2 ||
			r3 >= 0 && xi == r3 && c3 >= 0 && yi == c3 ||
			r3 >= 0 && xi == r3 && c4 >= 0 && yi == c4 ||
			r4 >= 0 && xi == r4 && c1 < n && yi == c1 ||
			r4 >= 0 && xi == r4 && c2 < n && yi == c2 ||
			r4 >= 0 && xi == r4 && c3 >= 0 && yi == c3 ||
			r4 >= 0 && xi == r4 && c4 >= 0 && yi == c4)
		{
			if (marked[xi][yi])
				cv++;
		}
	}
	return cv;
}

bool Board::KTBacktrackingSolution(
	int count, int prevx, int prevy,
	vector<int> unlabeled,
	vector<Pair>& compounLabel,
	vector<Pair> Dx,
	vector<Pair>& solution)
{
	bool done, found, result;
	int i, x, y;
	Pair v;
	vector<Pair> legal;

	if (count == 0)
	{
		v.row = prevx;
		v.col = prevy;
		marked[prevx][prevy] = true;
		if (!Dx.empty()) {
			vector<Pair>::iterator dit = find(
				Dx.begin(), Dx.end(), v);
			if (dit != Dx.end())
				Dx.erase(dit);
		}
		if (!unlabeled.empty()) {
			vector<int>::iterator uit = find(
				unlabeled.begin(), unlabeled.end(), 0);
			if (uit != unlabeled.end())
				unlabeled.erase(uit);
		}
		compoundLabel.push_back(v);
		count++;
	}
	legal = GetLegalKnightMoves(
		prevx, prevy);
	while (!legal.empty())
	{
		x = rand() % legal.size();
		for (i = 0; i < (int)unlabeled.size(); i++)
		{
			if (x == unlabeled[i])
			{
				y = i;
				break;
			}
		}
		vector<int>::iterator uit = find(
			unlabeled.begin(), unlabeled.end(), y);
		if (uit != unlabeled.end()) {
			unlabeled.erase(uit);
		}
		done = false;
		do {
			v = legal[x];
			if (!marked[v.row][v.col]) {
				marked[v.row][v.col] = true;
			}
			int cv = KTconstraintsViolated(v);
			if (cv > 0) {
				found = binary_search(compoundLabel.begin(), compoundLabel.end(), v);
				if (!found) {
					vector<Pair>::iterator cit = find(
						compoundLabel.begin(), compoundLabel.end(), v);
					if (cit == compoundLabel.end())
						compoundLabel.push_back(v);
					vector<Pair>::iterator lit = find(
						legal.begin(), legal.end(), v);
					if (lit != legal.end()) {
						legal.erase(lit);
					}
					vector<Pair>::iterator dit = find(
						Dx.begin(), Dx.end(), v);
					if (dit != Dx.end())
						Dx.erase(dit);
					marked[v.row][v.col] = true;
					prevx = v.row;
					prevy = v.col;
					done = true;
				}
				if (compoundLabel.size() == length)
				{
					for (i = 0; i < (int)compoundLabel.size(); i++)
						solution.push_back(compoundLabel[i]);
					return true;
				}
				result = KTBacktrackingSolution(
					count, prevx, prevy,
					unlabeled, compoundLabel, Dx, solution);
				if (result)
					return result;
			}
		} while (!Dx.empty() && !legal.empty() && !done);
	}
	Pair move;
	prevx = compoundLabel[compoundLabel.size() - 1].row;
	prevy = compoundLabel[compoundLabel.size() - 1].col;
	move.row = prevx;
	move.col = prevy;
	marked[prevx][prevy] = false;
	vector<Pair>::iterator cit = find(
		compoundLabel.begin(), compoundLabel.end(), move);
	if (cit != compoundLabel.end()) {
		compoundLabel.erase(cit);
	}
	return false;
}

void Board::DoBTSolution(
	int row, int col,
	HANDLE hBTMutex) {
	int count = 0;
	int prevx = row;
	int prevy = col;
	startingRow = row;
	startingCol = col;
	bool io = false, bt = false;
	auto start = high_resolution_clock::now();
	bt = KTBacktrackingSolution(
		count, prevx, prevy,
		unlabeled, compoundLabel,
		Dx, solution);
	auto stop = high_resolution_clock::now();
	auto duration = duration_cast<milliseconds>(stop - start);
	currentRuntime = (double)duration.count() / 1.0e3;
	if (bt)
		io = WriteBoardDataToFile(true, row, col);
	else
		io = WriteBoardDataToFile(false, row, col);
	cumulativeCount++;
	cumulativeRuntime += currentRuntime;
	ReleaseMutex(hBTMutex);
}

void Board::GetPossibleMoveCountMatrix() {
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			vector<Pair> possible = GetLegalKnightMoves(
				i, j);
			legalMoveCount[i][j] = possible.size();
		}
	}
}

void Board::PrintPossibleMoveCountMatrix()
{
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			std::cout << legalMoveCount[i][j] << " ";
		}
		std::cout << endl;
	}
}

void Board::PrintPairVector() {
	int** matrix = new int* [n];
	for (int i = 0; i < n; i++)
		matrix[i] = new int[n];
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
			matrix[i][j] = -1;
	}
	bool done = (int)solution.size() == n * n;
	if (done) {
		for (int i = 0, count = 0; i < (int)solution.size(); count++, i++)
		{
			Pair pair = solution[i];
			std::cout << "(" << setw(2) << pair.row << ", ";
			std::cout << setw(2) << pair.col << ")" << " ";
			if ((count + 1) % n == 0)
				std::cout << endl;
			matrix[pair.row][pair.col] = count;
		}
		for (int i = 0; i < n; i++)
		{
			for (int j = 0; j < n; j++)
				std::cout << setw(2) << matrix[i][j] << " ";
			std::cout << endl;
		}
	}
	else {
		for (int i = 0, count = 0; i < (int)compoundLabel.size(); count++, i++)
		{
			Pair pair = compoundLabel[i];
			std::cout << "(" << setw(2) << pair.row << ", ";
			std::cout << setw(2) << pair.col << ")" << " ";
			if ((i + 1) % n == 0)
				std::cout << endl;
		}
		std::cout << endl;
	}
	if (matrix != NULL) {
		for (int i = 0; i < n; i++)
			delete[] matrix[i];
		delete[] matrix;
	}
}

bool Board::WarnsdorffSolution(
	int count, int prevx, int prevy,
	vector<int> unlabeled,
	vector<Pair>& compoundLabel,
	vector<Pair> Dx,
	vector<Pair>& solution) {
	bool done, found, result;
	int i;
	Pair v;
	vector<Pair> legalMove, nextMoveVector, nextMove;

	if (count == 0)
	{
		prevx = startingRow;
		prevy = startingCol;
		v.row = prevx;
		v.col = prevy;
		marked[v.row][v.col] = true;
		if (!Dx.empty()) {
			vector<Pair>::iterator dit = find(
				Dx.begin(), Dx.end(), v);
			if (dit != Dx.end())
				Dx.erase(dit);
		}
		if (!unlabeled.empty()) {
			vector<int>::iterator uit = find(
				unlabeled.begin(), unlabeled.end(), 0);
			if (uit != unlabeled.end())
				unlabeled.erase(uit);
		}
		compoundLabel.push_back(v);
		count++;
	}
	legalMove = GetLegalKnightMoves(
		prevx, prevy);
	while (!legalMove.empty())
	{
		int min = INT_MAX;
		nextMoveVector.clear();
		for (int k = 0; k < (int)legalMove.size(); k++) {
			if (legalMoveCount[legalMove[k].row][legalMove[k].col] < min) {
				min = legalMoveCount[legalMove[k].row][legalMove[k].col];
				nextMoveVector.push_back(legalMove[k]);
			}
		}
		vector<int> index;
		index.clear();
		for (int k = 0; k < (int)nextMoveVector.size(); k++) {
			int number = legalMoveCount[nextMoveVector[k].row][nextMoveVector[k].col];
			if (number == min) {
				index.push_back(k);
			}
		}
		int r = rand() % index.size();
		v = nextMoveVector[r];
		done = false;
		do {
			if (!marked[v.row][v.col]) {
				marked[v.row][v.col] = true;
				nextMove.push_back(v);
			}
			int sum = 0;
			for (int i = 0; i < n; i++) {
				for (int j = 0; j < n; j++) {
					if (marked[i][j])
						sum++;
				}
			}
			int cv = KTconstraintsViolated(v);
			if (cv > 0) {
				found = binary_search(compoundLabel.begin(), compoundLabel.end(), v);
				if (!found) {
					vector<Pair>::iterator cit = find(
						compoundLabel.begin(), compoundLabel.end(), v);
					if (cit == compoundLabel.end())
						compoundLabel.push_back(v);
					vector<Pair>::iterator lit = find(
						legalMove.begin(), legalMove.end(), v);
					if (lit != legalMove.end()) {
						legalMove.erase(lit);
					}
					vector<Pair>::iterator dit = find(
						Dx.begin(), Dx.end(), v);
					if (dit != Dx.end())
						Dx.erase(dit);
					marked[v.row][v.col] = true;
					prevx = v.row;
					prevy = v.col;
					done = true;
				}
				if (compoundLabel.size() == length)
				{
					for (i = 0; i < (int)compoundLabel.size(); i++)
						solution.push_back(compoundLabel[i]);
					return true;
				}
				result = WarnsdorffSolution(
					count, prevx, prevy, unlabeled,
					compoundLabel, Dx, solution);
				if (result)
					return result;
			}
		} while (!Dx.empty() && !legalMove.empty() && !done);
	}
	Pair move;
	prevx = compoundLabel[compoundLabel.size() - 1].row;
	prevy = compoundLabel[compoundLabel.size() - 1].col;
	move.row = prevx;
	move.col = prevy;
	marked[prevx][prevy] = false;
	vector<Pair>::iterator cit = find(
		compoundLabel.begin(), compoundLabel.end(), move);
	if (cit != compoundLabel.end()) {
		compoundLabel.erase(cit);
	}
	return false;
}

void Board::DoWarnsdorffSolution(
	int row, int col,
	HANDLE hWDMutex) {
	int count = 0;
	int prevx = row;
	int prevy = col;
	startingRow = row;
	startingCol = col;
	bool io = false, ws = false;
	auto start = high_resolution_clock::now();
	Initialization();
	ws = WarnsdorffSolution(
		count, prevx, prevy,
		unlabeled, compoundLabel,
		Dx, solution);
	auto stop = high_resolution_clock::now();
	auto duration = duration_cast<milliseconds>(stop - start);
	currentRuntime = (double)duration.count() / 1.0e3;
	if (ws)
		io = WriteBoardDataToFile(true, row, col);
	else
		io = WriteBoardDataToFile(false, row, col);
	cumulativeCount++;
	cumulativeRuntime += currentRuntime;
	ReleaseMutex(hWDMutex);
}

bool Board::WriteBoardDataToFile(
	bool writeSolution, int row, int col) {
	char newLine[2] = { '\n', '\0' };
	errno_t err;
	if (strlen(buffer) == 0)
		strcpy_s(buffer, bufferLength, "Chessboard Rows and Columns = ");
	else
		strcat_s(buffer, bufferLength, "Chessboard Rows and Columns = ");
	char nBuf[64] = { '\0' };
	err = _itoa_s(n, nBuf, 10);
	strcat_s(buffer, bufferLength, nBuf);
	strcat_s(buffer, bufferLength, newLine);
	strcat_s(buffer, bufferLength, "Pseudo-random Number Generator Seed = ");
	err = _itoa_s((int)seed, nBuf, 10);
	strcat_s(buffer, bufferLength, nBuf);
	strcat_s(buffer, bufferLength, newLine);
	strcat_s(buffer, bufferLength, "Row = ");
	err = _itoa_s(row, nBuf, 10);
	strcat_s(buffer, bufferLength, nBuf);
	strcat_s(buffer, bufferLength, newLine);
	strcat_s(buffer, bufferLength, "Col = ");
	err = _itoa_s(col, nBuf, 10);
	strcat_s(buffer, bufferLength, nBuf);
	strcat_s(buffer, bufferLength, newLine);
	strcat_s(buffer, bufferLength, "Sequence Length = ");
	err = _itoa_s(length, nBuf, 10);
	strcat_s(buffer, bufferLength, nBuf);
	strcat_s(buffer, bufferLength, newLine);
	if (writeSolution && solution.size() == n * n)
	{
		int** matrix = new int* [n];
		for (int i = 0; i < n; i++)
			matrix[i] = new int[n];
		for (int i = 0; i < n; i++)
		{
			for (int j = 0; j < n; j++)
				matrix[i][j] = -1;
		}
		for (int i = 0, count = 0; i < (int)solution.size(); count++, i++)
		{
			Pair pair = solution[i];
			int r = pair.row;
			int c = pair.col;
			matrix[pair.row][pair.col] = count;
		}
		for (int i = 0; i < n; i++)
		{
			for (int j = 0; j < n; j++) {
				string space = " ";
				if (matrix[i][j] < 10)
					err = strcat_s(buffer, bufferLength, space.c_str());
				err = _itoa_s(matrix[i][j], nBuf, 10);
				err = strcat_s(buffer, bufferLength, nBuf);
				err = strcat_s(buffer, bufferLength, space.c_str());
			}
			strcat_s(buffer, bufferLength, newLine);
		}
		strcat_s(buffer, bufferLength, "Current Runtime in Seconds = ");
		string cr = to_string(currentRuntime);
		strcat_s(buffer, bufferLength, cr.c_str());
		strcat_s(buffer, bufferLength, newLine);
		if (matrix != NULL) {
			for (int i = 0; i < n; i++)
				delete[] matrix[i];
			delete[] matrix;
		}
	}
	return true;
}

typedef struct Data
{
	int row, col;
	Board* board;
	DWORD dwWaitMS;
} DATA, * PDATA;

DWORD WINAPI BTThreadCode(LPVOID lpParam)
{
	Data* pData = (Data*)lpParam;
	HANDLE hBTMutex = CreateMutex(NULL, FALSE, NULL);

	if (hBTMutex != NULL)
	{
		pData->board->DoBTSolution(
			pData->row,
			pData->col,
			hBTMutex);
		DWORD err = WaitForSingleObject(hBTMutex, pData->dwWaitMS);
		CloseHandle(hBTMutex);
	}
	return 0;
}

DWORD WINAPI WDThreadCode(LPVOID lpParam)
{
	Data* pData = (Data*)lpParam;
	HANDLE hWDMutex = CreateMutex(NULL, FALSE, NULL);

	if (hWDMutex != NULL)
	{
		pData->board->DoBTSolution(
			pData->row,
			pData->col,
			hWDMutex);
		DWORD err = WaitForSingleObject(
			hWDMutex, pData->dwWaitMS);
		CloseHandle(hWDMutex);
	}
	return 0;
}

int NQconstraintsViolated(vector<int>& Q) {
	int a, b, c = 0, i, j, n;

	n = Q.size();
	for (i = 0; i < n; i++) {
		a = Q[i];
		for (j = 0; j < n; j++) {
			b = Q[j];
			if (i != j && a != -1 && b != -1) {
				if (a == b) c++;
				if (i - j == a - b || i - j == b - a) c++;
			}
		}
	}
	return c;
}

bool NQBacktrack(
	vector<int> unlabeled,
	vector<int> compoundLabel,
	vector<int>& solution) {
	bool result;
	int i, j, v, x;
	vector<int> Dx(compoundLabel.size());

	if (unlabeled.empty()) {
		for (i = 0; i < (int)compoundLabel.size(); i++)
			solution[i] = compoundLabel[i];
		return true;
	}
	for (i = 0; i < (int)compoundLabel.size(); i++)
		Dx[i] = i;
	// pick a random variable from unlabeled array
	i = rand() % unlabeled.size();
	x = unlabeled[i];
	do {
		// pick a random value from domain of x
		j = rand() % Dx.size();
		v = Dx[j];
		vector<int>::iterator vIterator = find(Dx.begin(), Dx.end(), v);
		Dx.erase(vIterator);
		compoundLabel[x] = v;
		if (NQconstraintsViolated(compoundLabel) == 0) {
			vIterator = find(unlabeled.begin(), unlabeled.end(), x);
			unlabeled.erase(vIterator);
			result = NQBacktrack(unlabeled, compoundLabel, solution);
			if (result) return result;
			unlabeled.push_back(x);
		}
		else
			compoundLabel[x] = -1;
	} while (!Dx.empty());
	return false;
}

void PrintNQSolution(vector<int>& solution) {
	char hyphen[256] = { 0 };
	int column, i, i4, n = solution.size(), row;

	if (n <= 10) {
		for (i = 0; i < n; i++) {
			i4 = i * 4;
			hyphen[i4 + 0] = '-';
			hyphen[i4 + 1] = '-';
			hyphen[i4 + 2] = '-';
			hyphen[i4 + 3] = '-';
		}
		i4 = i * 4;
		hyphen[i4 + 0] = '-';
		hyphen[i4 + 1] = '\n';
		hyphen[i4 + 2] = '\0';
		for (row = 0; row < n; row++) {
			column = solution[row];
			std::cout << hyphen;
			for (i = 0; i < column; i++)
				std::cout << "|   ";
			std::cout << "| Q ";
			for (i = column + 1; i < n; i++)
				std::cout << "|   ";
			std::cout << '|' << endl;
		}
		std::cout << hyphen;
	}
	else
		for (row = 0; row < n; row++)
			std::cout << row << ' ' << solution[row] << endl;
	std::cout << endl;
}

int main(int argc, char* argv[]) {
	std::cout.setf(ios::fixed, ios::floatfield);
	std::cout.setf(ios::showpoint);
	while (true) {
		bool valid;
		int i, length, n, option;
		int srow = 0, scol = 0;
		int erow = 0, ecol = 0;
		int maxRuntimeSeconds = 0;
		unsigned int seed;
		std::cout << "Menu" << endl;
		std::cout << "1 N-Queens Problem (Backtracking)" << endl;
		std::cout << "2 Knight's Tour Problem (Backtracking)" << endl;
		std::cout << "3 Warnsdorff Heuristic Solution" << endl;
		std::cout << "4 Print Knight Move Count Matrix" << endl;
		std::cout << "5 Exit" << endl;
		std::cout << "Option (1 - 5): ";
		cin >> option;
		if (option >= 5)
			return 0;
		valid = false;
		while (!valid) {
			valid = true;
			std::cout << "n = ";
			cin >> n;
			while (n < 5 || n > 8) {
				std::cout << "n must be >= 5 and <= 8" << endl;
				std::cout << "n = ";
				cin >> n;
			}
			if (option != 4) {
				std::cout << "seed = ";
				cin >> seed;
				while (seed < 1) {
					std::cout << "seed must be >= 1" << endl;
					std::cout << "seed = ";
					cin >> seed;
				}
			}
			if (option == 2 || option == 3) {
				std::cout << "starting row = ";
				cin >> srow;
				while (srow < 0 || srow > 7) {
					std::cout << "starting row must be >= 0 and <= 9" << endl;
					std::cout << "starting row = ";
					cin >> srow;
				}
				std::cout << "starting col = ";
				cin >> scol;
				while (scol < 0 || scol > 7) {
					std::cout << "starting col must be >= 0 and <= 9" << endl;
					std::cout << "starting col = ";
					cin >> scol;
				}
				std::cout << "ending row = ";
				cin >> erow;
				while (erow < 0 || erow > 7) {
					std::cout << "starting row must be >= 0 and <= 7" << endl;
					std::cout << "starting row = ";
					cin >> erow;
				}
				std::cout << "ending col = ";
				cin >> ecol;
				while (ecol < 0 || ecol > 7) {
					std::cout << "ending col must be >= 0 and <= 7" << endl;
					std::cout << "ending col = ";
					cin >> ecol;
				}
			}
			if (option == 2 || option == 3) {
				std::cout << "length = ";
				cin >> length;
				while (length < 25 || length > 64)
				{
					std::cout << "length must be >= 25 and <= 64" << endl;
					std::cout << "length = ";
					cin >> length;
				}
				if (n == 5 && length < 25) {
					std::cout << "Invalid length value for given n" << endl;
					std::cout << "length  must be <= 25" << endl;
					valid = false;
				}
				else if (n == 6) {
					if (length > 36) {
						std::cout << "Invalid length value for given n" << endl;
						std::cout << "length must be <= 36" << endl;
						valid = false;
					}
				}
				else if (n == 7) {
					if (length > 49) {
						std::cout << "Invalid length value for given n" << endl;
						std::cout << "length must be <= 49" << endl;
						valid = false;
					}
				}
				else if (n == 8) {
					if (length > 64) {
						std::cout << "Invalid length value for given n" << endl;
						std::cout << "length must be <= 64" << endl;
						valid = false;
					}
				}
				std::cout << "maximum runtime per instance in seconds = ";
				cin >> maxRuntimeSeconds;
				while (maxRuntimeSeconds <= 0) {
					std::cout << "maximum runtime > 0" << endl;
					std::cout << "maximum runtime per instance in seconds = ";
					cin >> maxRuntimeSeconds;
				}
			}
		}
		if (option == 1) {
			auto start = high_resolution_clock::now();
			vector<int> compoundLabel(n, -1), solution(n, -1), unlabeled(n);
			for (i = 0; i < n; i++)
				unlabeled[i] = i;
			if (NQBacktrack(unlabeled, compoundLabel, solution))
				PrintNQSolution(solution);
			else
				std::cout << "Problem has no solution" << endl;
			auto stop = high_resolution_clock::now();
			auto duration = duration_cast<microseconds>(stop - start);
			std::cout << "Runtime: " << setprecision(6)
				<< (double)duration.count() / 1.0e6 << " seconds" << endl;
		}
		else if (option == 2) {
			string filename = "Chronological Backtracking Results.txt";
			ofstream file;
			file.open(filename.c_str(), ios::ate);
			Pair::n = n;
			for (int i = srow; i <= erow; i++) {
				for (int j = scol; j <= ecol; j++) {
					char lBuffer[lBufferLength] = { 0 };
					Board* board = new Board(
						true,
						true,
						lBuffer,
						lBufferLength,
						n,
						srow,
						scol,
						erow,
						ecol,
						length,
						maxRuntimeSeconds,
						seed);
					Data* pBTData = new Data();
					DWORD dwThreadId;
					pBTData->row = i;
					pBTData->col = j;
					pBTData->board = board;
					pBTData->dwWaitMS = maxRuntimeSeconds * 1000;
					HANDLE hBTThread = CreateThread(
						NULL,			    // default security attributes
						0,					// use default stack size  
						BTThreadCode,		// thread function name
						pBTData,			// argument to thread function 
						0,			        // use default creation flags 
						&dwThreadId);		// returns the thread identifier
					if (hBTThread != NULL)
					{
						WaitForSingleObject(
							hBTThread, maxRuntimeSeconds * 1000);
						CloseHandle(hBTThread);
					}
					file << lBuffer;
					std::cout << lBuffer;
				}
			}
		}
		else if (option == 3) {
			string filename = "Warnsdorff Heuristic Results.txt";
			ofstream file;
			file.open(filename.c_str(), ios::ate);
			Pair::n = n;
			for (int i = srow; i <= erow; i++) {
				for (int j = scol; j <= ecol; j++) {
					char lBuffer[lBufferLength] = { 0 };
					Board* board = new Board(
						true,
						true,
						lBuffer,
						lBufferLength,
						n,
						srow,
						scol,
						erow,
						ecol,
						length,
						maxRuntimeSeconds,
						seed);
					Data* pWDData = new Data();
					DWORD dwThreadId;
					pWDData->row = i;
					pWDData->col = j;
					pWDData->board = board;
					pWDData->dwWaitMS = maxRuntimeSeconds * 1000;
					HANDLE hWDThread = CreateThread(
						NULL,			    // default security attributes
						0,					// use default stack size  
						WDThreadCode,		// thread function name
						pWDData,			// argument to thread function 
						0,			        // use default creation flags 
						&dwThreadId);		// returns the thread identifier
					if (hWDThread != NULL)
					{
						WaitForSingleObject(
							hWDThread, maxRuntimeSeconds * 1000);
						CloseHandle(hWDThread);
					}
					file << lBuffer;
					std::cout << lBuffer;
				}
			}
		}
		else if (option == 4) {
			Board board(
				true,
				true,
				NULL,
				0,
				n,
				0,
				0,
				0,
				0,
				0,
				0,
				1);
			board.GetPossibleMoveCountMatrix();
			board.PrintPossibleMoveCountMatrix();
		}
	}
	return 0;
}

Sorting C++ Source Code by James Pate Williams, Jr.

#include "SampleStatistics.h"
#include "Sorting.h"
#include <stdlib.h>
#include <iostream>
#include <iomanip>
#include <vector>
#include <chrono>
using namespace std::chrono;
using namespace std;

// https://www.geeksforgeeks.org/measure-execution-time-function-cpp/

char GetSubOption(char title[])
{
	char line[128];
	cout << title << endl;
	cout << "1 Insertion Sort" << endl;
	cout << "2 Selection Sort" << endl;
	cout << "3 Heap Sort" << endl;
	cout << "4 Quick Sort" << endl;
	cout << "5 Singleton's Sort" << endl;
	cout << "6 Tree Sort 3" << endl;
	cout << "7 Exit Submenu" << endl;
	cout << "8 Exit Program" << endl;
	cout << "Suboption Number: ";
	cin.getline(line, 128);
	return line[0];
}

char DoSort(char option) {
	char subOption, line[128];
	char title[128] = { 0 };
	char name[6][64] = { { 0 } };
	strcpy_s(name[0], 64, "Insertion Sort");
	strcat_s(name[1], 64, "Selection Sort");
	strcat_s(name[2], 64, "Heap Sort");
	strcat_s(name[3], 64, "Quick Sort");
	strcat_s(name[4], 64, "Singleton's Sort");
	strcat_s(name[5], 64, "Tree Sort 3");
	vector<double> runtime;
	if (option == '1')
	{
		strcpy_s(title, 128, "Single Sorting Tests Submenu");
	}
	else {
		strcpy_s(title, 128, "Statistical Sorting Tests Submenu");
	}
	subOption = GetSubOption(title);
	if (subOption < '1' || subOption > '8') {
		cout << "Illegal Suboption Exiting Application!" << endl;
		exit(-1);
	}
	if (subOption == '7' || subOption == '8') {
		return subOption;
	}
	int index = (int)(subOption - '1');
	cout << name[index] << endl;
	cout << "PRNG Seed             = ";
	cin.getline(line, 128);
	int seed = atoi(line);
	cout << "Number of Samples     = ";
	cin.getline(line, 128);
	int n = atoi(line);
	cout << "Maximum Sample Value  = ";
	cin.getline(line, 128);
	int maximum = atoi(line);
	srand(seed);
	int number;
	if (option == '2') {
		cout << "Number of Experiments = ";
		cin.getline(line, 128);
		number = atoi(line);
	}
	else {
		number = 1;
	}
	double* sample = new double[n];
	double* shadow = new double[n];
	for (int i = 0; i < number; i++) {
		for (int j = 0; j < n; j++) {
			sample[j] = rand() % maximum;
			if (option == '1') {
				shadow[j] = sample[j];
			}
		}
		auto start = high_resolution_clock::now();
		switch (subOption) {
		case '1':
			Sorting::InsertionSort(sample, n);
			break;
		case '2':
			Sorting::SelectionSort(sample, n);
			break;
		case '3':
			Sorting::HeapSort(sample, n);
			break;
		case '4':
			Sorting::QuickSort(sample, n);
			break;
		case '5':
			Sorting::SingletonsSort(sample, n);
			break;
		case '6':
			Sorting::TreeSort3(sample, n);
			break;
		default:
			cout << "Unknown Suboption" << endl;
			break;
		}
		auto stop = high_resolution_clock::now();
		auto duration = duration_cast<microseconds>(stop - start);
		double rt = (double)duration.count();
		runtime.push_back(rt);
	}
	if (option == '1') {
		for (int i = 0; i < n; i++) {
			cout << setw(2) << shadow[i] << ' ';
			cout << setw(2) << sample[i] << endl;
		}
	}
	else if (option == '2') {
		cout << "Runtimes in Microseconds" << endl;
		cout << "Minimum Runtime       = " 
			<< setw(4) << fixed << setprecision(0)
			<< SampleStatistics::Min(runtime) << endl;
		cout << "Maximum Runtime       = " 
			<< setw(4) << fixed << setprecision(0)
			<< SampleStatistics::Max(runtime) << endl;
		cout << "Mean Runtime          = " 
			<< setw(4) << fixed << setprecision(0)
			<< SampleStatistics::Mean(runtime) << endl;
		cout << "Median Runtime        = "
			<< setw(4) << fixed << setprecision(0) <<
			SampleStatistics::Median(runtime) << endl;
	}
	delete[] sample;
	delete[] shadow;
	return subOption;
}

int main() {
	while (true) {
		char line[128];
		cout << "Menu" << endl;
		cout << "1 Single Sortimg Tests" << endl;
		cout << "2 Statistical Tests" << endl;
		cout << "3 Exit" << endl;;
		cout << "Option number: ";
		cin.getline(line, 128);
		if (line[0] == '3')
			break;
		if (line[0] == '1' || line[0] == '2') {
			while (true)
			{
				char doSort = DoSort(line[0]);
				if (doSort == '7')
					break;
				else if (doSort == '8')
					exit(1);
			}
		}
	}
	return 0;
}
#pragma once
#include <limits.h>
#include <vector>
using namespace std;

class SampleStatistics {
public:
	inline static double Mean(vector<double> x) {
		double sum = 0;

		for (int i = 0; i < (int)x.size(); i++)
			sum += (double)x[i];

		return sum / (int)x.size();
	};
	inline static double Max(vector<double> x)
	{
		double max = DBL_MIN;

		for (int i = 0; i < (int)x.size(); i++)
			if (x[i] > max)
				max = x[i];

		return max;
	};
	inline static double Min(vector<double> x)
	{
		double min = DBL_MAX;

		for (int i = 0; i < (int)x.size(); i++)
			if (x[i] < min)
				min = x[i];

		return min;
	};
	inline static double Median(vector<double> x)
	{
		int n = x.size(), n2 = n / 2 - 1;

		if (n % 2 == 1)
			return x[(n + 1) / 2 - 1];

		return 0.5 * (x[n2] + x[n2 + 1]);
	};
	inline static double Variance(vector<double> x)
	{
		double mean = Mean(x), sumSq = 0;
		int n = x.size();
		int n1 = n - 1;

		for (int i = 0; i < n; i++)
		{
			double delta = x[i] - mean;

			sumSq += delta * delta;
		}

		return sumSq / n1;
	};
};
#pragma once
// Insertion sort, heap sort, and quick sort algorithms
// From "Algortihms" by Thomas H. Cormen, et. al.
// Selection sort, Singleton's sort, and Tree Sort 3
// From "Sorting and Sort Systems" by Harold Lorin
class Sorting {
public:
	static void InsertionSort(double a[], int n);
	static void SelectionSort(double a[], int n);
	static void HeapSort(double a[], int n);
	static void QuickSort(double a[], int n);
	static void SingletonsSort(double a[], int n);
	static void TreeSort3(double a[], int n);
private:
	// five helper methods for Heap Sort
	static int Parent(int i);
	static int Left(int i);
	static int Right(int i);
	static void BuildHeap(double a[], int n, int& heapSize);
	static void Heapify(double a[], int i, int n, int& heapSize);
	// helper methods for Quick Sort
	static int Partition(double a[], int n, int lo, int hi);
	static void DoQuickSort(double a[], int n, int p, int r);
	// helper function for Singleton's sort
	static void DoSingletonsSort(double a[], int ii, int jj);
	// helper function for TreeSort3
	static void SiftUp(double a[], int i, int n);
};
#include "Sorting.h"

void Sorting::InsertionSort(double a[], int n) {
    for (int j = 1; j < n; j++)
    {
        double key = a[j];
        int i = j - 1;

        while (i >= 0 && a[i] > key)
        {
            a[i + 1] = a[i];
            i--;
        }

        a[i + 1] = key;
    }
}

void Sorting::SelectionSort(double a[], int n) {
    for (int i = 0; i < n - 1; i++)
    {
        for (int j = i + 1; j < n; j++)
        {
            if (a[i] > a[j])
            {
                double t = a[i];

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

int Sorting::Parent(int i) {
    return i / 2;
}

int Sorting::Left(int i) {
    return 2 * i;
}

int Sorting::Right(int i) {
    return 2 * i + 1;
}

void Sorting::Heapify(double a[], int i, int n, int& heapSize) {
    int l = Left(i);
    int r = Right(i);
    int largest = -1;

    if (l < heapSize && a[l] > a[i])
        largest = l;
    else
        largest = i;

    if (r < heapSize && a[r] > a[largest])
        largest = r;

    if (largest != i)
    {
        double t = a[i];

        a[i] = a[largest];
        a[largest] = t;

        Heapify(a, largest, n, heapSize);
    }
}

void Sorting::BuildHeap(double a[], int n, int& heapSize) {
    for (int i = n / 2; i >= 0; i--)
        Heapify(a, i, n, heapSize);
}

void Sorting::HeapSort(double a[], int n) {
    int heapSize = n;

    BuildHeap(a, n, heapSize);

    for (int i = n - 1; i >= 0; i--)
    {
        double t = a[0];

        a[0] = a[i];
        a[i] = t;

        heapSize--;

        Heapify(a, 0, n, heapSize);
    }
}

int Sorting::Partition(double a[], int n, int lo, int hi) {
    int pivotIndex = lo + (hi - lo) / 2;
    double x = a[pivotIndex];
    double t = x;

    a[pivotIndex] = a[hi];
    a[hi] = t;

    int storeIndex = lo;

    for (int i = lo; i < hi; i++)
    {
        if (a[i] < x)
        {
            t = a[i];
            a[i] = a[storeIndex];
            a[storeIndex++] = t;
        }
    }

    t = a[storeIndex];
    a[storeIndex] = a[hi];
    a[hi] = t;

    return storeIndex;
}

void Sorting::DoQuickSort(double a[], int n, int p, int r) {
    if (p < r)
    {
        int q = Partition(a, n, p, r);

        DoQuickSort(a, n, p, q - 1);
        DoQuickSort(a, n, q + 1, r);
    }
}

void Sorting::QuickSort(double a[], int n) {
    DoQuickSort(a, n, 0, n - 1);
}

void Sorting::DoSingletonsSort(double a[], int ii, int jj) {
    int m = 0;
    int i = ii;
    int j = jj;
    int iu[50] = { 0 };
    int il[50] = { 0 };
    int ij = 0, l = 1, k = 0;
    double t = 0.0, tt = 0.0;
    Label1:
        if (i >= j)
            goto Label8;
        goto Label2;
    Label2:
        k = i;
        ij = (j + i) / 2;
        t = a[ij];
        if (a[i] <= t)
            goto Label3;
        a[ij] = a[i];
        a[i] = t;
        t = a[ij];
        goto Label3;
    Label3:
        l = j;
        if (a[j] >= t)
            goto Label5;
        a[ij] = a[j];
        a[j] = t;
        t = a[ij];
        if (a[i] <= t)
            goto Label5;
        a[ij] = a[i];
        a[i] = t;
        t = a[ij];
        goto Label5;
    Label4:
        a[l] = a[k];
        a[k] = tt;
        goto Label5;
    Label5:
        l--;
        if (a[l] > t)
            goto Label5;
        tt = a[l];
        goto Label6;
    Label6:
        k++;
        if (a[k] < t)
            goto Label6;
        if (k <= l)
            goto Label4;
        if (l - i <= j - k)
            goto Label7;
        il[m] = i;
        iu[m] = l;
        i = k;
        m++;
        goto Label9;
    Label7:
        il[m] = k;
        iu[m] = j;
        j = l;
        m++;
        goto Label9;
    Label8:
        m--;
        if (m == -1)
            return;
        i = il[m];
        j = iu[m];
        goto Label9;
    Label9:
        if (j - i > 10)
            goto Label2;
        if (i == ii)
            goto Label1;
        i--;
        goto Label10;
    Label10:
        i++;
        if (i == j)
            goto Label8;
        t = a[i + 1];
        if (a[i] <= t)
            goto Label10;
        k = i;
        goto Label11;
    Label11:
        a[k + 1] = a[k];
        k--;
        if (t < a[k])
            goto Label11;
        a[k + 1] = t;
        goto Label10;
}

void Sorting::SingletonsSort(double a[], int n) {
    DoSingletonsSort(a, 0, n - 1);
}

void Sorting::SiftUp(double a[], int i, int n) {
    double copy = a[i];
    int j;

    while (true)
    {
        j = i + i;

        if (j <= n)
        {
            if (j < n)
            {
                if (a[j + 1] > a[j])
                    j++;
            }

            if (a[j] > copy)
            {
                a[i] = a[j];
                i = j;
                continue;
            }
        }

        break;
    }

    a[i] = copy;
}

void Sorting::TreeSort3(double a[], int n)
{
    for (int i = n / 2; i >= 1; i--)
        SiftUp(a, i, n - 1);

    for (int i = n - 1; i >= 1; i--)
    {
        SiftUp(a, 0, i);

        double t = a[0];

        a[0] = a[i];
        a[i] = t;
    }
}

Python Root Finder for Linear, Quadratic, Cubic, and Quartic Equations by James Pate Williams, Jr.

z1 = complex(0, 0)
z2 = complex(0, 0)
z3 = complex(0, 0)
z4 = complex(0, 0)

def complex_quadratic_formula(a, b, c):
    
    import cmath

    z1 = complex((-b + cmath.sqrt(b * b - 4 * a * c)) / (2 * a))
    z2 = complex((-b - cmath.sqrt(b * b - 4 * a * c)) / (2 * a))
    
    print("z1 = ", z1)
    print("z2 = ", z2)

def real_quadratic_formula(a, b, c):
    if b * b - 4 * a * c >= 0:
        
        import math
    
        z1 = float((-b + math.sqrt(b * b - 4 * a * c)) / (2 * a))
        z2 = float((-b - math.sqrt(b * b - 4 * a * c)) / (2 * a))

    else:
        
        import cmath

        z1 = complex((-b + cmath.sqrt(b * b - 4 * a * c)) / (2 * a))
        z2 = complex((-b - cmath.sqrt(b * b - 4 * a * c)) / (2 * a))
    
    print("z1 = ", z1)
    print("z2 = ", z2)

def cubic_equation(A, B, C, D):
    # https://en.wikipedia.org/wiki/Cubic_function
    # http://www.wolframalpha.com/widgets/view.jsp?id=3f4366aeb9c157cf9a30c90693eafc55

    import cmath
    import math

    A2 = float(A * A)
    B2 = float(B * B)
    B3 = float(B * B * B)
    C2 = float(C * C)
    C3 = float(C * C * C)
    D2 = float(D * D)
    Delta = float(18 * A * B * C * D - 4.0 * B3 * D
                + B2 * C2 - 4.0 * A * C3 - 27.0 * A2 * D2)
    Delta0 = float(B2 - 3.0 * A * C)
    Delta1 = float(2.0 * B3 - 9.0 * A * B * C + 27.0 * A2 * D)
    r = cmath.sqrt(-27.0 * A2 * Delta)
    c = complex(0.5 * (Delta1 + r)) ** (1.0 / 3.0)
    d = complex(0.5 * (Delta1 - r)) ** (1.0 / 3.0)

    if c.real == 0 and c.imag == 0:
        c = d

    i = complex(0.0, 1.0)
    chi = complex(-0.5, 0) + i * 0.5 * math.sqrt(3)

    import array

    ar = array.array('f', [0,0,0])
    ai = array.array('f', [0,0,0])

    for k in range(0, 3, 1):
        chik = chi ** k
        z = -(B + chik * c + Delta0 / (chik * c)) / (3.0 * A)
        ar[k] = z.real
        ai[k] = z.imag;

    z1 = ar[0] + ai[0] * complex(0, 1);
    z2 = ar[1] + ai[1] * complex(0, 1);
    z3 = ar[2] + ai[2] * complex(0, 1);
   
    print("z1 = ", z1)
    print("z2 = ", z2)
    print("z3 = ", z3)

def quartic_equation(A, B, C, D, E):
    # https://en.wikipedia.org/wiki/Quartic_function
    # https://keisan.casio.com/exec/system/1181809416

    import cmath

    A2 = float(A * A)
    A3 = float(A * A2)
    A4 = float(A2 * A2)
    B2 = float(B * B)
    B3 = float(B * B2)
    B4 = float(B2 * B2)
    C2 = float(C * C)
    C3 = float(C * C2)
    C4 = float(C2 * C2)
    D2 = float(D * D)
    D3 = float(D * D2)
    D4 = float(D2 * D2)
    E2 = float(E * E)
    E3 = float(E * E2)
    E4 = float(E2 * E2)
    Delta = float(256 * A3 * E3 - 192 * A2 * B * D * E2
                  - 128 * A2 * C2 * E2 + 144 * A2 * C * D2 * E
                  - 27 * A2 * D4 + 144 * A * B2 * C * E2
                  - 6 * A * B2 * D2 * E - 80 * A * B * C2 * D * E
                  + 18 * A * B * C * D3 + 16 * A * C4 * E
                  - 4 * A * C3 * D2 - 27 * B4 * E2 + 18 * B3 * C * D * E
                  - 4 * B3 * D3 - 4 * B2 * C3 * E + B2 * C2 * D2)
    Delta0 = float(C2 - 3 * B * D + 12 * A * E)
    Delta1 = float(2 * C3 - 9 * B * C * D
                   + 27 * B2 * E + 27 * A * D2 - 72 * A * C * E)
    p = float((8 * A * C - 3 * B2) / (8 * A2))
    q = float((B3 - 4 * A * B * C + 8 * A2 * D) / (8 * A3))
    R = cmath.sqrt(Delta1 * Delta1 - 4 * Delta0 * Delta0 * Delta0)
    P = cmath.sqrt(-27 * Delta)
    Q = complex(0.5 * (Delta1 + P)) ** complex(1.0 / 3.0)
    S = 0.5 * cmath.sqrt(((-2.0 * p) / 3.0)
                         + (Q + Delta0 / Q) / (3 * A))
    a = complex(-B / (4 * A), 0)
    s = 0.5 * cmath.sqrt(-4 * S * S - 2 * p + q / S)
    z1 = a - S + s
    z2 = a - S - s
    s = 0.5 * cmath.sqrt(-4 * S * S - 2 * p - q / S)
    z3 = a + S + s
    z4 = a + S - s

    print("z1 = ", z1)
    print("z2 = ", z2)
    print("z3 = ", z3)
    print("z4 = ", z4)

option = int(1)

print("Menu")
print("1 Solve a linear equation")
print("2 Solve a quadratic eqution")
print("3 Solve a cubic equation")
print("4 Solve a quartic equation")
print("5 Quit")

option = int(input("Option 1 to 5: "))

while option != 5:
    if option >= 1:
        a = float(input("a = "))
        b = float(input("b = "))
    if option >= 2:
        c = float(input("c = "))
    if option >= 3:
        d = float(input("d = "))
    if option == 4:
        e = float(input("e = "))

    if option == 1:
        z1 = -b / a
        print("z1 = ", z1)

    elif option == 2:
        real_quadratic_formula(a, b, c)

    elif option == 3:
        cubic_equation(a, b, c, d)

    elif option == 4:
        quartic_equation(a, b, c, d, e)

    option = int(input("Option 1 to 5: "))
Root Finder Output for Degrees 1, 2, and 3

Root Finder Output for Degrees 3 and 4

Windows Mutual Exclusion and Threading Example in C++ by James Pate Williams, Jr.

// James Pate Williams, Jr. (c) 10/04/2022
// Second version, no longer dependent on
// a global object
#include <windows.h>
#include <synchapi.h>
#include <chrono>
#include <iostream>
#include <vector>
using namespace std::chrono;
using namespace std;

// https://learn.microsoft.com/en-us/windows/win32/procthread/creating-threads
// https://www.geeksforgeeks.org/measure-execution-time-function-cpp/

long long counter;

class Test
{
public:
	inline void CountUp(
		HANDLE mutex,
		long long number)
	{
		counter = 0;

		for (long long i = 0; i < number; i++)
		{
			counter++;
		}

		ReleaseMutex(mutex);
	}
};

typedef struct Data
{
	long long number;
	DWORD dwWaitMS;
	Test* test;
} DATA, *PDATA;

DWORD WINAPI ThreadCode(LPVOID lpParam)
{
	Data* pData = (Data*)lpParam;
	HANDLE mutex = CreateMutex(NULL, FALSE, NULL);

	if (mutex != NULL)
	{
		pData->test->CountUp(
			mutex, pData->number);

		DWORD err = WaitForSingleObject(mutex, pData->dwWaitMS);
		CloseHandle(mutex);
	}

	return 0;
}

int main()
{
	long long number = 1000000000000;

	for (DWORD dwWaitMS = 30000;
		dwWaitMS <= 240000; dwWaitMS += 30000)
	{
		cout << "Number  = " << number << endl;
		cout << "Wait MS = " << dwWaitMS << endl;

		DWORD dwThreadId;
		Data* pData = new Data;

		pData->number = number;
		pData->dwWaitMS = dwWaitMS;
		pData->test = new Test();
		counter = 0;
		
		HANDLE hThread = NULL;

		hThread = CreateThread(
					NULL,			    // default security attributes
					0,					// use default stack size  
					ThreadCode,			// thread function name
					pData,				// argument to thread function 
					0,			        // use default creation flags 
					& dwThreadId);		// returns the thread identifier

		if (hThread != NULL)
		{
			auto start = high_resolution_clock::now();
			WaitForSingleObject(hThread, dwWaitMS);
			auto stop = high_resolution_clock::now();
			auto duration = duration_cast<milliseconds>(stop - start);
			double runtime = (double)duration.count() / 1.0e3;

			cout << "Counter = " << counter << endl;
			cout << "Runtime = " << runtime << endl;
			
			CloseHandle(hThread);
		}
	}
	return 0;
}
Number  = 1000000000000
Wait MS = 30000
Counter = 19226393400
Runtime = 30.006
Number  = 1000000000000
Wait MS = 60000
Counter = 7248012054
Runtime = 60.003
Number  = 1000000000000
Wait MS = 90000
Counter = 5264969696
Runtime = 90.011
Number  = 1000000000000
Wait MS = 120000
Counter = 10312016711
Runtime = 120
Number  = 1000000000000
Wait MS = 150000
Counter = 15342215540
Runtime = 150.011
Number  = 1000000000000
Wait MS = 180000
Counter = 20399182381
Runtime = 180
Number  = 1000000000000
Wait MS = 210000
Counter = 25358449910
Runtime = 210.009
Number  = 1000000000000
Wait MS = 240000
Counter = 30917682739
Runtime = 240.017

C:\Users\james\source\repos\MutexTest\Debug\MutexTest.exe (process 49516) exited with code 0.
Press any key to close this window . . .

Simulated Annealing and Other Algorithms to Optimally or Non-Optimally Solve the Traveling Salesperson Problem (TSP) by James Pate Williams, Jr.

In the Fall Semester of 2000 I took an industrial engineering course at Auburn University. The course was Industrial Systems 6700 and was taught by Professor and Chair Alice E. Smith. The course was entitled “Search Methods for Optimization”. The first two papers that Professor Smith handed out dealt with physical and computerized simulating annealing and were: “Equation of State Calculations by Fast Computing Machines” by Nicholas MetropolisArianna W. RosenbluthMarshall N. Rosenbluth, and Augusta H. Teller Los Alamos Scientific Laboratory, Los Alamos, New Mexico Edward Teller Notice that the inventor of the hydrogen bomb Edward Teller was involved with the paper and Optimization by Simulated Annealing | Science. The first paper dates from the year I was born namely 1953 and the second paper 1983, the latter year was the year that I flunked out of the Georgia Institute of Technology.

Our first assignment dealt with optimization via computerized simulating annealing (SA) and Professor Smith gave the class her implementation of the algorithm. By the way, over the last 22 years I have reused Smith’s algorithm in many of my own computer applications (C, C++, and C# computer languages).

SA gained the favor of Federal Express in 1983 and the Science paper mentioned above details the delivery service corporations use of the algorithm. The state of the art in solving the traveling salesperson problem (TSP) is the Concorde program written by scholars at the University of Waterloo in Canada: The traveling salesman problem | Mathematics | University of Waterloo (uwaterloo.ca). I believe Concorde uses linear programming.

With my software I am only able to solve small in size instances of the TSP. My implementations in a single computer application to solve toy random instances of the TSP for 3 to 25 number of locations are: brute force (3 to 9 cities), steepest descent local search (SDLS) from http://www.cs.bham.ac.uk/~wbl/biblio/gecco2007/docs/p1226.pdf, tabu search, evolutionary hill climber (EHC), genetic algorithm (GA), simulated annealing, some hybrid methods, minimum spanning tree, etc.

Steepest Descent Local Search

Tabu (Taboo) Search

Evolutionary Hill Climber

Genetic Algorithm

Simulated Annealing (Probably Optimal)

Back to Trial Division Using the C Programming Language by James Pate Williams, Jr.

I used Henri Cohen’s trial division algorithm from A Course in Computational Algebraic Number Theory. Trial division is the most primitive factoring method. I used a prime number sieve of Eratosthenes of all prime numbers less than or equal 100,000,000. There are 50,847,534 such primes. Also I used the C data type long long which yields a number of 64-bits. The vanilla C source is shown below as PDF file along with a run of the application. [Note: Somehow I lost this project, but fortunately I was working on an application that uses Brent’s modification of Pollard’s rho method. I have included the source code from that project below. The original code had a bug ridden Miller-Rabin probabilistic prime number test.]

Brent’s Method to Factor 2^63-3

Trial division is still a somewhat useful means of factoring large integers on modern massively parallel super computers. You can partition the factor base among the processors and run the algorithm in parallel.

p – 1 Method Missed Factor of 2

Trial Division

Sorting C++ Source Code and Testing Application by James Pate Williams, Jr.

I read the treatise Sorting and Sort Systems by Harold Lorin in the summer of 1979. That excellent intellectual work inspired me to implement the code and and algorithms from that book. Lorin used Algol, FORTRAN and IBM’s language PL/I. I used the Data General Eclipse minicomputer at LaGrange College from which I had graduated with a Bachelor of Arts degree in Chemistry. I utilized the computer language BASIC (Beginner’s All-purpose Symbolic Instruction Code) and Data General’s Advanced Operating System Macro-assembly language. Later in my academic life at Auburn University I learned about Heap Sort. I wrote about six sorting algorithms on Facebook that I implemented in the C# computer language back in 2015. It took me two days to create a similar application using the C++ computer language. Here are some results and the source code in PDF file format.

Insertion Sort Results
Selection Sort Results
Quick Sort Results
Heap Sort Results
Singleton’s Sort Results

Knight’s Tour C++ Application by James Pate Williams, Jr.

Using the chronological backtracking C++ code that I developed from the backtracking algorithm applied to the N-Queens constraint satisfaction problem, I developed an algorithm to solve the Knight’s Tour Problem. The general backtracking algorithm is from Foundations of Constraint Satisfaction by E. P. K. Tsang Chapter 2 page 37. I display some executions of the program in this blog.

Chronological Backtracking Closed Tour Solution for a 6×6 Chess Board
Warnsdorff Closed Tour Solution for a 6×6 Chess Board
Warnsdorff Open Tour Solution for an 8×8 Chess Board
Knight Move Count Matrix Used by Warnsdorff Heuristic

Data from another of my Knight’s Tour applications.

Eigenfunctions of a Nonlinear Second Order Ordinary Differential Equation Initial Value Eigenvalue Problem by James Pate Williams, Jr.

Solutions of a nonlinear second order ordinary differential equation initial value eigenvalue problem:

y”(x) = x * x * y(x) + n * n * y(x) * y(x) for all x in [0, 1]
y(0) = 0
y'(0) = 1
n in [0, 1, 2, …]

I graphed the first five eigenfunctions.

Graphs for n = 0 to n = 4

Planetary Precession by James Pate Williams, Jr.

There are three classic theoretical tests of Albert Einstein’s  Theory of General Relativity: the perihelion precession of Mercury, the other Solar System planets, and the planetoid Pluto, the bending of light by massive bodies, and the gravitational red shift. I recently wrote a C# program for displaying the exaggerated Rosette motion of theoretical planets (Schwarzschild’s solution to Einstein’s general relativity field equation that admit the existence of black holes). I also wrote a C++ program to calculate planetary precession values that agree with experimental results.

Precession.cpp (c) James Pate Williams, Jr. August 2022

This program calculates the planetary precessions of the planets in our solar system. Some of the  equations and data are from “General Relativity” by Hans Stephani 1982 page 103 and the following websites. Also, two calculations of the mass of the Sun are exhibited, along with my weight on different planets:

https://nssdc.gsfc.nasa.gov/planetary/factsheet/

https://farside.ph.utexas.edu/teaching/celestial/Celestialhtml/node44.html

https://imagine.gsfc.nasa.gov/features/yba/CygX1_mass/gravity/sun_mass.html

https://socratic.org/questions/how-do-you-calculate-the-mass-of-the-sun-m-sun-using-kepler-s-third-law-t-2-4-pi

https://en.wikipedia.org/wiki/Surface_gravity

https://www.schoolsobservatory.org/discover/quick/weight/https://physicscalc.com/physics/escape-velocity-calculator/#:~:text=Steps%20to%20Find%20Escape%20Velocity%201%20Obtain%20the,the%20double%20the%20result%20is%20the%20escape%20velocity

Planetary Precession C++ Program Output

Rosette Motion 1

Rosette Motion 2
Rosette Motion 3
Rosette Motion 4
Rosette Motion 5

Rosette Motion 6

Rosette Motion 7