Blog Entry © Saturday, November 28, 2025, by James Pate Williams, Jr. Fair Coin and Fair Dice Probabilities

Blog Entry © Tuesday 25, 2025, by James Pate Williams, Jr. Enumeration of the Tic-Tac-Toe Game Tree

/*
 *  Board.h
 *  Permutation
 *
 *  Created by James Pate Williams, Jr. on 7/21/08.
 *	Modified on Monday 11/24/2025 
 *  Copyright 2008 James Pate Williams, Jr. All rights reserved.
 *
 */

#ifndef _Board_
#define _Board_

#include <fstream>
#include <vector>

class Board {
	private:
		
		char state[9];
		std::vector<std::vector<char>> board;

	public:
		
		Board(void);
		Board(const Board &board);
		Board(const char b[3][3]);
		Board(const std::vector<char> &v);
		
		char GetState(int n) const;
		void PrintBoard(int n) const;
		void WriteBoard(std::ofstream &out) const;
		
		friend bool operator == (const Board &b1, const Board &b2);
		friend bool operator != (const Board &b1, const Board &b2);
		friend bool Found(const Board &initial, const Board &board);
		friend bool Win(const Board b, char winner[8]);
		friend bool LegalBoard(
			bool& legal, bool& won, char& whoWon, const Board b);
		friend void PrintBoard(char board[3][3], int number);
		
		static const char CharX, CharO, CharS;
};

#endif

/*
 *  Board.cpp
 *  Permutation
 *
 *  Created by James Pate Williams, Jr. on 7/21/08.
 *  Copyright 2008 James Pate Williams, Jr. All rights reserved.
 *
 */

#include <algorithm>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <vector>
#include "Board.h"

const char Board::CharO = 'o';
const char Board::CharX = 'x';
const char Board::CharS = '_';

Board::Board(void) {
	board.resize(3, std::vector<char>(3));
	for (int i = 0; i < 9; i++)
		state[i] = '\0';
}

Board::Board(const Board &b) {
	board.resize(3, std::vector<char>(3));
	for (int i = 0; i < 3; i++)
		for (int j = 0; j < 3; j++)
			state[3 * i + j] = board[i][j] = b.state[3 * i + j];
}

Board::Board(const char b[3][3]) {
	board.resize(3, std::vector<char>(3));
	for (int i = 0; i < 3; i++)
		for (int j = 0; j < 3; j++)
			state[3 * i + j] = board[i][j] = b[i][j];
}

Board::Board(const std::vector<char> &v) {
	board.resize(3, std::vector<char>(3));
	for (int i = 0; i < 3; i++)
		for (int j = 0; j < 3; j++)
			state[3 * i + j] = board[i][j] = v[3LL * i + j];
}

char Board::GetState(int n) const {
	return state[n];
}

void Board::PrintBoard(int n) const  {
	std::cout << n << std::endl;
	
	for (int i = 0; i < 3; i++) {
	
		for (int j = 0; j < 3; j++)
			std::cout << state[3 * i + j];
				
		std::cout << std::endl;
	}
	
	std::cout << std::endl;
	std::cout << std::endl;
}

void Board::WriteBoard(std::ofstream &out) const {
	for (int i = 0; i < 9; i++) 
		out << state[i];
				
	out << std::endl;
}

bool operator == (const Board &b1, const Board &b2) {
	int i, equal = 0;
	
	for (i = 0; i < 9; i++)
		if (b1.state[i] == b2.state[i])
			equal++;

	return equal == 9;
}

bool operator != (const Board &b1, const Board &b2) {
	for (int i = 0; i < 9; i++)
		if (b1.state[i] == b2.state[i])
			return false;
	
	return true;
}

bool Found(const Board &initial, const Board &board) {

	// reflections
	
	if (initial.state[0] == board.state[6] && initial.state[1] == board.state[7] &&
		initial.state[2] == board.state[8] && initial.state[3] == board.state[3] &&
		initial.state[4] == board.state[4] && initial.state[5] == board.state[5] &&
		initial.state[6] == board.state[0] && initial.state[7] == board.state[1] &&
		initial.state[8] == board.state[2])
		return true;
		
	if (initial.state[0] == board.state[2] && initial.state[1] == board.state[1] &&
		initial.state[2] == board.state[0] && initial.state[3] == board.state[5] &&
		initial.state[4] == board.state[4] && initial.state[5] == board.state[3] &&
		initial.state[6] == board.state[8] && initial.state[7] == board.state[7] &&
		initial.state[8] == board.state[6])
		return true;
	
	if (initial.state[0] == board.state[0] && initial.state[1] == board.state[3] &&
		initial.state[2] == board.state[6] && initial.state[3] == board.state[1] &&
		initial.state[4] == board.state[4] && initial.state[5] == board.state[7] &&
		initial.state[6] == board.state[2] && initial.state[7] == board.state[5] &&
		initial.state[8] == board.state[8])
		return true;

	if (initial.state[0] == board.state[8] && initial.state[1] == board.state[5] &&
		initial.state[2] == board.state[2] && initial.state[3] == board.state[7] &&
		initial.state[4] == board.state[4] && initial.state[5] == board.state[1] &&
		initial.state[6] == board.state[6] && initial.state[7] == board.state[3] &&
		initial.state[8] == board.state[0])
		return true;

	// rotations
	
	if (initial.state[0] == board.state[6] && initial.state[1] == board.state[3] &&
		initial.state[2] == board.state[0] && initial.state[3] == board.state[7] &&
		initial.state[4] == board.state[4] && initial.state[5] == board.state[1] &&
		initial.state[6] == board.state[8] && initial.state[7] == board.state[5] &&
		initial.state[8] == board.state[2])
		return true;
			
	if (initial.state[0] == board.state[8] && initial.state[1] == board.state[7] &&
		initial.state[2] == board.state[6] && initial.state[3] == board.state[5] &&
		initial.state[4] == board.state[4] && initial.state[5] == board.state[3] &&
		initial.state[6] == board.state[2] && initial.state[7] == board.state[1] &&
		initial.state[8] == board.state[0])
		return true;
	
	if (initial.state[0] == board.state[2] && initial.state[1] == board.state[5] &&
		initial.state[2] == board.state[8] && initial.state[3] == board.state[1] &&
		initial.state[4] == board.state[4] && initial.state[5] == board.state[7] &&
		initial.state[6] == board.state[0] && initial.state[7] == board.state[3] &&
		initial.state[8] == board.state[6])
		return true;
			
	// equal

	return initial == board;
}

bool Win(const Board b, char winner[8]) {
	bool win = false;
	
	for (int i = 0; i < 8; i++)
		winner[i] = Board::CharS;

	if (b.board[0][0] == b.board[0][1] && b.board[0][0] == 
		b.board[0][2] && b.board[0][0] != Board::CharS) {
		winner[0] = b.board[0][0];
		win = true;
	}
	
	if (b.board[1][0] == b.board[1][1] && b.board[1][0] ==
		b.board[1][2] && b.board[1][0] != Board::CharS)  {
		winner[1] = b.board[1][0];
		win = true;
	}
	
	if (b.board[2][0] == b.board[2][1] && b.board[2][0] == 
		b.board[2][2] && b.board[2][0] != Board::CharS)  {
		winner[2] = b.board[2][0];
		win = true;
	}
	
	if (b.board[0][0] == b.board[1][0] && b.board[0][0] == 
		b.board[2][0] && b.board[0][0] != Board::CharS)  {
		winner[3] = b.board[0][0];
		win = true;
	}
	
	if (b.board[0][1] == b.board[1][1] && b.board[0][1] == 
		b.board[2][1] && b.board[0][1] != Board::CharS)  {
		winner[4] = b.board[0][1];
		win = true;
	}
	
	if (b.board[0][2] == b.board[1][2] && b.board[0][2] == 
		b.board[2][2] && b.board[0][2] != Board::CharS)  {
		winner[5] = b.board[0][2];
		win = true;
	}
	
	if (b.board[0][0] == b.board[1][1] && b.board[0][0] == 
		b.board[2][2] && b.board[0][0] != Board::CharS)  {
		winner[6] = b.board[0][0];
		win = true;
	}                
	
	if (b.board[0][2] == b.board[1][1] && b.board[0][2] == 
		b.board[2][0] && b.board[0][2] != Board::CharS)  {
		winner[7] = b.board[0][2];
		win = true;
	}

	return win;
}
	
void PrintBoard(char board[3][3], int number) {
	
	std::cout << "number = " << number << std::endl;
	std::cout << std::endl;

	for (int i = 0; i < 3; i++) {
	
		for (int j = 0; j < 3; j++)
			std::cout << board[i][j];
		
		std:: cout << std::endl;
	}

	std::cout << std::endl;
}

bool LegalBoard(
	bool &legal, bool &won, char &whoWon, const Board b) {
	bool winO, winX;
	char winner[8];
	int countO = 0, countS = 0, countX = 0;
	int i, j, oWins = 0, xWins = 0;
	
	legal = won = false;
	
	for (i = 0; i < 3; i++) {
	
		for (j = 0; j < 3; j++) {
			if (b.board[i][j] == Board::CharO)
				countO++;
			else if (b.board[i][j] == Board::CharX)
				countX++;
			else if (b.board[i][j] == Board::CharS)
				countS++;
		}
	}

	if ((countX >= 3 || countO >= 3) && Win(b, winner)) {
	
		for (i = 0; i < 8; i++) {
			if (winner[i] == Board::CharX)
				xWins++;
			else if (winner[i] == Board::CharO)
				oWins++;
		}
		
		winX = (xWins  >= 1 && oWins == 0) && 
			   (countX <= 5 && countO <= 4 && countX == countO + 1);
		
		winO = (xWins  == 0 && oWins >= 1) &&
			   (countO <= 5 && countX <= 4 && countO == countX + 1);
		
		legal = won = winX || winO;

		if (won && winX)
			whoWon = Board::CharX;
		else if (won && winO)
			whoWon = Board::CharO;
	}
	else
		legal = (countX <= 5 && countO <= 4 && countX == countO + 1) ||
			    (countO <= 5 && countX <= 4 && countO == countX + 1);
			  
	return legal;
}

#define _WriteBV_
#include <algorithm>

#ifdef _WriteBV_
#include <fstream>
#endif

#include <iomanip>
#include <iostream>
#include <vector>
#include "Board.h"

static int Compare(const void *vPtr1, const void *vPtr2) {
	char *cPtr1 = (char *) vPtr1;
	char *cPtr2 = (char *) vPtr2;

	return strcmp(cPtr1, cPtr2);
}

static void EnumerateTicTacToe(void) {
	bool legal, won;
	char b[3][3], whoWon;
	char tree[3][3][3][3][3][3][3][3][3];
	int i0, i1, i2, i3, i4, i5, i6, i7, i8;
	int legals = 0, number = 0, wCount = 0;
	int xCount = 0, oCount = 0;
	int bvSize, i, j;

	std::vector<Board> boardVector;
	
	for (i0 = 0; i0 < 3; i0++) {
	for (i1 = 0; i1 < 3; i1++) {
	for (i2 = 0; i2 < 3; i2++) {
	for (i3 = 0; i3 < 3; i3++) {
	for (i4 = 0; i4 < 3; i4++) {
	for (i5 = 0; i5 < 3; i5++) {
	for (i6 = 0; i6 < 3; i6++) {
	for (i7 = 0; i7 < 3; i7++) {
	for (i8 = 0; i8 < 3; i8++) {
		if (i0 == 0)
			b[0][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i0 == 1)
			b[0][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i0 == 2)
			b[0][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;
				
		if (i1 == 0)
			b[0][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i1 == 1)
			b[0][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i1 == 2)
			b[0][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;

		if (i2 == 0)
			b[0][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i2 == 1)
			b[0][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i2 == 2)
			b[0][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;

		if (i3 == 0)
			b[1][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i3 == 1)
			b[1][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i3 == 2)
			b[1][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;

		if (i4 == 0)
			b[1][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i4 == 1)
			b[1][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i4 == 2)
			b[1][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;
			
		if (i5 == 0)
			b[1][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i5 == 1)
			b[1][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i5 == 2)
			b[1][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;
			
		if (i6 == 0)
			b[2][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i6 == 1)
			b[2][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i6 == 2)
			b[2][0] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;
			
		if (i7 == 0)
			b[2][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i7 == 1)
			b[2][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i7 == 2)
			b[2][1] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;
			
		if (i8 == 0)
			b[2][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharX;
		else if (i8 == 1)
			b[2][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharO;
		else if (i8 == 2)
			b[2][2] = tree[i0][i1][i2][i3][i4][i5][i6][i7][i8] = Board::CharS;
		
		Board bb(b);

		if (LegalBoard(legal, won, whoWon, bb)) {
			Board board(b);

			legals++;
		
			if (won) {
				wCount++;
				
				if (whoWon == Board::CharX)
					xCount++;
				else if (whoWon == Board::CharO)
					oCount++;
			}
			
			boardVector.push_back(board);
		}
		
		number++;

	}}}}}}}}}

	for (i = 0; i < boardVector.size() - 1; i++)
		for (j = i; j < boardVector.size(); j++)
			if (Found(boardVector[i], boardVector[j]))
				boardVector.erase(boardVector.begin() + j);

	bvSize = (int)boardVector.size();
	std::cout << std::setw(4) << "legal total = " << legals << std::endl;
	std::cout << std::setw(4) << "xWins total = " << xCount << std::endl;
	std::cout << std::setw(4) << "oWins total = " << xCount << std::endl;
	std::cout << std::setw(4) << "oxWin total = " << xCount + oCount << std::endl;
	std::cout << std::setw(4) << "unique size = " << boardVector.size();
	std::cout << std::endl;
	
#ifdef _WriteBV_
 
	std::cout << std::endl;
	std::ofstream outFile("Tic TacToe Game Tree Enumeration.txt");
	std::cout << std::endl;
	
	for (i = 0; i < bvSize; i++) {
		boardVector[i].PrintBoard(i + 1);
		boardVector[i].WriteBoard(outFile);
	}

#endif
}

int main (int argc, char * const argv[]) {
	EnumerateTicTacToe();
    return 0;
}

Blog Entry © Monday, November 24, 2025, by James Pate Williams, Jr., A* Informed Search Application to Solve the 15-Tile Puzzle

Blog Entry © Sunday, November 23, 2025, by James Pate Williams, Jr. Modification of My A* Informed Search Solver for the 8-Tile Puzzle

Blog Entry © Monday, November 17, 2025, by James Pate Williams, Jr. An Elitist Evolutionary Hill Climber to Solve the 8-Tile Puzzle

Blog Entry © Monday, September 8, 2025, by James Pate Williams, Jr., Comparison of Two Applications to Find the Minima of Classical Objective Functions

#pragma once
#include "pch.h"

class Functions
{

public:

    // Ackley's function

    static double lbx1[3];
    static double ubx1[3];
    static double fBest1;
    static double xBest1[3];

    static double f1(
        int n, std::vector<double> z);

    // Beale's function

    static double lbx2[3];
    static double ubx2[3];
    static double fBest2;
    static double xBest2[3];

    static double f2(
        int n, std::vector<double> z);
    
    // Booth's function

    static double lbx3[3];
    static double ubx3[3];
    static double fBest3;
    static double xBest3[3];

    static double f3(
        int n, std::vector<double> z);

    // Rosenbrock function

    static double lbx4[3];
    static double ubx4[3];
    static double fBest4;
    static double xBest4[3];

    static double f4(
        int n, std::vector<double> z);

    // Holder table function

    static double lbx5[3];
    static double ubx5[3];
    static double fBest5;
    static double xBest5[3];

    static double f5(
        int n, std::vector<double> z);

    // McCormick function

    static double lbx6[3];
    static double ubx6[3];
    static double fBest6;
    static double xBest6[3];

    static double f6(
        int n, std::vector<double> z);

    static double rosenbrock(
        int n, std::vector<double> x,
        std::vector<double>& g);

    static double mccormick(
        int n, std::vector<double> z,
        std::vector<double>& g);
};

#include "pch.h"
#include "Functions.h"

// Ackley's function

double Functions::lbx1[3] = { 0, -5, -5 };
double Functions::ubx1[3] = { 0, +5, +5 };
double Functions::fBest1 = 0.0;
double Functions::xBest1[3] = { 0, 0.0, 0.0 };

double Functions::f1(int n, std::vector<double> z)
{
	double e = exp(1.0), pi = 4.0 * atan(1.0);
	double x = z[1], y = z[2], pi2 = 2.0 * pi;

	return -20 * exp(-0.2 * sqrt(0.5 * (x * x + y * y))) -
		exp(0.5 * (cos(pi2 * x) + cos(pi2 * y))) + e + 20;
}

// Beale's function

double Functions::lbx2[3] = { 0, -4.5, -4.5 };
double Functions::ubx2[3] = { 0, +4.5, +4.5 };
double Functions::fBest2 = 0.0;
double Functions::xBest2[3] = { 0, 3.0, 2.5 };

double Functions::f2(int n, std::vector<double> z)
{
	double x = z[1], y = z[2];

	return pow(1.5 - x + x * y, 2) + pow(2.25 - x + x * y * y, 2) +
		pow(2.625 - x + x * y * y * y, 2);
}

// Booth's function

double Functions::lbx3[3] = {0, -10, -10};
double Functions::ubx3[3] = {0, +10, +10};
double Functions::fBest3 = 0.0;
double Functions::xBest3[3] = {0, 1.0, 3.0};

double Functions::f3(int n, std::vector<double> z)
{
	double x = z[1], y = z[2];

	return pow(x + 2 * y - 7, 2) + pow(2 * x + y - 5, 2);
}

// Rosenbrock function

double Functions::lbx4[3] = {0, -100.0, -100.0};
double Functions::ubx4[3] = {0, +100.0, +100.0};
double Functions::fBest4 = 0.0;
double Functions::xBest4[3] = {0, 1.0, 1.0};

double Functions::f4(int n, std::vector<double> x)
{
	double temp = x[2] - x[1] * x[1];

	return temp * temp * 100.0 + (1.0 - x[1]) * (1.0 - x[1]);
}

// Holder table function

double Functions::lbx5[3] = {0, -10.0, -10.0};
double Functions::ubx5[3] = {0, +10.0, +10.0};
double Functions::fBest5 = -19.2085;
double Functions::xBest5[3] = {0, 8.05502, 9.66459};

double Functions::f5(int n, std::vector<double> z)
{
	double pi = 4.0 * atan(1.0);
	double x = z[1], y = z[2];

	return -fabs(sin(x) * cos(y) *
		exp(fabs(1.0 - sqrt(x * x + y * y) / pi)));
}

// McCormick function

double Functions::lbx6[3] = {0, -1.5, -3};
double Functions::ubx6[3] = {0, +4, +4};
double Functions::fBest6 = -1.9133;
double Functions::xBest6[3] = {0, -0.54719, -1.54719};

double Functions::f6(int n, std::vector<double> z)
{
	double x = z[1], y = z[2];

	return sin(x + y) + pow(x - y, 2) - 1.5 * x + 2.5 * y + 1;
}

double Functions::rosenbrock(
	int n, std::vector<double> x,
	std::vector<double>& g)
{
	double temp;

	temp = x[2] - x[1] * x[1];
	g[1] = (-temp * 400.0 + 2.0) * x[1] - 2.0;
	g[2] = temp * 200.0;
	return temp * temp * 100.0 + (1.0 - x[1]) * (1.0 - x[1]);
}

double Functions::mccormick(
	int n, std::vector<double> z,
	std::vector<double>& g)
{
	double x = z[1], y = z[2];
	double c = cos(x + y);
	double p = 2.0 * (x - y);

	g[1] = c + p - 1.5;
	g[2] = c - p + 2.5;

	return sin(x + y) + pow(x - y, 2) - 1.5 * x + 2.5 * y + 1;
}

The PRAXIS and FLEMIN C source code can be found in the handbook “A Numerical Library in C for Scientists and Engineers” (c) 1996 by H. T. Lau, Ph.D. I translated the C code to C++ using Standard Template Library vectors.

Blog Entry © Wednesday, September 3, 2025, by James Pate Williams, Jr. Solution of Exercise 1.11 in Modern Quantum Chemistry an Introduction to Advanced Electronic Structure Theory © 1996 by Attila Szabo and Neil S. Ostlund

// EigenVV2x2.cpp (c) Wednesday, September 3, 2025
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Solution to Exercise 1.11 in "Quantum Chemistry
// an Introduction to Advanced Electronic Structure
// Theory (c) 1996 by Attila Szabo and Neil S. Ostlund 

#include "pch.h"
#include "framework.h"
#include "EigenVV2x2.h"

#define MAX_LOADSTRING 100

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
std::wstring text;                              // output wide string

// Forward declarations of functions included in this code module:
ATOM                MyRegisterClass(HINSTANCE hInstance);
BOOL                InitInstance(HINSTANCE, int);
LRESULT CALLBACK    WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    About(HWND, UINT, WPARAM, LPARAM);

int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    UNREFERENCED_PARAMETER(hPrevInstance);
    UNREFERENCED_PARAMETER(lpCmdLine);

    // TODO: Place code here.

    // Initialize global strings
    LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
    LoadStringW(hInstance, IDC_EIGENVV2X2, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_EIGENVV2X2));

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, nullptr, 0, 0))
    {
        if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
        {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }

    return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
    WNDCLASSEXW wcex;

    wcex.cbSize = sizeof(WNDCLASSEX);

    wcex.style          = CS_HREDRAW | CS_VREDRAW;
    wcex.lpfnWndProc    = WndProc;
    wcex.cbClsExtra     = 0;
    wcex.cbWndExtra     = 0;
    wcex.hInstance      = hInstance;
    wcex.hIcon          = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_EIGENVV2X2));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_EIGENVV2X2);
    wcex.lpszClassName  = szWindowClass;
    wcex.hIconSm        = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

    return RegisterClassExW(&wcex);
}

//
//   FUNCTION: InitInstance(HINSTANCE, int)
//
//   PURPOSE: Saves instance handle and creates main window
//
//   COMMENTS:
//
//        In this function, we save the instance handle in a global variable and
//        create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   hInst = hInstance; // Store instance handle in our global variable

   HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

#define IDC_STATIC1         1000
#define IDC_STATIC2         1010
#define IDC_STATIC3         1020
#define IDC_STATIC4         1030

#define IDC_EDIT_A11        2000
#define IDC_EDIT_A12        2010
#define IDC_EDIT_A21        2020
#define IDC_EDIT_A22        2030
#define IDC_EDIT_MULTILINE  3000

#define IDC_BUTTON_COMPUTE  4000
#define IDC_BUTTON_CANCEL   4010

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    static HFONT hFont = { };
    static HWND hEditMultiline = { };
    static HWND hEditA11 = { };
    static HWND hEditA12 = { };
    static HWND hEditA21 = { };
    static HWND hEditA22 = { };

    switch (message)
    {
    case WM_CREATE:
    {
        hFont = CreateFont(16, 0, 0, 0, FW_NORMAL, FALSE, FALSE, FALSE,
            ANSI_CHARSET, OUT_DEFAULT_PRECIS, CLIP_DEFAULT_PRECIS,
            DEFAULT_QUALITY, DEFAULT_PITCH | FF_SWISS, L"Courier New Bold");

        CreateWindowEx(0, L"STATIC", L"A[1][1]:", WS_CHILD | WS_VISIBLE,
            10, 10, 80, 20, hWnd, (HMENU)IDC_STATIC1, hInst, NULL);
        CreateWindowEx(0, L"STATIC", L"A[1][2]:", WS_CHILD | WS_VISIBLE,
            10, 40, 80, 20, hWnd, (HMENU)IDC_STATIC2, hInst, NULL);
        CreateWindowEx(0, L"STATIC", L"A[2][1]:", WS_CHILD | WS_VISIBLE,
            10, 70, 80, 20, hWnd, (HMENU)IDC_STATIC3, hInst, NULL);
        CreateWindowEx(0, L"STATIC", L"A[2][2]:", WS_CHILD | WS_VISIBLE,
            10, 100, 80, 20, hWnd, (HMENU)IDC_STATIC4, hInst, NULL);

        hEditA11 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
            100, 10, 200, 20, hWnd, (HMENU)IDC_EDIT_A11, hInst, NULL);
        hEditA12 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
            100, 40, 200, 20, hWnd, (HMENU)IDC_EDIT_A12, hInst, NULL);
        hEditA21 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
            100, 70, 200, 20, hWnd, (HMENU)IDC_EDIT_A21, hInst, NULL);
        hEditA22 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
            100, 100, 200, 20, hWnd, (HMENU)IDC_EDIT_A22, hInst, NULL);

        hEditMultiline = CreateWindowEx(
            WS_EX_CLIENTEDGE,                       // Extended style for sunken border
            TEXT("EDIT"),                           // Class name
            TEXT(""),                               // Initial text (can be blank)
            WS_CHILD | WS_VISIBLE | WS_VSCROLL | ES_AUTOHSCROLL |
            ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL | WS_HSCROLL | WS_VSCROLL,
            310, 10, 300, 300,                      // Position and size
            hWnd,                                   // Parent window handle
            (HMENU)IDC_EDIT_MULTILINE,              // Unique control ID
            hInst,                                  // Application instance
            NULL                                    // Extra parameter
        );

        CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            10, 130, 80, 30, hWnd, (HMENU)IDC_BUTTON_COMPUTE, hInst, NULL);
        CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            220, 130, 80, 30, hWnd, (HMENU)IDC_BUTTON_CANCEL, hInst, NULL);

        SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, TRUE);
        SetDlgItemText(hWnd, IDC_EDIT_A11, L"3");
        SetDlgItemText(hWnd, IDC_EDIT_A12, L"1");
        SetDlgItemText(hWnd, IDC_EDIT_A21, L"1");
        SetDlgItemText(hWnd, IDC_EDIT_A22, L"3");
    }
    break;
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);
            // Parse the menu selections:
            switch (wmId)
            {
            case IDM_ABOUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
                break;
            case IDC_BUTTON_COMPUTE:
            {
                WCHAR line[128] = L"";

                text = L"";

                // eigenvalues
                std::vector<double> omega(3);
                // matrix
                std::vector<std::vector<double>> A(3,
                    std::vector<double>(3));
                std::vector<std::vector<double>> B(3,
                    std::vector<double>(3));
                // eigenvectors
                std::vector<std::vector<double>> c(3,
                    std::vector<double>(3));

                GetWindowText(hEditA11, line, 128);
                std::wstring A11Str(line);
                A[1][1] = std::stod(A11Str);
                GetWindowText(hEditA12, line, 128);
                std::wstring A12Str(line);
                A[1][2] = std::stod(A12Str);
                GetWindowText(hEditA21, line, 128);
                std::wstring A21Str(line);
                A[2][1] = std::stod(A21Str);
                GetWindowText(hEditA22, line, 128);
                std::wstring A22Str(line);
                A[2][2] = std::stod(A22Str);
                
                double term1 = A[1][1] + A[2][2];
                double term2 = pow(A[2][2] - A[1][1], 2.0);
                double term3 = 4.0 * A[1][2] * A[2][1];

                // compute eigenvalues

                omega[1] = 0.5 * (term1 - sqrt(term2 + term3));
                omega[2] = 0.5 * (term1 + sqrt(term2 + term3));

                swprintf_s(line, L"Eigenvalues:\r\n\r\n");
                text += std::wstring(line);
                swprintf_s(line, L"omega 1 = %13.10lf\r\n", omega[1]);
                text += std::wstring(line);
                swprintf_s(line, L"omega 2 = %13.10lf\r\n\r\n", omega[2]);
                text += std::wstring(line);

                // compute eigenvalues using a unitary transformation 
                // matrix A must be symmetric

                double theta0 = 0.0;

                if (A[1][1] == A[2][2])
                {
                    theta0 = 0.5 * acos(0.0);
                }

                else
                {
                    theta0 = 0.5 * atan(2.0 * A[1][2] /
                        (A[1][1] - A[2][2]));
                }

                // compute the eigenvalues
                
                omega[1] = A[1][1] * pow(cos(theta0), 2.0) +
                    A[2][2] * pow(sin(theta0), 2.0) +
                    A[1][2] * sin(2.0 * theta0);
                omega[2] = A[1][1] * pow(cos(theta0), 2.0) +
                    A[2][2] * pow(sin(theta0), 2.0) -
                    A[1][2] * sin(2.0 * theta0);
                
                swprintf_s(line, L"Eigenvalues:\r\n\r\n");
                text += std::wstring(line);
                swprintf_s(line, L"omega 1 = %13.10lf\r\n", omega[1]);
                text += std::wstring(line);
                swprintf_s(line, L"omega 2 = %13.10lf\r\n\r\n", omega[2]);
                text += std::wstring(line);

                // compute eigenvectors

                c[1][1] = cos(theta0);
                c[1][2] = sin(theta0);
                c[2][1] = sin(theta0);
                c[2][2] = -cos(theta0);

                swprintf_s(line, L"Eigenvectors:\r\n\r\n");
                text += std::wstring(line);
                swprintf_s(line, L"c 11 = %13.10lf\r\n", c[1][1]);
                text += std::wstring(line);
                swprintf_s(line, L"c 12 = %13.10lf\r\n", c[1][2]);
                text += std::wstring(line);
                swprintf_s(line, L"c 21 = %13.10lf\r\n", c[2][1]);
                text += std::wstring(line);
                swprintf_s(line, L"c 22 = %13.10lf\r\n", c[2][2]);
                text += std::wstring(line);
                SetWindowText(hEditMultiline, text.c_str());
                break;
            }
            case IDC_BUTTON_CANCEL:
            {
                PostQuitMessage(0);
                break;
            }
            break;
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
        {
            PAINTSTRUCT ps;
            HDC hdc = BeginPaint(hWnd, &ps);
            // TODO: Add any drawing code that uses hdc here...
            EndPaint(hWnd, &ps);
        }
        break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

Blog Entry (c) Wednesday, August 13, 2025, by James Pate Williams, Jr. Exercises from an Online Textbook

#include <complex>
#include <vector>

class CmpLinearAlgebra
{
public:
    static void CmpPrintMatrix(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& Ac);
    static void CmpAddition(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpAnticommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpCommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpAdjoint(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& Ac,
        std::vector<std::vector<std::complex<double>>>& Ad);
    static std::complex<double> CmpDeterminant(
        bool& failure, int n,
        std::vector<std::vector<std::complex<double>>>& A);
    static bool CmpGaussianElimination(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpGaussianFactor(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<size_t>& pivot);
    static bool CmpGaussianSolution(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpSubstitution(
        int m, int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpInverse(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& Mi);
    static void CmpCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<std::complex<double>>>& C,
        std::vector<std::vector<std::complex<double>>>& I,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& adjoint,
        std::vector<std::complex<double>>& a);
    static void CmpMatrixKernel(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& X,
        size_t& r);
    static void CmpMatrixImage(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& N,
        std::vector<std::vector<std::complex<double>>>& X,
        int rank);
};
#include <vector>

class DblLinearAlgebra
{
public:
    static void DblPrintMatrix(
        int m, int n, std::vector<std::vector<double>>& A);
    static void DblAddition(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblAnticommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblCommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static double DblDeterminant(
        bool& failure, int n,
        std::vector<std::vector<double>>& A);
    static bool DblGaussianElimination(
        int m, int n, std::vector<std::vector<double>>& A,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblGaussianFactor(
        int n, std::vector<std::vector<double>>& M,
        std::vector<size_t>& pivot);
    static bool DblGaussianSolution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblSubstitution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblInverse(
        int n, std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& A);
    static void DblCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<double>>& C,
        std::vector<std::vector<double>>& I,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& adjoint,
        std::vector<double>& a);
    static void DblMatrixKernel(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& X,
        size_t& r);
    static void DblMatrixImage(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& N,
        std::vector<std::vector<double>>& X,
        int rank);
};
// Exercises from "Modern Quantum Chemistry An Introduction to Advanced
// Electronic Structure Theory" by Attila Szabo and Neil S. Ostlund
// https://chemistlibrary.wordpress.com/wp-content/uploads/2015/02/modern-quantum-chemistry.pdf
// Program (c) Tuesday, August 12, 2025 by James Pate Williams, Jr.

#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"

int main()
{
	double AArcb[3][3] = { { 2, 3, -1 }, { 4, 4, -3 }, { -2, 3, -1 } };
	double AArso[3][3] = { { 1, 1, 0 }, { 1, 2, 2 }, { 0, 2, -1 } };
	double BBrso[3][3] = { { 1, -1, 1 }, { -1, 0, 0 }, { 1, 0, 1} };
	double BBr[3][3] = { { 1, -1, 1 }, { -1 , 0, 0 }, { 1, 0, 1 } };
	double AAcr[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	double AAci[3][3] = { { 1, 1, 2 }, { 3, 0, 1 }, { 0, 2, 4 } };
	double BBcr[3][3] = { { 1, 0, 1 }, { 1 , 1, 0 }, { 0, 1, 1 } };
	double BBci[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	int m = 3, n = 3, p = 3;

	std::vector<double> br(3);
	std::vector<size_t> pivot(3);
	std::vector<std::vector<double>> Arcb(3, std::vector<double>(3));
	std::vector<std::vector<double>> Arso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Brso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Br(3, std::vector<double>(3));
	std::vector<std::vector<double>> Cr(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ai(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ari(3, std::vector<double>(3));
	std::vector<std::vector<std::complex<double>>> Ac(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Bc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Cc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Dc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Ec(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Fc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Gc(3,
		std::vector<std::complex<double>>(3));
	
	for (int i = 0; i < m; i++)
	{
		for (int j = 0; j < p; j++)
		{
			Arcb[i][j] = AArcb[i][j];
			Arso[i][j] = AArso[i][j];
			Brso[i][j] = BBrso[i][j];
			Ac[i][j]._Val[0] = AAcr[i][j];
			Ac[i][j]._Val[1] = AAci[i][j];
		}
	}

	for (int i = 0; i < p; i++)
	{
		for (int j = 0; j < n; j++)
		{
			Br[i][j] = BBr[i][j];
			Bc[i][j]._Val[0] = BBcr[i][j];
			Bc[i][j]._Val[1] = BBci[i][j];
		}
	}
	
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Br, Cr);
	std::cout << "Ar * Br = Cr Conte & de Boor" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << "Ac * Bc = Cc" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
	std::cout << std::endl;
	// Exercise 1.2
	std::cout << "Exercise 1.2 page 5 Commutator" << std::endl;
	DblLinearAlgebra::DblCommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	std::cout << "Exercise 1.2 page 5 Anticommutator" << std::endl;
	DblLinearAlgebra::DblAnticommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Cc, Dc);
	std::cout << "Exercise 1.3 page 6 Cc adjoint" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Dc);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Ac, Ec);
	CmpLinearAlgebra::CmpAdjoint(3, 3, Bc, Fc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Fc, Ec, Gc);
	std::cout << "Exercise 1.3 page 6 Bc adjoint * Ac adjoint" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Gc);
	std::cout << std::endl;
	std::cout << "Ar matrix" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Arcb);
	bool inv = DblLinearAlgebra::DblInverse(n, Arcb, Ai);
	std::cout << std::endl;
	std::cout << "Ar Conte & de Boor inverse = " << inv << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ai);
	std::cout << std::endl;
	std::cout << "Ar * Ar inverse" << std::endl;
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Ai, Ari);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ari);
	std::cout << std::endl;
	std::cout << "Ac" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
	std::cout << std::endl;
	inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
	std::cout << "Ac inverse = " << inv << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Bc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << std::endl;
	std::cout << "Ac * AC inverse" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
}

Blog Entry © Sunday, July 27, 2025, A Bit of Programming Nostalgia Prime Number Related Programs by James Williams, Jr.

/*
  Author:  Pate Williams c 1995

  The following program is a solution to problem 18.15
  in Pascalgorithms by Edwin D. Reilly and Francis D.
  Federighi page 627. The program uses Simpson's rule
  to calculate the number of primes less than or equal
  a given number.
*/

#include <math.h>
#include <stdio.h>

typedef double real;

static real f(real x)
{
	return(1.0 / log(x));
}

static real simpson(int n, real a, real b)
{
	int i;
	real evensum, h, oddsum, twoh, x;

	if (n % 2 == 1) n = n - 1;
	h = (b - a) / n;
	twoh = h + h;
	x = h + a;
	oddsum = 0.0;
	for (i = 1; i <= n / 2; i++)
	{
		oddsum += f(x);
		x = twoh + x;
	}
	x = twoh + a;
	evensum = 0.0;
	for (i = 1; i <= n / 2 - 1; i++)
	{
		evensum += f(x);
		x = twoh + x;
	}
	return(h / 3.0 * (f(a) + f(b) + 4.0 * oddsum + 2.0 * evensum));
}

int main(void)
{
	int i, n, Nmaximum = 0, Nminimum = 0, Nstep;

	printf("n = ");			scanf_s("%d", &n);
	printf("N minimum = "); scanf_s("%d", &Nminimum);
	printf("N maximum = "); scanf_s("%d", &Nmaximum);
	printf("N step = ");	scanf_s("%d", &Nstep);
	printf("\n");
	printf("----------------------------------------\n");
	printf("Min\t\tMax\t\tprimes\n");
	printf("----------------------------------------\n");
	for (i = Nminimum; i <= Nmaximum; i += Nstep)
	{
		printf("%8d\t%8d\t%8.0lf\n", Nminimum, i + Nstep,
			simpson(n, Nminimum, i + Nstep));
	}
	printf("----------------------------------------\n");
	return(0);
}
n = 1024
N minimum = 0
N maximum = 10000000
N step = 1000000

----------------------------------------
Min             Max             primes
----------------------------------------
       0         1000000           78551
       0         2000000          148923
       0         3000000          216788
       0         4000000          283122
       0         5000000          348361
       0         6000000          412754
       0         7000000          476461
       0         8000000          539590
       0         9000000          602224
       0        10000000          664424
       0        11000000          726239
----------------------------------------

D:\PrimeCounter\x64\Release\PrimeCounter.exe (process 51884) exited with code 0 (0x0).
Press any key to close this window . . .
/*
  Author:  Pate Williams c 1995

  The following is a translation of the Pascal program
  sieve found in Pascalgorithms by Edwin D. Reilly and
  Francis D. Federighi page 652. This program uses sets
  to represent the sieve (see C Programming Language An
  Applied Perspective by Lawrence Miller and Alec Qui-
  lici pages 160 - 162).
*/

#include <math.h>
#include <stdio.h>

#define _WORD_SIZE 32
#define _VECT_SIZE 524288
#define SET_MIN    0
#define SET_MAX    16777215

typedef unsigned long SET[_VECT_SIZE];
typedef long ELEMENT;
typedef unsigned long LONG;

SET set;

static int get_bit_pos(int* long_ptr, int* bit_ptr,
	ELEMENT element)
{
	*long_ptr = element / _WORD_SIZE;
	*bit_ptr = element % _WORD_SIZE;
	return(element >= SET_MIN && element <= SET_MAX);
}

static void set_bit(ELEMENT element, int inset)
{
	int bit, word;

	if (get_bit_pos(&word, &bit, element))
	{
		if (inset > 0)
			set[word] |= (01 << bit);
		else
			set[word] &= ~(01 << bit);
	}
}

static int get_bit(ELEMENT element)
{
	int bit, word;

	return(get_bit_pos(&word, &bit, element) ?
		(set[word] >> bit) & 01 : 0);
}

static void set_Add(ELEMENT element)
{
	set_bit(element, 1);
}

static void set_Del(ELEMENT element)
{
	set_bit(element, 0);
}

static int set_Mem(ELEMENT element)
{
	return get_bit(element);
}

static void primes(long n)
{
	long c, i, inc, k;
	double x;

	set_Add(2);
	for (i = 3; i <= n; i++)
		if ((i + 1) % 2 == 0)
			set_Add(i);
		else
			set_Del(i);
	c = 3;
	do
	{
		i = c * c;
		inc = c + c;
		while (i <= n)
		{
			set_Del(i);
			i = i + inc;
		}
		c += 2;
		while (set_Mem(c) == 0) c += 1;
	} while (c * c <= n);
	k = 0;
	for (i = 2; i <= n; i++)
		if (set_Mem(i) == 1) k++;
	x = n / log(n) - 5.0;
	x = x + exp(1.0 + 0.15 * log(n) * sqrt(log(n)));
	printf("%8ld\t%8ld\t%8.0lf\n", n, k, x);
}

int main(void)
{
	long n = 100L;

	printf("----------------------------------------\n");
	printf("n\t\tprimes\t\ttheory\n");
	printf("----------------------------------------\n");
	do
	{
		primes((int)n);
		n = 10L * n;
	} while (n < (long)SET_MAX);
	printf("----------------------------------------\n");
	return(0);
}
----------------------------------------
n               primes          theory
----------------------------------------
     100              25              29
    1000             168             181
   10000            1229            1261
  100000            9592            9634
 1000000           78498           78396
10000000          664579          665060
----------------------------------------

D:\Sieve\x64\Release\Sieve.exe (process 60092) exited with code 0 (0x0).
Press any key to close this window . . .

Blog Entry (c) Tuesday, June 3, 2025, Sorting Algorithms by James Pate Williams, Jr.

First, we make sure the sorts actually work as expected and then we do some timing tests.

^^ Integer Data ^^
** Options Menu **
1 Single Sorting Tests
2 Statistical Tests
3 Create log file and log events
4 Exit
Option number: 1
-- Single Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 1
Insertion Sort
PRNG Seed 1
Number of Samples 20
Maximum Sample Value 20
11 -17
-16 -16
9 -14
-17 -12
-14 -7
17 0
10 0
0 3
7 6
3 7
8 7
0 8
-7 9
11 10
-12 11
13 11
7 12
12 13
16 16
6 17
-- Single Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 2
std::qsort
PRNG Seed 1
Number of Samples 20
Maximum Sample Value 20
11 -17
-16 -16
9 -14
-17 -12
-14 -7
17 0
10 0
0 3
7 6
3 7
8 7
0 8
-7 9
11 10
-12 11
13 11
7 12
12 13
16 16
6 17
-- Single Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 3
Singleton's FORTRAN Sort
PRNG Seed 1
Number of Samples 20
Maximum Sample Value 20
11 -17
-16 -16
9 -14
-17 -12
-14 -7
17 0
10 0
0 3
7 6
3 7
8 7
0 8
-7 9
11 10
-12 11
13 11
7 12
12 13
16 16
6 17
-- Single Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number:
^^ Integer Data ^^
** Options Menu **
1 Single Sorting Tests
2 Statistical Tests
3 Create log file and log events
4 Exit
Option number: 2
-- Statistical Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 1
Insertion Sort
PRNG Seed 1
Number of Samples 1000
Maximum Sample Value 1000
Number of Experiments 100
Runtimes in microseconds
Minimum runtime 524
Maximum runtime 1751
Mean runtime 802
Median runtime 668
-- Statistical Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 2
std::qsort
PRNG Seed 1
Number of Samples 1000
Maximum Sample Value 1000
Number of Experiments 100
Runtimes in microseconds
Minimum runtime 115
Maximum runtime 1751
Mean runtime 481
Median runtime 391
-- Statistical Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 3
Singleton's FORTRAN Sort
PRNG Seed 1
Number of Samples 1000
Maximum Sample Value 1000
Number of Experiments 100
Runtimes in microseconds
Minimum runtime 93
Maximum runtime 1751
Mean runtime 363
Median runtime 174
-- Statistical Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number:
; Copilot and James Pate Williams, Jr.
; 2/8/2025 - 2/9/2025
; We use the eax register for array indices
; The array base in register ecx
; The register ebx is general purpose;
;
;class SortingCPP {
;public:
;	static void InsertionSort(std::vector<T>& a)
;	{
;		for (size_t j = 1; j < a.size(); j++)
;		{
;			T key = a[j];
;			int i = j - 1;
;
;			while (i >= 0 && a[i] > key)
;			{
;				a[i + 1] = a[i];
;				i--;
;			}
;
;			a[i + 1] = key;
;		}
;	};

.MODEL FLAT, C
.STACK 4096

.DATA
    ; Allocate space for uninitialized variables
    i DWORD ?
    j DWORD ?
    key DWORD ?
    n DWORD ?
    t DwORD ?

.CODE
InsertionSortASM PROC

    ; Parameters:
    ; array = [esp + 8]
    ; n = [esp + 12]

    push ebp
    mov ebp, esp
    sub esp, 16                 ; Allocate space for local variables

    mov ecx, [ebp + 8]          ; base of array
    mov eax, [ebp + 12]         ; n number of array elements
    mov [n], eax                ; store n

    ; Initialize variables
    mov dword ptr [i], 0        ; i = 0
    mov dword ptr [j], 1        ; j = 1

for_loop:

    mov eax, [j]                ; load j into register
    mov ebx, [n]                ; load n
    cmp eax, ebx                ; compare j to n
    je  Exit                    ; we are done

    mov ebx, [ecx + eax * 4]    ; ebx = a[j]
    mov [key], ebx              ; key = a[j]
    dec eax                     ; j = j - 1 
    mov [i], eax;               ; i = j - 1
    inc eax                     ; increment
    inc eax                     ; j = j + 1
    mov [j], eax                ; store j

while_loop:

    mov eax, [i]                ; load i into register
    cmp eax, -1                 ; is i == -1 ?
    jz  end_while               ; end the while loop
    mov ebx, [ecx + eax * 4]    ; load a[i]
    mov eax, [key]              ; load key into register
    cmp ebx, eax                ; compare a[i] to key
    jle end_while               ; end the while loop
    mov eax, [i]                ; load i
    mov ebx, [ecx + eax * 4]    ; load a[i]
    inc eax                     ; eax = i + 1
    mov edx, [ecx + eax * 4]    ; load a[i + 1]
    ;mov [t], ebx               ; t = a[i]
    mov edx, ebx                ; edx = a[i]
    mov eax, [i]                ; load i again
    inc eax                     ; i + 1
    mov [ecx + eax * 4], edx    ; a[i + 1] = a[i]
    dec eax                     ; i--
    dec eax                     ; i--
    mov [i], eax                ; store updated i
    jmp while_loop              ; continue while

end_while:

    mov eax, [i]                ; load i
    inc eax                     ; eax = i + 1
    mov ebx, [key]              ; ebx = key
    mov [ecx + eax * 4], ebx    ; a[i + 1] = key

    jmp for_loop                ; continue for loop

Exit:

    mov esp, ebp
    pop ebp
    ret

InsertionSortASM ENDP
END