Tag: programming

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

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;
}