Category: Uncategorized

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

Exercises Szabo OstlundDownload
#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 Wednesday, July 10, 2024, © James Pate Williams, Jr. My Dual Interests in Cryptography and Number Theory

I became fascinated with secret key cryptography as a child. Later, as an adult, in around 1979, I started creating crude symmetric cryptographic algorithms. I became further enthralled with cryptography and number theory in 1996 upon reading Applied Cryptography, Second Edition: Protocols, Algorithms, and Source Code in C by Bruce Schneier and later the Handbook of Applied Cryptography by Alfred J. Menezes, Paul C. van Oorschot, and Scott A. Vanstone. After implementing many of the algorithms in both tomes, I communicated my results to two of the authors namely Bruce Schneier and Professor Alfred J. Menezes. In 1997 I developed a website devoted to constraint satisfaction problems and their solutions, cryptography, and number theory. I posted legal C and C++ source code. Professor Menezes advertised my website along with his treatise. See the following blurb:

In the spirit of my twin scientific infatuations, I offer yet another C integer factoring implementation utilizing the Free Large Integer Package (known more widely as lip) which was created by Arjen K. Lenstra (now a Professor Emeritus). This implementation includes Henri Cohen’s Trial Division algorithm, the Brent-Cohen-Pollard rho method, the Cohen-Pollard p – 1 stage 1 method, and the Lenstra lip Elliptic Curve Method. If I can get the proper authorization, I will later post the source code.

total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
1
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.014000 seconds:
2
5 ^ 2
41
397
2113
enter number below:
0

total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
2
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 1.531000 seconds:
2
5 ^ 2
41
397
2113
415878438361
3630105520141
enter number below:
0

total time required for initialization: 0.055000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
3
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.066000 seconds:
2
5 ^ 2
41
838861
415878438361
3630105520141
enter number below:
0

total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
4
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.013000 seconds:
2
5
205
838861
415878438361
3630105520141
enter number below:
0

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

Blog Entry (c) Friday, June 21, 2024, by James Pate Williams, Jr. Comparison of Two Prime Number Sieves

First the C++ results:

Limit = 1000000
Number of primes <= 1000000 78498
Milliseconds taken by Sieve of Atkin: 12
Number of primes <= 1000000 78498
Milliseconds taken by Sieve of Eratosthenes: 14
Limit = 10000000
Number of primes <= 10000000 664579
Milliseconds taken by Sieve of Atkin: 159
Number of primes <= 10000000 664579
Milliseconds taken by Sieve of Eratosthenes: 204
Limit = 100000000
Number of primes <= 100000000 5761455
Milliseconds taken by Sieve of Atkin: 1949
Number of primes <= 100000000 5761455
Milliseconds taken by Sieve of Eratosthenes: 2343
Limit = 0

Next, we have the Java results:

C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 1000000 0
number of primes less than equal 1000000 = 78498
total computation time in seconds = 0.008

C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 10000000 0
number of primes less than equal 10000000 = 664579
total computation time in seconds = 0.098

C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 1000000 0
number of primes less than equal 1000000 = 78498
total computation time in seconds = 0.011

C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 10000000 0
number of primes less than equal 10000000 = 664579
total computation time in seconds = 0.151

C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 100000000 0
number of primes less than equal 100000000 = 5761455
total computation time in seconds = 1.511

C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 100000000 0
number of primes less than equal 100000000 = 5761455
total computation time in seconds = 1.995

Notice that the Java application outperforms the C++ application.

// PrimeSieveComparison.cpp (c) Friday, June 21, 2024
// by James Pate Williams, Jr.
//
//  SieveOfAtkin.java
//  SieveOfAtkin
//
//  Created by James Pate Williams, Jr. on 9/29/07.
//  Copyright (c) 2007 James Pate Williams, Jr. All rights reserved.
//
//  SieveOfEratosthenes.java
//  SieveOfEratosthenes
//
//  Created by James Pate Williams, Jr. on 9/29/07.
//  Copyright (c) 2007 James Pate Williams, Jr. All rights reserved.
//

#include <math.h>
#include <iostream>
#include <chrono>
using namespace std::chrono;
using namespace std;

const int Maximum = 100000000;
bool sieve[Maximum + 1];

void SieveOfAtkin(int limit)
{
	auto start = high_resolution_clock::now();
	int e, k, n, p, x, xx3, xx4, y, yy;
	int primeCount = 2, sqrtLimit = (int)sqrt(limit);

	for (n = 5; n <= limit; n++)
		sieve[n] = false;

	for (x = 1; x <= sqrtLimit; x++) {
		xx3 = 3 * x * x;
		xx4 = 4 * x * x;
		for (y = 1; y <= sqrtLimit; y++) {
			yy = y * y;
			n = xx4 + yy;
			if (n <= limit && (n % 12 == 1 || n % 12 == 5))
				sieve[n] = !sieve[n];
			n = xx3 + yy;
			if (n <= limit && n % 12 == 7)
				sieve[n] = !sieve[n];
			n = xx3 - yy;
			if (x > y && n <= limit && n % 12 == 11)
				sieve[n] = !sieve[n];
		}
	}

	for (n = 5; n <= sqrtLimit; n++) {
		if (sieve[n]) {
			e = 1;
			p = n * n;
			while (true) {
				k = e * p;
				if (k > limit)
					break;
				sieve[k] = false;
				e++;
			}
		}
	}
	
	for (n = 5; n <= limit; n++)
		if (sieve[n])
			primeCount++;

	auto stop = high_resolution_clock::now();
	auto duration = duration_cast<milliseconds>(stop - start);

	std::cout << "Number of primes <= " << limit << ' ';
	std::cout << primeCount << endl;
	std::cout << "Milliseconds taken by Sieve of Atkin: "
		<< duration.count() << endl;
}

void SieveOfEratosthenes(int limit)
{
	auto start = high_resolution_clock::now();
	int i = 0, k = 0, n = 0, nn = 0;
	int primeCount = 0, sqrtLimit = (int)sqrt(limit);

	// initialize the prime number sieve

	for (n = 2; n <= limit; n++)
		sieve[n] = true;

	// eliminate the multiples of n

	for (n = 2; n <= sqrtLimit; n++)
		for (i = 2; i <= n - 1; i++)
			sieve[i * n] = false;

	// eliminate squares

	for (n = 2; n <= sqrtLimit; n++) {
		if (sieve[n]) {
			k = 0;
			nn = n * n;
			i = nn + k * n;
			while (i <= limit) {
				sieve[i] = false;
				i = nn + k * n;
				k++;
			}
		}
	}

	primeCount = 0;

	for (n = 2; n <= limit; n++)
		if (sieve[n])
			primeCount++;

	auto stop = high_resolution_clock::now();
	auto duration = duration_cast<milliseconds>(stop - start);

	std::cout << "Number of primes <= " << limit << ' ';
	std::cout << primeCount << endl;
	std::cout << "Milliseconds taken by Sieve of Eratosthenes: "
		<< duration.count() << endl;
}

int main()
{
	while (true)
	{
		int limit = 0;
		std::cout << "Limit = ";
		cin >> limit;

		if (limit == 0)
			break;

		SieveOfAtkin(limit);
		SieveOfEratosthenes(limit);
	}

	return 0;
}

Factorizations of Some Fibonacci Sequence Numbers, Lucas Sequence Numbers and Some Other Numbers Using Arjen K. Lenstra’s Free Large Integer Package and the Elliptic Curve Method (c) January 28, 2024, by James Pate Williams, Jr.

All of the following computations were performed on a late 2015 Dell XPS 8900 personal computer with a 64-bit Intel Core I7 processor @ 4.0GHz with 16GB of DDR2 RAM.

Factorization of Six Fibonacci Sequence Numbers:

 

Fibonacci 500

# digits 105

5 ^ 2 p # digits 1

15 c # digits 2

101 p # digits 3

401 p # digits 3

1661 c # digits 4

3001 p # digits 4

10291 c # digits 5

570601 p # digits 6

112128001 p # digits 9

1353439001 p # digits 10

28143378001 p # digits 11

5465167948001 p # digits 13

84817574770589638001 p # digits 20

158414167964045700001 p # digits 21

Runtime (s) = 1.206000



Fibonacci 505

# digits 106

5 p # digits 1

743519377 p # digits 9

44614641121 p # digits 11

770857978613 p # digits 12

960700389041 p # digits 12

12588421794766514566269164716286291055826556238643852856601641 p # digits 62

Runtime (s) = 1.959000



Fibonacci 510

# digits 107

2 ^ 3 p # digits 1

11 p # digits 2

61 p # digits 2

1021 p # digits 4

1597 p # digits 4

3469 p # digits 4

3571 p # digits 4

9521 p # digits 4

53551 p # digits 5

95881 p # digits 5

142445 c # digits 6

1158551 p # digits 7

3415914041 p # digits 10

20778644396941 p # digits 14

20862774425341 p # digits 14

81358225616651 c # digits 14

162716451241291 p # digits 15

Runtime (s) = 2.682000



Fibonacci 515

# digits 108

5 p # digits 1

519121 p # digits 6

5644193 p # digits 7

512119709 p # digits 9

84388938382141 p # digits 14

300367026458796424297447559250634818495937628065437243817852436228914621 p # digits 72

Runtime (s) = 7.861000



Fibonacci 520

# digits 109

131 p # digits 3

451 c # digits 3

521 p # digits 3

2081 p # digits 4

2161 p # digits 4

3121 p # digits 4

24571 p # digits 5

90481 p # digits 5

2519895 c # digits 7

21183761 p # digits 8

57089761 p # digits 8

102193207 p # digits 9

1932300241 p # digits 10

14736206161 p # digits 11

5836312049326721 p # digits 16

42426476041450801 p # digits 17

Runtime (s) = 5.155000



Fibonacci 525

# digits 110

2 p # digits 1

5 p # digits 1

421 p # digits 3

701 p # digits 3

3001 p # digits 4

3965 c # digits 4

4201 p # digits 4

141961 p # digits 6

2553601 p # digits 7

230686501 p # digits 9

8288823481 p # digits 10

82061511001 p # digits 11

19072991752501 c # digits 14

8481116649425701 p # digits 16

17231203730201189308301 p # digits 23

Runtime (s) = 2.026000

Factorization of Six Lucas Sequence Numbers



Lucas 340

113709744839525149336680459091826532688903186653162057995534262332121127

# digits 72

7 p # digits 1

2161 p # digits 4

5441 p # digits 4

897601 p # digits 6

23230657239121 p # digits 14

17276792316211992881 p # digits 20

3834936832404134644974961 p # digits 25

Runtime (s) = 109.103000



Lucas 345

# digits 73

2 ^ 2 p # digits 1

31 p # digits 2

461 p # digits 3

1151 p # digits 4

1529 c # digits 4

324301 p # digits 6

686551 p # digits 6

1485571 p # digits 7

4641631 p # digits 7

19965899801 c # digits 11

117169733521 p # digits 12

3490125311294161 p # digits 16

Runtime (s) = 0.032000



Lucas 350

13985374084677485786380981408251904922622980674054858121032362563653278123

# digits 74

3 p # digits 1

401 p # digits 3

2801 p # digits 4

11521 c # digits 5

28001 p # digits 5

570601 p # digits 6

12317523121 p # digits 11

248773766357061401 p # digits 18

7358192362316341243805801 p # digits 25

Runtime (s) = 21.047000



Lucas 355

69362907070206748494476200566565775354902428015845969798000696945226974645

# digits 74

5 p # digits 1

4261 p # digits 4

6673 p # digits 4

75309701 p # digits 8

309273161 p # digits 9

46165371073 p # digits 11

9207609261398081 p # digits 16

49279722643391864192801 p # digits 23

Runtime (s) = 40.726000



Lucas 360

769246427201094785080787978422393713094534885688979999504447628313150135520

# digits 75

2 ^ 5 p # digits 1

3 ^ 2 p # digits 1

23 p # digits 2

41 p # digits 2

105 c # digits 3

107 p # digits 3

241 p # digits 3

2161 p # digits 4

2521 p # digits 4

3439 c # digits 4

8641 p # digits 4

20641 p # digits 5

103681 p # digits 6

109441 p # digits 6

191306797 c # digits 9

10783342081 p # digits 11

13373763765986881 p # digits 17

Runtime (s) = 0.032000



Lucas 365

19076060504701386559675231910437330047906343529583769121365013189782992678011

# digits 77

11 p # digits 2

151549 p # digits 6

514651 p # digits 6

7015301 p # digits 7

8942501 p # digits 7

9157663121 p # digits 10

11899937029 p # digits 11

3252336525249736694804553589211 p # digits 31

The following two numbers were first factorized by J. M. Pollard on an 8-bit Phillips P2012 personal computer with 64 KB RAM and two 640 KB disc drives. The times required by Pollard were 41 and 47 hours.

2^144-3
22300745198530623141535718272648361505980413

# digits 44

492729991333 p # digits 12

45259565260477899162010980272761 p # digits 32

Runtime (s) = 0.086000

2^153+3
11417981541647679048466287755595961091061972995

# digits 47

5 p # digits 1

11 p # digits 2

600696432006490087537 p # digits 21

345598297796034189382757 p # digits 24

Runtime (s) = 0.676000

Partial factorization of the Twelfth Fermat Number 2^4096+1
# digits 1234

114689 p # digits 6

26017793 p # digits 8

63766529 p # digits 8

190274191361 p # digits 12

Runtime (s) = 1532.878000

My Exploration of the Computer-Generated Random Song Space (c) January 19, 2024, by James Pate Williams, Jr.

I have been interested in creating “music” using computers since I read an article by Martin Gardner in “Scientific American” way back in the 1970s. Gardner mentioned three types of noise found in nature: Brownian, Fractal, and White noises. My first feeble attempts at computer generated music happened in the early 1980s when I used LaGrange College’s Personal Computers (generic not IBM PCs). Later after my father bought me a Commodore Amiga 2000 on Saturday, April 30, 1988, I explored the sound producing space of that brand of Personal Computers and the Amiga BASIC computer language which was created by Microsoft. I designed and implemented a keyboard emulator program and a stochastic “music” generating program that used Brownian, Fractal, and White noise.

Skip forward to the early 2000s. While I was a graduate student in software engineering and computer science in the period 2000 to 2005, I translated Sun Java’s MIDI (Musical Instrument Digital Interface) sequencer to the Apple and PC C languages. I also created a MIDI language assembler and compiler. MIDI is well-known as a gestural music language. Sometime in the late 2000s, I designed and implemented an application in C# to create a random song from either musical chords or scales. Nowadays I still reuse a lot of my decade’s old software.

User interface of my “Random Song from Scale” C# application:

After creating a random song, I upload (import) the MIDI file to my SONAR Platinum Digital Audio Workstation then use a synthesizer to create a sound (audio) track to export as a wave, MP3, etc. file.

On January 17, 2024, I purchased two bargain basement Universal Audio VST plug-ins: a Moog Mini-Moog synthesizer emulator and a Waterfall Hammond B3 emulator complete with a software Type 147 amplifier with Doppler effect Lesley rotating speaker. Each plug-in cost me a cool $49.00 USD. Now the prices have jumped to $199.00. The real 1970s analog Mini-Moog is available for approximately $3,500 to $24,000 USD. I know an individual with a Voyager Moog synthesizer that was hand signed by the inventor Robert Moog.

https://www.musiciansfriend.com/pro-audio/universal-audio-moog-minimoog–uad-instrument-mac-windows/m00313000000000?rNtt=Universal%20Audio%20Minimoog%20plug-in&index=5

My 1988 Commodore Amiga 2000 Is Still Functioning by James Pate Williams, Jr.

On Friday September 8, 2023, I setup my thirty-five-year-old personal computer which is a Commodore Amiga 2000. My father bought the computer for me on Saturday, April 30, 1988. It still works although the Commodore 1084 RGB color display has a non-functional power button. My remedy for the problem was to Scotch tape the power button in the “On” position. I have an Amiga BASIC by Microsoft manual and software, a Modula-2 compiler, Motorola MC68000 macro-assembly language software, and a Pecan UCSD Pascal compiler. The Amiga was the first multimedia personal computer. In May 1988 I created two Amiga BASIC programs: a keyboard emulator and a primitive computer-generated music program that used three types of noise namely Brownian, fractal, and white noise. I used the computer extensively in the years 1988 to 1994. In December 1994 my mother and older sister purchased a mom-and-pop store Microsoft-Intel personal computer.

Numerical Two-Dimensional Quadrature (Integration) by James Pate Williams, Jr.

Back in 2022 I performed some numerical experiments using the two-dimensional quadrature methods: Simpson’s Rule and Monte Carlo. I implemented both sequential and multitasking techniques. The integration algorithms were applied to electron-electron repulsion integrals and a two-dimensional exponential function. The two-dimensional exponential function has an exact value of pi divided by 4. The multitasking speed-up is about two-fold in the Simpson’s Rule case and almost three-fold in the Monte Carlo case.