Tag: programming

Rice-Golomb Encoder and Decoder Copyright (c) Thursday, April 3, 2025, to Sunday, April 6, 2025, by James Pate Williams, Jr. BA, BS, Master of Software Engineering, Doctor of Philosophy Computer Science

Online references:

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

Click to access 19790014634.pdf

// Rice-Golomb Encoder and Decoder
// Copyright (c) Thursday, April 3, 2025
// by James Pate Williams, Jr.
// BA, BS, Master of Software Engineering
// Doctor of Philosophy Computer Science
// Online references:
// https://en.wikipedia.org/wiki/Golomb_coding
// https://ntrs.nasa.gov/api/citations/19790014634/downloads/19790014634.pdf

#include <iostream>
#include <string>
#include <vector>
//#include <stdlib.h>

bool Encode(const char* NChars, size_t NCharsCount,
    long M, long long& N, std::vector<char>& qBits,
    std::vector<char>& rBits, unsigned int& qSize, unsigned int& rSize,
    long long& q, long long& r, unsigned int& NSize) {
    N = NChars[0] - (long long)'0';
    for (unsigned int i = 1; i < NCharsCount; i++) {
        N = 10 * N + (long long)NChars[i] - (long long)'0';
    }
    q = N / M;
    r = N % M;
    qSize = 0;
    while (qSize < q) {
        qBits.push_back('1');
        qSize++;
    }
    qBits.push_back('0');
    qSize++;
    rSize = 0;
    unsigned int b = (unsigned int)floor(log2(M));
    if (b > 62) {
        return false;
    }
    long long p = (long long)pow(2, b + 1);
    if (r < p - M) {
        long long rr = r;
        while (rr > 0) {
            long long digit = (rr & 1) == 1 ? 1 : 0;
            rBits.push_back((char)digit + '0');
            rSize++;
            rr >>= 1;
        }
        rBits.push_back('0');
        rSize++;
    }
    else {
        long long rr = r + p - M;
        while (rSize < b + 1) {
            long long digit = rr & 1 ? 1 : 0;
            rBits.push_back((char)digit + '0');
            rSize++;
            rr >>= 1;
        }
    }
    long long rValue = rBits[0];
    for (size_t i = 1; i < rSize; i++) {
        rValue = rValue * 2 + rBits[i];
    }
    long long NBitCount = 0;
    while (N > 0) {
        N >>= 1;
        NBitCount++;
    }
    std::cout << "q-bits size = " << qSize << std::endl;
    std::cout << "r-bits size = " << rSize << std::endl;
    std::cout << "N-bits size = " << qSize + rSize << std::endl;
    std::cout << "N-Chars * 8-Bits per Char = " << NCharsCount * 8 << std::endl;
    std::cout << "% Compression = " << 100.0 * (1.0 - (qSize + rSize) /
        (NCharsCount * 8.0)) << std::endl;
    return true;
}

void Decode(long long M, long long& N,
    std::vector<char> qBits, std::vector<char> rBits,
    unsigned int& qSize, unsigned int& rSize,
    long long& q, long long& r) {
    int count = 0;
    while (qBits[count] != '0') {
        count++;
    }
    q = count;
    int c = (int)rSize - 1;
    unsigned int b = (unsigned int)floor(log2(M));
    long long p = (long long)pow(2, b + 1);
    long long s = 0;
    r = rBits[c--] - (long long)'0';
    do {
        r = 2 * r + rBits[c] - (long long)'0';
        c--;
    } while (c >= 0);
    if (r < p - M) {
        s = r;
    }
    else {
        s = r + p - M;
        c = 1;
        r = rBits[0] - (long long)'0';
        while (c < (int)(b + 1)) {
            r = 2 * r + rBits[c] - (long long)'0';
            c++;
        }
        s = r;
    }
    r = s;
    N = q * M + r;
}

int main() {
    char line[128] = { };
    size_t NSize = 0, qSize = 0, rSize = 0;
    long long M = 10, N = 42, q = -1, r = -1;
    std::vector<char> qBits, rBits;
    std::cout << "M = ";
    std::cin.getline(line, 127);
    std::string str1(line);
    M = std::stoi(str1);
    std::cout << "N = ";
    std::cin.getline(line, 127);
    std::string str2(line);
    Encode(str2.c_str(), strlen(str2.c_str()), M, N,
        qBits, rBits, qSize, rSize, q, r, NSize);
    std::cout << "q = " << q << std::endl;
    std::cout << "r = " << r << std::endl;
    std::cout << "q-size = " << qSize << std::endl;
    std::cout << "r-size = " << rSize << std::endl;
    std::cout << "q ";
    for (unsigned int i = 0; i < qSize; i++) {
        std::cout << qBits[i] << ' ';
    }
    std::cout << std::endl;
    std::cout << "r ";
    for (int i = (int)rSize - 1; i >= 0; i--) {
        std::cout << rBits[i] << ' ';
    }
    std::cout << std::endl;
    Decode(M, N, qBits, rBits, qSize, rSize, q, r);
    std::cout << "q = " << q << std::endl;
    std::cout << "r = " << r << std::endl;
    std::cout << "q-size = " << qSize << std::endl;
    std::cout << "r-size = " << rSize << std::endl;
    std::cout << "q ";
    for (unsigned int i = 0; i < qSize; i++) {
        std::cout << qBits[i] << ' ';
    }
    std::cout << std::endl;
    std::cout << "r ";
    for (int i = rSize - 1; i >= 0; i--) {
        std::cout << rBits[i] << ' ';
    }
    std::cout << std::endl;
    std::cout << "N = " << N << std::endl;
    return 0;
}

M = 64
N = 1027
q-bits size = 17
r-bits size = 3
N-bits size = 20
N-Chars * 8-Bits per Char = 32
% Compression = 37.5
q = 16
r = 3
q-size = 17
r-size = 3
q 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
r 0 1 1
q = 16
r = 3
q-size = 17
r-size = 3
q 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
r 0 1 1
N = 1027

Chapter One Straight-Line Program Interpreter from “Modern Compiler Implementation in Java Second Edition” (c) 2002 by Andrew W. Appel, Translation to C++ by James Pate Williams, Jr. on Thursday, April 3, 2025

Wayback in the Spring Semester of 2006, after I was awarded my Doctor of Philosophy Degree in Computer Science, I partially audited a Compiler Design Course. Due to my mental aberrations, I was unable to complete the course. The instructor was on a Sabbatical from the United States Air Force Academy in Colorado Springs, Colorado. The textbook we used, and I still have a copy, was “Modern Compiler Implementation in Java Second Edition” © 2002 by Andrew W. Appel. Below is a translation from Java to C++ that I just completed.

// Chapter One program translated from Java to C++ by
// James Pate Williams, Jr. (c) Wednesday April 3, 2025
// Reference: "Modern Complier Implementation in Java
// Second Edition" (c) 2002 by Andrew W. Appel

#ifndef _SLPInterpreter_H
#include <iostream>
#include <stack>
#include <string>
#include <vector>

class TableEntry {
public:
	std::string symbol, value;
	TableEntry(std::string symbol, std::string value) {
		this->symbol = symbol;
		this->value = value;
	}
};

std::stack<std::string> sStack;
std::vector<TableEntry> symbolTable;

class Exp {
public:
	Exp() { };
	virtual ~Exp() { };
};

std::stack<Exp> eStack;

class ExpList {
public:
	ExpList() { };
	virtual ~ExpList() { };
};

class Stm {
public:
	Stm() { };
	virtual ~Stm() { };
};

class CompoundStm : public Stm {
public:
	Stm stm1, stm2;
	CompoundStm(Stm stm1, Stm stm2) {
		this->stm1 = stm1;
		this->stm2 = stm2;
	};
};

class AssignStm : public Stm {
public:
	std::string id;
	Exp exp;
	AssignStm(std::string id, Exp exp) {
		this->id = id;
		this->exp = exp;
		bool found = false;
		for (int i = 0; !found && i < (int)symbolTable.size(); i++) {
			if (symbolTable[i].symbol == id) {
				found = true;
			}
		}
		if (!found) {
			symbolTable.push_back(TableEntry(id, ""));
		}
	};
	void Print() {
		std::cout << this->id << ' ';
	};
};

class PrintStm : public Stm {
public:
	ExpList exps;
	PrintStm(ExpList exps) {
		this->exps = exps;
	};
};

class IdExp : public Exp {
public:
	std::string id;
	IdExp(std::string id) {
		this->id = id;
		Print();
		TableEntry te(id, "");
	};
	void Print() {
		std::cout << id << ' ';
	};
};

class NumExp : public Exp {
public:
	int num;
	NumExp(int num) {
		this->num = num;
		Print();
		char buffer[128] = { };
		_itoa_s(num, buffer, 127, 10);
		sStack.push(std::string(buffer));
	};
	void Print() {
		std::cout << num << ' ';
	};
};

enum class ArithmeticOp {
	Plus, Minus, Times, Div
};

class OpExp : public Exp {
public:
	Exp left, right;
	ArithmeticOp op;
	OpExp(Exp left, ArithmeticOp op, Exp right) {
		this->left = left;
		this->op = op;
		this->right = right;
		std::string ops = "";
		switch (op) {
		case ArithmeticOp::Plus:
			ops = "+";
			break;
		case ArithmeticOp::Minus:
			ops = "-";
			break;
		case ArithmeticOp::Times:
			ops = "*";
			break;
		case ArithmeticOp::Div:
			ops = "/";
			break;
		};
		std::cout << ops << std::endl;
		eStack.push(left);
		eStack.push(right);
		sStack.push(ops);
	};
};

class EseqExp : public Exp {
public:
	Stm stm; Exp exp;
	EseqExp(Stm stm, Exp exp) {
		this->stm = stm;
		this->exp = exp;
	};
};

class PairExpList : public ExpList {
public:
	Exp head;
	ExpList tail;
	PairExpList(Exp head, ExpList tail) {
		this->head = head;
		this->tail = tail;
	};
};

class LastExpList : public ExpList {
public:
	Exp head;
	LastExpList(Exp head) {
		this->head = head;
	};
};

#endif _SLPInterpreter_H

int main() {
	int a = 0, b = 0;
	Stm prog(CompoundStm(AssignStm("a",
		OpExp(NumExp(5), ArithmeticOp::Plus, NumExp(3))),
		CompoundStm(AssignStm("b",
			EseqExp(PrintStm(PairExpList(IdExp("a"),
				LastExpList(OpExp(IdExp("a"),
					ArithmeticOp::Minus, NumExp(1))))),
				OpExp(NumExp(10), ArithmeticOp::Times, IdExp("a")))),
			PrintStm(LastExpList(IdExp("b"))))));
	bool first = true;
	int result = 0;
	//sStack.push("0");
	while (!sStack.empty()) {
		std::string lts, ops, rts;
		if (first) {
			ops = sStack.top();
			sStack.pop();
			lts = sStack.top();
			sStack.pop();
			rts = sStack.top();
			sStack.pop();
			first = false;
		}
		else {
			lts = sStack.top();
			sStack.pop();
			ops = sStack.top();
			sStack.pop();
			rts = sStack.top();
			sStack.pop();
		}
		int lvi = std::stoi(lts);
		int rvi = std::stoi(rts);
		if (ops == "+") {
			result = lvi + rvi;
		}
		else if (ops == "-") {
			result = lvi - rvi;
		}
		else if (ops == "*") {
			result = lvi * rvi;
		}
		else if (ops == "/") {
			result = lvi / rvi;
		}
		char ascii[128] = { };
		_itoa_s(result, ascii, 10);
		if (sStack.size() != 0) {
			sStack.push(std::string(ascii));
		}
	}
	std::cout << "Result = " << result << std::endl;
	return 0;
}

Blog Entry © Sunday, March 29, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Slater Determinant Coefficients for Z = 2 to 4

Enter the atomic number Z (2 to 6 or 0 to quit): 2
2       1       1       +       a(1)b(2)
1       0       0       -       a(2)b(1)
# Even Permutations = 1
Enter the atomic number Z (2 to 6 or 0 to quit): 3
6       3       1       +       a(1)b(2)c(3)
5       2       0       -       a(1)b(3)c(2)
4       2       0       -       a(2)b(1)c(3)
3       1       1       +       a(2)b(3)c(1)
2       1       1       +       a(3)b(1)c(2)
1       0       0       -       a(3)b(2)c(1)
# Even Permutations = 3
Enter the atomic number Z (2 to 6 or 0 to quit): 4
24      12      0       +       a(1)b(2)c(3)d(4)
23      11      1       -       a(1)b(2)c(4)d(3)
22      11      1       -       a(1)b(3)c(2)d(4)
21      10      0       +       a(1)b(3)c(4)d(2)
20      10      0       +       a(1)b(4)c(2)d(3)
19      9       1       -       a(1)b(4)c(3)d(2)
18      9       1       -       a(2)b(1)c(3)d(4)
17      8       0       +       a(2)b(1)c(4)d(3)
16      8       0       +       a(2)b(3)c(1)d(4)
15      7       1       -       a(2)b(3)c(4)d(1)
14      7       1       -       a(2)b(4)c(1)d(3)
13      6       0       +       a(2)b(4)c(3)d(1)
12      6       0       +       a(3)b(1)c(2)d(4)
11      5       1       -       a(3)b(1)c(4)d(2)
10      5       1       -       a(3)b(2)c(1)d(4)
9       4       0       +       a(3)b(2)c(4)d(1)
8       4       0       +       a(3)b(4)c(1)d(2)
7       3       1       -       a(3)b(4)c(2)d(1)
6       3       1       -       a(4)b(1)c(2)d(3)
5       2       0       +       a(4)b(1)c(3)d(2)
4       2       0       +       a(4)b(2)c(1)d(3)
3       1       1       -       a(4)b(2)c(3)d(1)
2       1       1       -       a(4)b(3)c(1)d(2)
1       0       0       +       a(4)b(3)c(2)d(1)
# Even Permutations = 12
Enter the atomic number Z (2 to 6 or 0 to quit):

// AOPermutations.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Saturday, March 29, 2025
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Signs of the atomic orbitals in a Slater Determinant

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>

int main()
{
    char alpha[] = { 'a', 'b', 'c', 'd', 'e', 'f' }, line[128] = {};
    int factorial[7] = { 1, 1, 2, 6, 24, 120, 720 };

    while (true)
    {
        int col = 0, counter = 0, row = 0, sign = 1, t = 0, Z = 0, zfact = 0;
        int numberEven = 0;
        std::cout << "Enter the atomic number Z (2 to 6 or 0 to quit): ";
        std::cin.getline(line, 127);
        std::string str(line);
        Z = std::stoi(str);

        if (Z == 0)
        {
            break;
        }

        if (Z < 2 || Z > 6)
        {
            std::cout << "Illegal Z, please try again" << std::endl;
            continue;
        }

        zfact = factorial[Z];

        std::vector<char> orb(Z);
        std::vector<int> tmp(Z), vec(Z);

        for (int i = 0; i < Z; i++)
        {
            orb[i] = alpha[i];
            vec[i] = i + 1;
        }

        do
        {
            for (int i = 0; i < (int)vec.size(); i++)
            {
                tmp[i] = vec[i];
            }

            t = 0;

            do
            {
                t++;
            } while (std::next_permutation(tmp.begin(), tmp.end()));

            std::cout << t << '\t' << t / 2 << '\t';
            std::cout << (t / 2 & 1) << '\t';

            if (Z == 2 || Z == 3)
            {
                if ((t / 2 & 1) == 0)
                {
                    std::cout << "-\t";
                }

                else
                {
                    std::cout << "+\t";
                    numberEven++;
                }
            }

            else
            {
                if ((t / 2 & 1) == 1)
                {
                    std::cout << "-\t";
                }

                else
                {
                    std::cout << "+\t";
                    numberEven++;
                }
            }

            for (int i = 0; i < Z; i++)
            {
                std::cout << orb[i] << '(' << vec[i] << ')';
            }

            row++;
            std::cout << std::endl;

            if (zfact != 2 && row == zfact)
            {
                std::cout << std::endl;
                break;
            }

            row %= Z;
        } while (std::next_permutation(vec.begin(), vec.end()));

        std::cout << "# Even Permutations = ";
        std::cout << numberEven << std::endl;
    }

    return 0;
}

Blog Entry © Thursday, January 23, 2025, by James Pate Williams, Jr. Ackermann’s Super-Exponential Recursive Function in Vanilla C Programming Language

i = 2
j = 1
a(2, 1) =
4
# decimal digits = 1
enter another set (n to quit)? y
i = 2
j = 2
a(2, 2) =
16
# decimal digits = 2
enter another set (n to quit)? y
i = 2
j = 3
a(2, 3) =
65536
# decimal digits = 5
enter another set (n to quit)? y
i = 2
j = 4
a(2, 4) =
200352993040684646497907235156025575044782547556975141926501697371089\
405955631145308950613088093334810103823434290726318182294938211881266886\
950636476154702916504187191635158796634721944293092798208430910485599057\
015931895963952486337236720300291696959215610876494888925409080591145703\
767520850020667156370236612635974714480711177481588091413574272096719015\
183628256061809145885269982614142503012339110827360384376787644904320596\
037912449090570756031403507616256247603186379312648470374378295497561377\
098160461441330869211810248595915238019533103029216280016056867010565164\
...
506264233788565146467060429856478196846159366328895429978072254226479040\
061601975197500746054515006029180663827149701611098795133663377137843441\
619405312144529185518013657555866761501937302969193207612000925506508158\
327550849934076879725236998702356793102680413674571895664143185267905471\
716996299036301554564509004480278905570196832831363071899769915316667920\
895876857229060091547291963638167359667395997571032601557192023734858052\
112811745861006515259888384311451189488055212914577569914657753004138471\
712457796504817585639507289533753975582208777750607233944558789590571915\
6736
# decimal digits = 19729
enter another set (n to quit)?
/* 
** Computation of Akermann's super
** exponential function by James
** Pate Williams, Jr. (c) Tuesday,
** August 27, 2024 lip version
*/

#include <stdio.h>
#include "lip.h"

verylong Ackermann(verylong zi, verylong zj) {
	verylong a = 0;
	if (zscompare(zi, 1) == 0) {
		verylong ztwo = 0;
		zintoz(2, &ztwo);
		zexp(ztwo, zj, &a);
		return a;
	}
	else if (zscompare(zj, 1) == 0)
	{
		verylong ztwo = 0, ziminus1 = 0;
		zintoz(2, &ztwo);
		zsadd(zi, -1, &ziminus1);
		return Ackermann(ziminus1, ztwo);
	}
	else if (
		zscompare(zi, 2) >= 0 &&
		zscompare(zj, 2) >= 0) {
		verylong ziminus1 = 0;
		verylong zjminus1 = 0;
		verylong temp = 0;
		zsadd(zi, -1, &ziminus1);
		zsadd(zj, -1, &zjminus1);
		if (zscompare(ziminus1, 1) >= 0 &&
			zscompare(zjminus1, 1) >= 0) {
			return
				Ackermann(ziminus1, Ackermann(zi, zjminus1));
		}
	}
	return 0;
}

int DigitCount(verylong za) {
	int count = 0;
	while (zscompare(za, 0) > 0) {
		zsdiv(za, 10, &za);
		count++;
	}
	return count;
}

int main(void) {
	for (;;) {
		char buffer[256] = { '\0' };
		int i = 0, j = 0, number = 0;
		verylong za = 0, zi = 0, zj = 0;
		buffer[0] = '\0';
		printf_s("i = ");
		scanf_s("%d", &i);
		printf_s("j = ");
		scanf_s("%d", &j);
		zintoz(i, &zi);
		zintoz(j, &zj);
		printf_s("a(%d, %d) = \n", i, j);
		za = Ackermann(zi, zj);
		zwriteln(za);
		number = DigitCount(za);
		printf_s("# decimal digits = %d\n",
			number);
		printf_s("enter another set (n to quit)? ");
		scanf_s("%s", buffer, sizeof(buffer));
		zfree(&za);
		if (buffer[0] == 'n')
			break;
	}
	return 0;
}

Blog Entry © Thursday, January 23, 2025, by James Pate Williams, Jr. Merge Sort Verus Quick Sort

== Menu ==
1 Side-by-Side Tests
2 Timing Comparisons
3 Exit
Enter an option: 1
Enter a PRNG Seed >= 1: 1
0.001251        0.001251        0.001251
0.563585        0.014985        0.014985
0.193304        0.174108        0.174108
0.808740        0.193304        0.193304
0.585009        0.303995        0.303995
0.479873        0.350291        0.350291
0.350291        0.479873        0.479873
0.895962        0.513535        0.513535
0.822840        0.563585        0.563585
0.746605        0.585009        0.585009
0.174108        0.710501        0.710501
0.858943        0.746605        0.746605
0.710501        0.808740        0.808740
0.513535        0.822840        0.822840
0.303995        0.858943        0.858943
0.014985        0.895962        0.895962
== Menu ==
1 Side-by-Side Tests
2 Timing Comparisons
3 Exit
Enter an option:

== Menu ==
1 Side-by-Side Tests
2 Timing Comparisons
3 Exit
Enter an option: 2
Enter a PRNG Seed >= 1: 1
merge sort mean runtime (ms) = 523
quick sort mean runtime (ms) = 435
merge sort std dev (s) = 0.033457
quick sort std dev (s) = 0.027816
== Menu ==
1 Side-by-Side Tests
2 Timing Comparisons
3 Exit
Enter an option:
// MergeVersusQuick.c : This file contains the 'main' function.
// Program execution begins and ends there.
// mergesort is from "A Numerical Library in C for Scientists
// and Engineers" by H. T. Lau Translated from ALGOL NUMAL
// QuickSort is from "Introduction to Algorithms" by Thomas H.
// Cormen, Charles E. Leiserson, and Ronald L. Rivest

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

#define LENGTH1			17		// static side-by-side test
#define LENGTHM1		16		// upper index
#define LENGTH2		  4097		// time test array length
#define LENGTHM2      4096		// upper inndex
#define NUMBER_TESTS  4096		// number of timing tests

void mergesort(float a[], int p[], int low, int up)
{
	int* allocate_integer_vector(int, int);
	void free_integer_vector(int*, int);
	void merge(int, int, int, int[], float[], int[]);
	int i, lo, step, stap, umlp1, umsp1, rest, restv, * hp;

	hp = allocate_integer_vector(low, up);
	for (i = low; i <= up; i++) p[i] = i;
	restv = 0;
	umlp1 = up - low + 1;
	step = 1;
	do {
		stap = 2 * step;
		umsp1 = up - stap + 1;
		for (lo = low; lo <= umsp1; lo += stap)
			merge(lo, step, step, p, a, hp);
		rest = up - lo + 1;
		if (rest > restv && restv > 0)
			merge(lo, rest - restv, restv, p, a, hp);
		restv = rest;
		step *= 2;
	} while (step < umlp1);
	free_integer_vector(hp, low);
}

int* allocate_integer_vector(int l, int u)
{
	/*  Allocates an integer vector of range [l..u].  */

	void system_error(char*);
	int* p;

	p = (int*)malloc((unsigned)(u - l + 1) * sizeof(int));
	if (!p) system_error("Failure in allocate_integer_vector().");
	return p - l;
}

void free_integer_vector(int* v, int l)
{
	/*  Frees an integer vector of range [l..u].  */

	free((char*)(v + l));
}

void system_error(char error_message[])
{
	void exit(int);

	printf("%s", error_message);
	exit(1);
}

void merge(int lo, int ls, int rs, int p[], float a[], int hp[])
{
	/* this procedure is used internally by MERGESORT */

	int l, r, lout, rout, i, pl, pr;

	l = lo;
	r = lo + ls;
	lout = rout = 0;
	i = lo;
	do {
		pl = p[l];
		pr = p[r];
		if (a[pl] > a[pr]) {
			hp[i] = pr;
			r++;
			rout = (r == lo + ls + rs);
		}
		else {
			hp[i] = pl;
			l++;
			lout = (l == lo + ls);
		}
		i++;
	} while (!(lout || rout));
	if (rout) {
		for (i = lo + ls - 1; i >= l; i--) p[i + rs] = p[i];
		r = l + rs;
	}
	for (i = r - 1; i >= lo; i--) p[i] = hp[i];
}

int Partition(float a[], int lo, int hi)
{
	int pivotIndex = lo + (hi - lo) / 2;
	float x = a[pivotIndex];
	float 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 DoQuickSort(float a[], int n, int p, int r)
{
	if (p < r)
	{
		int q = Partition(a, p, r);

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

void QuickSort(float a[], int n)
{
	DoQuickSort(a, n, 1, n);
}

double runtime[2][NUMBER_TESTS];
float A[LENGTH1] = { 0 }, a[LENGTH1] = { 0 };
float AA[LENGTH2] = { 0 }, aa[LENGTH2] = { 0 };
float b[LENGTH1] = { 0 };
int n = NUMBER_TESTS;
int pp[LENGTH2];

int main()
{
	while (1)
	{
		float x;
		int i, j, option = 0, p[LENGTH1], seed = 1;

		printf("== Menu ==\n");
		printf("1 Side-by-Side Tests\n");
		printf("2 Timing Comparisons\n");
		printf("3 Exit\n");
		printf("Enter an option: ");
		scanf_s("%d", &option);

		if (option == 3)
			break;

		if (option == 1 || option == 2)
		{
			printf("Enter a PRNG Seed >= 1: ");
			scanf_s("%u", &seed);
			srand(seed);

			if (option == 1)
			{
				for (i = 1; i <= LENGTHM1; i++)
				{
					x = (float)rand() / RAND_MAX;
					A[i] = x;
					a[i] = x;
					b[i] = x;
				}

				mergesort(A, p, 1, LENGTHM1);
				QuickSort(a, LENGTHM1);

				for (i = 1; i <= LENGTHM1; i++)
					printf("%f\t%f\t%f\n", b[i], A[p[i]], a[i]);
			}

			else if (option == 2)
			{
				double mean[2] = { 0 }, median[2] = { 0 }, stdDev[2] = { 0 };
				double Sx[2] = { 0 };
				
				for (j = 0; j < n; j++)
				{
					for (i = 1; i <= LENGTHM2; i++)
					{
						x = (float)rand() / RAND_MAX;
						AA[i] = x;
						aa[i] = x;
					}
					
					clock_t start1 = clock();
					mergesort(AA, pp, 1, LENGTHM2);
					clock_t finis1 = clock();
					clock_t start2 = clock();
					QuickSort(aa, LENGTHM2);
					clock_t finis2 = clock();
					runtime[0][j] = ((double)finis1 - start1) / 
						CLOCKS_PER_SEC;
					runtime[1][j] = ((double)finis2 - start2) / 
						CLOCKS_PER_SEC;
					mean[0] += runtime[0][j];
					mean[1] += runtime[1][j];
				}

				for (j = 0; j < n; j++)
				{
					Sx[0] = pow((double)runtime[0][j] - mean[0], 2.0) / ((double)n - 1);
					Sx[1] = pow((double)runtime[1][j] - mean[1], 2.0) / ((double)n - 1);
				}

				stdDev[0] = (float)sqrt(Sx[0]);
				stdDev[1] = (float)sqrt(Sx[1]);

				printf("merge sort mean runtime (ms) = %3.0lf\n", 1.0e6 * mean[0] / n);
				printf("quick sort mean runtime (ms) = %3.0lf\n", 1.0e6 * mean[1] / n);
				printf("merge sort std dev (s) = %f\n", stdDev[0]);
				printf("quick sort std dev (s) = %f\n", stdDev[1]);
			}
		}
	}

	return 0;
}

Blog Entry (c) Wednesday, November 6, 2024, by James Pate Williams, Jr. Small Angular Momentum Quantum Numbers Gaunt Coefficients

// GauntCoefficients.cpp (c) Monday, November 4, 2024
// by James Pate Williams, Jr., BA, BS, MSWE, PhD
// Computes the Gaunt angular momentum coefficients
// Also the Wigner-3j symbols are calculated 
// https://en.wikipedia.org/wiki/3-j_symbol
// https://doc.sagemath.org/html/en/reference/functions/sage/functions/wigner.html#
// https://www.geeksforgeeks.org/factorial-large-number/
#include <iostream>
using namespace std;
typedef long double real;
real pi;
// iterative n-factorial function
real Factorial(int n)
{
    real factorial = 1;

    for (int i = 2; i <= n; i++)
        factorial *= i;
    if (n < 0)
        factorial = 0;
    return factorial;
}
real Delta(int lt, int rt)
{
    return lt == rt ? 1.0 : 0.0;
}
real Wigner3j(
    int j1, int j2, int j3,
    int m1, int m2, int m3)
{
    real delta = Delta(m1 + m2 + m3, 0) * 
        powl(-1.0, j1 - j2 - m3);
    real fact1 = Factorial(j1 + j2 - j3);
    real fact2 = Factorial(j1 - j2 + j3);
    real fact3 = Factorial(-j1 + j2 + j3);
    real denom = Factorial(j1 + j2 + j3 + 1);
    real numer = delta * sqrt(
        fact1 * fact2 * fact3 / denom);
    real fact4 = Factorial(j1 - m1);
    real fact5 = Factorial(j1 + m1);
    real fact6 = Factorial(j2 - m2);
    real fact7 = Factorial(j2 + m2);
    real fact8 = Factorial(j3 - m3);
    real fact9 = Factorial(j3 + m3);
    real sqrt1 = sqrtl(
        fact4 * fact5 * fact6 * fact7 * fact8 * fact9);
    real sumK = 0;
    int K = (int)fmaxl(0, fmaxl((real)j2 - j3 - m1,
        (real)j1 - j3 + m2));
    int N = (int)fminl((real)j1 + j2 - j3, 
        fminl((real)j1 - m1, (real)j2 + m2));
    for (int k = K; k <= N; k++)
    {
        real f0 = Factorial(k);
        real f1 = Factorial(j1 + j2 - j3 - k);
        real f2 = Factorial(j1 - m1 - k);
        real f3 = Factorial(j2 + m2 - k);
        real f4 = Factorial(j3 - j2 + m1 + k);
        real f5 = Factorial(j3 - j1 - m2 + k);
        sumK += powl(-1.0, k) / (f0 * f1 * f2 * f3 * f4 * f5);
    }
    return numer * sqrt1 * sumK;
}
real GauntCoefficient(
    int l1, int l2, int l3, int m1, int m2, int m3)
{
    real factor = sqrtl(
        (2.0 * l1 + 1.0) *
        (2.0 * l2 + 1.0) *
        (2.0 * l3 + 1.0) /
        (4.0 * pi));
    real wigner1 = Wigner3j(l1, l2, l3, 0, 0, 0);
    real wigner2 = Wigner3j(l1, l2, l3, m1, m2, m3);
    return factor * wigner1 * wigner2;
}
int main()
{
    pi = 4.0 * atanl(1.0);
    cout << "Gaunt(1, 0, 1, 1, 0, 0)  = ";
    cout << GauntCoefficient(1, 0, 1, 1, 0, 0);
    cout << endl;
    cout << "Gaunt(1, 0, 1, 1, 0, -1) = ";
    cout << GauntCoefficient(1, 0, 1, 1, 0, -1);
    cout << endl;
    real number = -1.0 / 2.0 / sqrtl(pi);
    cout << "-1.0 / 2.0 / sqrt(pi)    = ";
    cout << number << endl;
    return 0;
}

Blog Entry (c) Friday, October 18, 2024, by James Pate Williams, Jr. Ab Initio Quantum Chemical Calculation

On Wednesday, October 16, 2024, I bought an Amazon Kindle book named “Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory” by Attila Szabo and Neil S. Ostlund. It cost me $10.69 which is a real bargain. In Appendix B there is a listing for a FORTRAN program to perform an ab initio Hartree-Fock Self Consistent Field calculation for a two-electron heteronuclear molecule namely the helium-hydrogen cation. I successfully translated the program from FORTRAN to C++. I had to remember that FORTRAN stores matrices in column major order and C/C++ stores matrices in row major order. I took the transposes of two FORTRAN COMMON matrices to get the correct C++ storage. The authors of the book did an extensive treatment of the test calculation. The application is only 823 lines of monolithic C++ source code. I used FORTRAN like array indexing starting at 1 instead of the C initial beginning index of 0. I think I will try to get in touch with the authors to get permission to post the source code and results on my blog. 

P. S. I got permission from Dover Books to publish my source code and results. I think I will reconsider posting the C++ source code. The actual ground state energy of the cation is -2.97867. My calculation and the book’s computation are in percentage errors of about 4%. The book’s value is a little closer to the exact value than my result. The book calculation was done in FORTRAN double precision on a Digital Equipment Corporation PDP-10 minicomputer. My recreation of the book’s endeavor was executed on an Intel Itanium Core 7 and Windows 10 Professional machine using Win32 C++. The C++ compiler was from Microsoft Visual Studio 2019 Community Version.

Note I added a calculation for a homonuclear molecule, namely, the hydrogen diatomic molecule.

Ab Initio Quantum Chemical ResultsDownload

Blog Entry Friday, July 19, 2024, Easy Internet Math “Puzzle” (c) James Pate Williams, Jr.

Internet Math PuzzleDownload
#include <math.h>
#include <iostream>
using namespace std;

long double f(long double x)
{
	return powl(8.0, x) - powl(2.0, x) -
		2.0 * (powl(6.0, x) - powl(3.0, x));
}

long double g(long double x)
{
	return powl(8.0, x) * logl(8.0) - powl(2.0, x) * logl(2.0) -
		2.0 * (powl(6.0, x) * logl(6.0) - powl(3.0, x) * logl(3.0));
}

long double Newton(long double x, int maxIts, int& iterations)
{
	long double x0 = x;
	long double x1 = 0.0;
	
	iterations = 0;

	while (true) {
		long double dx = 0.0;
		long double fx = f(x0);
		long double gx = g(x0);
		x1 = x0 - fx / gx;
		dx = fabsl(x1 - x0);
		iterations++;
		if (dx < 1.0e-15)
			break;
		if (fabsl(fx) < 1.0e-15)
			break;
		if (iterations == maxIts)
			break;
		x0 = x1;
	}

	return x1;
}

int main() {
	int iterations = 0, maxIts;
	long double x0 = 0.0, x1 = 0.0;

	while (true) {
		cout << "x0 = ";
		cin >> x0;
		if (x0 == 0)
			break;
		cout << "maximum iterations = ";
		cin >> maxIts;
		x1 = Newton(x0, maxIts, iterations);
		cout << "x1 = " << x1 << endl;
		cout << "iterations = ";
		cout << iterations << endl;
	}

	return 0;
}

Blog Entry Monday, June 24, 2024 (c) James Pate Williams, Jr. Computing Binomial Coefficients and Pascal’s Triangle in the C Language

Enter n (<= 18) below:
5

Enter k (<= 18) below:
0

1       1

Enter n (<= 18) below:
5

Enter k (<= 18) below:
1

5       5

Enter n (<= 18) below:
5

Enter k (<= 18) below:
2

10      10

Enter n (<= 18) below:
0
Enter n (<= 18) below:
0

Pascal's Triangle:

    1
    1     1
    1     2     1
    1     3     3     1
    1     4     6     4     1
    1     5    10    10     5     1
    1     6    15    20    15     6     1
    1     7    21    35    35    21     7     1
    1     8    28    56    70    56    28     8     1
    1     9    36    84   126   126    84    36     9     1
    1    10    45   120   210   252   210   120    45    10     1
    1    11    55   165   330   462   462   330   165    55    11     1
    1    12    66   220   495   792   924   792   495   220    66    12     1
    1    13    78   286   715  1287  1716  1716  1287   715   286    78    13     1
    1    14    91   364  1001  2002  3003  3432  3003  2002  1001   364    91    14     1
    1    15   105   455  1365  3003  5005  6435  6435  5005  3003  1365   455   105    15     1
    1    16   120   560  1820  4368  8008 11440 12870 11440  8008  4368  1820   560   120    16     1
    1    17   136   680  2380  6188 12376 19448 24310 24310 19448 12376  6188  2380   680   136    17     1
    1    18   153   816  3060  8568 18564 31824 43758 48620 43758 31824 18564  8568  3060   816   153    18     1

C:\Users\james\source\repos\BinomialCoefficeint\Debug\BinomialCoefficeint.exe (process 40028) exited with code 0.
Press any key to close this window . . .
// BinomialCoefficient.c (c) Monday, June 24, 2024
// by James Pate Williams, Jr. BA, BS, MSwE, PhD

#include <stdio.h>
#include <stdlib.h>
typedef long long ll;

ll** Binomial(ll n)
{
    ll** C = (ll**)calloc(n + 1, sizeof(ll*));

    if (C == NULL)
        exit(-1);

    for (int i = 0; i < n + 1; i++)
    {
        C[i] = (ll*)calloc(n + 1, sizeof(ll));

        if (C[i] == NULL)
            exit(-1);
    }

    if (n >= 0)
    {
        C[0][0] = 1;
    }

    if (n >= 1)
    {
        C[1][0] = 1;
        C[1][1] = 1;
    }

    if (n >= 2)
    {
        for (int i = 2; i <= n; i++)
        {
            for (int j = 2; j <= n; j++)
            {
                C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
            }
        }
    }

    return C;
}

ll Factorial(ll n)
{
    ll fact = 1;

    if (n > 1)
    {
        for (int i = 2; i <= n; i++)
            fact = i * fact;
    }

    return fact;
}

ll BC(ll n, ll k)
{
    return Factorial(n) / (Factorial(n - k) * Factorial(k));
}

int main()
{
    int i = 0, j = 0;
    ll** C = Binomial(20);

    while (1)
    {
        char buffer[256] = { '\0' };
        
        printf_s("Enter n (<= 18) below:\n");
        scanf_s("%s", buffer, sizeof(buffer));
        printf_s("\n");

        ll n = atoll(buffer);

        if (n == 0)
            break;

        printf_s("Enter k (<= 18) below:\n");
        scanf_s("%s", buffer, sizeof(buffer));
        printf_s("\n");

        ll k = atoll(buffer);
                
        printf_s("%lld\t%lld\n\n", C[n + 2][k + 2], BC(n, k));
    }

    printf_s("Pascal's Triangle:\n\n");

    for (i = 2; i <= 20; i++)
    {
        for (j = 2; j <= 20; j++)
            if (C[i][j] != 0)
                printf_s("%5lld ", C[i][j]);

        printf_s("\n");
    }

    for (i = 0; i <= 20; i++)
        free(C[i]);

    free(C);
}

Blog Entry Sunday, June 23, 2024 (c) James Pate Williams, Jr.

The object of this C Win32 application is to find a multiple of 9 with its digits summing to a multiple of 9 also. The first column below is a multiple of 9 whose digits sum to 9 also. The second column is the sum of digits found in the column one number. The last column is the first column divided by 9.

Enter PRNG seed:
1
Enter number of bits (4 to 16):
4
    9       9       1
Enter number of bits (4 to 16):
5
   27       9       3
Enter number of bits (4 to 16):
6
   45       9       5
Enter number of bits (4 to 16):
7
  117       9      13
Enter number of bits (4 to 16):
8
  252       9      28
Enter number of bits (4 to 16):
0

C:\Users\james\source\repos\CProductOf9Console\Debug\CProductOf9Console.exe (process 23280) exited with code 0.
Press any key to close this window . . .
Enter PRNG seed:
1
Enter number of bits (4 to 16):
9
  369      18      41
Enter number of bits (4 to 16):
10
  846      18      94
Enter number of bits (4 to 16):
11
 1080       9     120
Enter number of bits (4 to 16):
12
 3015       9     335
Enter number of bits (4 to 16):
13
 5040       9     560
Enter number of bits (4 to 16):
14
10350       9    1150
Enter number of bits (4 to 16):
15
30870      18    3430
Enter number of bits (4 to 16):
16
57798      36    6422
Enter number of bits (4 to 16):
0
// CProductOf9Console.c (c) Sunday, June 23, 2024
// by James Pate Williams, Jr., BA, BS, MSwE, PhD

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

char nextStr[256], numbStr[256];

void ConvertToString(int number, int radix)
{
	int i = 0;

	while (number > 0)
	{
		nextStr[i++] = (char)(number % radix + '0');
		number /= radix;
	}

	nextStr[i++] = '\0';
	_strrev(nextStr);
}

int Sum(int next)
{
	long sum = 0;

	ConvertToString(next, 10);

	for (int i = 0; i < (int)strlen(nextStr); i++)
		sum += (long)nextStr[i] - '0';

	if (sum % 9 == 0 && sum != 0)
		return sum;

	return -1;
}

long GetNext(int numBits, int* next)
{
	long hi = 0, lo = 0, nine = 0;

	nextStr[0] = '\0';
	numbStr[0] = '\0';

	if (numBits == 4)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 16);

			if (*next != 0 && *next >= 8 && *next < 16)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 5)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 32);

			if (*next >= 16 && *next < 32)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 6)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 64);

			if (*next >= 32 && *next < 64)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 7)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 128);

			if (*next >= 64 && *next < 128)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 8)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 256);

			if (*next >= 128 && *next < 256)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 9)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 512);

			if (*next >= 256 && *next < 512)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 10)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 1024);

			if (*next >= 512 && *next < 1024)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 11)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 2048);

			if (*next >= 1024 && *next < 2048)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 12)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 4096);

			if (*next >= 2048 && *next < 4096)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 13)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 8192);

			if (*next >= 4096 && *next < 8192)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 14)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 16384);

			if (*next >= 8192 && *next < 16384)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 15)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 32768);

			if (*next >= 16384 && *next < 32768)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 16)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 65536);

			if (*next >= 32768 && *next < 65536)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	return -1;
}

int main()
{
	char buffer[256] = { '\0' };
	long seed = 0;

	printf_s("Enter PRNG seed:\n");
	scanf_s("%s", buffer, sizeof(buffer));
	seed = atol(buffer);
	srand((unsigned int)seed);

	while (1)
	{
		int next = 0, nine = 0, numberBits = 0;

		printf_s("Enter number of bits (4 to 16):\n");
		scanf_s("%s", buffer, sizeof(buffer));
		numberBits = atol(buffer);

		if (numberBits == 0)
			break;

		if (numberBits < 4 || numberBits > 16)
		{
			printf_s("illegal number of bits must >= 4 and <= 16\n");
			continue;
		}

		nine = GetNext(numberBits, &next);

		if (nine == -1)
		{
			printf_s("illegal result, try again\n");
			continue;
		}

		printf_s("%5ld\t%5ld\t%5ld\n", next, nine, next / 9);
	}

	return 0;
}