// 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;
}
Tag: programming
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.
Blog Entry Friday, July 19, 2024, Easy Internet Math “Puzzle” (c) James Pate Williams, Jr.
#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;
}
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;
}