Author: jamespatewilliamsjr
My whole legal name is James Pate Williams, Jr. I was born in LaGrange, Georgia approximately 70 years ago. I barely graduated from LaGrange High School with low marks in June 1971. Later in June 1979, I graduated from LaGrange College with a Bachelor of Arts in Chemistry with a little over a 3 out 4 Grade Point Average (GPA). In the Spring Quarter of 1978, I taught myself how to program a Texas Instruments desktop programmable calculator and in the Summer Quarter of 1978 I taught myself Dayton BASIC (Beginner's All-purpose Symbolic Instruction Code) on LaGrange College's Data General Eclipse minicomputer. I took courses in BASIC in the Fall Quarter of 1978 and FORTRAN IV (Formula Translator IV) in the Winter Quarter of 1979. Professor Kenneth Cooper, a genius poly-scientist taught me a course in the Intel 8085 microprocessor architecture and assembly and machine language. We would hand assemble our programs and insert the resulting machine code into our crude wooden box computer which was designed and built by Professor Cooper. From 1990 to 1994 I earned a Bachelor of Science in Computer Science from LaGrange College. I had a 4 out of 4 GPA in the period 1990 to 1994. I took courses in C, COBOL, and Pascal during my BS work. After graduating from LaGrange College a second time in May 1994, I taught myself C++. In December 1995, I started using the Internet and taught myself client-server programming. I created a website in 1997 which had C and C# implementations of algorithms from the "Handbook of Applied Cryptography" by Alfred J. Menezes, et. al., and some other cryptography and number theory textbooks and treatises.
Capo 2nd Joseph (Joe) Gay Lead Guitar and James Pate Williams Jr Rythm Guitar a la May 26, 2009 or Earlier
I used a Fender 12-string acoustic electric guitar.
Vector Analysis by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition© 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Pages 48 and 50
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider
#include <iostream>
#include <vector>
static double InnerProduct(
std::vector<double> A,
std::vector<double> B,
int n)
{
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += A[i] * B[i];
return sum;
}
static void VectorProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double>&C)
{
C.resize(3);
C[0] = A[1] * B[2] - A[2] * B[1];
C[1] = A[2] * B[0] - A[0] * B[2];
C[2] = A[0] * B[1] - A[1] * B[0];
}
static double TripleProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double> C)
{
double sum0 = 0.0, sum1 = 0.0;
sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
return sum0 - sum1;
}
static void Exercises_Section_1_13_1_Triple_Products()
{
// triple products (a)
std::vector<double> A1 = { 2, 0, 0 };
std::vector<double> B1 = { 0, 3, 0 };
std::vector<double> C1 = { 0, 0, 5 };
double tp_1 = TripleProduct(A1, B1, C1);
std::cout << "Exercise 1 (a)" << '\t';
std::cout << tp_1 << std::endl;
// triple products (b)
std::vector<double> A2 = { 1, 1, 1 };
std::vector<double> B2 = { 3, 1, 0 };
std::vector<double> C2 = { 0, -1, 5 };
std::cout << "Exercise 1 (b)" << '\t';
double tp_2 = TripleProduct(A2, B2, C2);
std::cout << tp_2 << std::endl;
// triple products (c)
std::vector<double> A3 = { 2, -1, 1 };
std::vector<double> B3 = { 1, 1, 1 };
std::vector<double> C3 = { 2, 0, 3 };
double tp_3 = TripleProduct(A3, B3, C3);
std::cout << "Exercise 1 (c)" << '\t';
std::cout << tp_3 << std::endl;
// triple products (d)
std::vector<double> A4 = { 0, 0, 1 };
std::vector<double> B4 = { 1, 0, 0 };
std::vector<double> C4 = { 0, 1, 0 };
double tp_4 = TripleProduct(A4, B4, C4);
std::cout << "Exercise 1 (d)" << '\t';
std::cout << tp_4 << std::endl;
// volume of a parallelpipped
std::vector<double> A5 = { 3, 4, 1 };
std::vector<double> B5 = { 2, 3, 4 };
std::vector<double> C5 = { 0, 0, 5 };
double tp_5 = TripleProduct(A5, B5, C5);
std::cout << "Exercise 2" << '\t';
std::cout << tp_5 << std::endl;
// volume of a parallelpipped
std::vector<double> A6 = { 3, 2, 1 };
std::vector<double> B6 = { 4, 2, 1 };
std::vector<double> C6 = { 0, 1, 4 };
std::vector<double> D6 = { 0, 0, 7 };
std::vector<double> AB = { -1, 0, 0 };
std::vector<double> AC = { 3, 1, -3 };
std::vector<double> AD = { 3, 2, -6 };
std::cout << "Exercise 3" << '\t';
double tp_6 = TripleProduct(AB, AC, AD);
std::cout << tp_6 << std::endl;
// volume of a tetrahedron
std::vector<double> AB1 = { 1, 1, 0 };
std::vector<double> AC1 = { 1, -1, 0 };
std::vector<double> AD1 = { 0, 0, 2 };
std::cout << "Exercise 4" << '\t';
double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
std::cout << fabs(tp_7) << std::endl;
std::vector<double> P1 = { 0, 0, 0 };
std::vector<double> P2 = { 1, 1, 0 };
std::vector<double> P3 = { 3, 4, 0 };
std::vector<double> P4 = { 4, 5, 0 };
std::vector<double> P5 = { 0, 0, 1 };
std::vector<double> Q1 = { 1, 1, 0 };
std::vector<double> Q2 = { 2, 3, 0 };
std::vector<double> Q3 = { 1, 1, 1 };
std::cout << "Exercise 5" << '\t';
double tp_8 = TripleProduct(Q1, Q2, Q3);
std::cout << fabs(tp_8) << std::endl;
std::vector<double> A10 = { 1, 1, 1 };
std::vector<double> B10 = { 2, 4, -1 };
std::vector<double> C10 = { 1, 1, 3 };
double tp_10 = TripleProduct(A10, B10, C10);
std::vector<double> D10 = { 0 };
VectorProduct(A10, B10, D10);
double magnitude = sqrt(InnerProduct(D10, D10, 3));
std::cout << "Exercise 10" << '\t';
std::cout << tp_10 / magnitude << '\t';
std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
std::vector<double> A11 = { 1, 1, 1 };
std::vector<double> B11 = { 2, 4, -1 };
std::vector<double> C11 = { 1, 1, 3 };
std::vector<double> D11 = { 3, 2, 1 };
std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
VectorProduct(A11, B11, AB11);
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, BCCA11);
double Q11 = InnerProduct(AB11, BCCA11, 3);
double A_x = A11[0], A_y = A11[1], A_z = A11[2];
double B_x = B11[0], B_y = B11[1], B_z = B11[2];
double C_x = C11[0], C_y = C11[1], C_z = C11[2];
double term1 = +(A_y * B_z - A_z * B_y) * (B_z * C_x - B_x * C_z) * (C_x * A_y - C_y * A_x);
double term2 = -(A_y * B_z - A_z * B_y) * (B_x * C_y - B_y * C_x) * (C_z * A_x - C_x * A_z);
double term3 = +(A_z * B_x - A_x * B_z) * (B_x * C_y - B_y * C_x) * (C_y * A_z - C_z * A_y);
double term4 = -(A_z * B_x - A_x * B_z) * (B_y * C_z - B_z * C_y) * (C_x * A_y - C_y * A_x);
double term5 = +(A_x * B_y - A_y * B_x) * (B_y * C_z - B_z * C_y) * (C_z * A_x - C_x * A_z);
double term6 = -(A_x * B_y - A_y * B_x) * (B_z * C_x - B_x * C_z) * (C_y * A_z - C_z * A_y);
double P11 = term1 + term2 + term3 + term4 + term5 + term6;
std::cout << "Q = (A x B) . (B x C) x (C x A) = " << Q11 << std::endl;
std::cout << "P = (A x B) . (B x C) x (C x A) = " << P11 << std::endl;
}
static void Exercises_Section_1_14_Vector_Identities()
{
std::vector<double> A11 = { 1, 1, 1 };
std::vector<double> B11 = { 2, 4, -1 };
std::vector<double> C11 = { 1, 1, 3 };
std::vector<double> D11 = { 3, 2, 1 };
std::cout << "Section 1.14 page 50 Exercises Exercise 1" << std::endl;
std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
std::cout << "TPI1 = (A x B) x (C x D) = [A, B, D]C - [A, B, C]D = " << std::endl;
double TP1411a = TripleProduct(A11, B11, D11);
double TP1411b = TripleProduct(A11, B11, C11);
std::cout << "[A, B, D] = " << TP1411a << std::endl;
std::cout << "[A, B, C] = " << TP1411b << std::endl;
std::cout << "TPI1_x = " << TP1411a * C11[0] << std::endl;
std::cout << "TPI1_y = " << TP1411a * C11[1] << std::endl;
std::cout << "TPI1_z = " << TP1411a * C11[2] << std::endl;
std::cout << "TPI2_x = " << TP1411b * D11[0] << std::endl;
std::cout << "TPI2_y = " << TP1411b * D11[1] << std::endl;
std::cout << "TPI2_z = " << TP1411b * D11[2] << std::endl;
std::cout << "RHS1 = [A, B, D]C - [A, B, C]D = " << std::endl;
std::vector<double> RHS1(3);
RHS1[0] = TP1411a * C11[0] - TP1411b * D11[0];
RHS1[1] = TP1411a * C11[1] - TP1411b * D11[1];
RHS1[2] = TP1411a * C11[2] - TP1411b * D11[2];
std::cout << "RHS1_x = " << RHS1[0] << std::endl;
std::cout << "RHS1_y = " << RHS1[1] << std::endl;
std::cout << "RHS1_z = " << RHS1[2] << std::endl;
std::vector<double> CD11(3), TD11(3);
std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
VectorProduct(A11, B11, AB11);
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, BCCA11);
VectorProduct(A11, B11, AB11);
VectorProduct(C11, D11, CD11);
VectorProduct(AB11, CD11, TD11);
std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
std::cout << "A x B = " << AB11[0] << '\t' << AB11[1] << '\t' << AB11[2] << std::endl;
std::cout << "C x D = " << CD11[0] << '\t' << CD11[1] << '\t' << CD11[2] << std::endl;
std::cout << "TD11 = (A x B) x (C x D) = " << std::endl;
std::cout << "TD11_x = " << TD11[0] << std::endl;
std::cout << "TD11_y = " << TD11[1] << std::endl;
std::cout << "TD11_z = " << TD11[2] << std::endl;
VectorProduct(B11, C11, BC11);
VectorProduct(C11, A11, CA11);
VectorProduct(BC11, CA11, D11);
VectorProduct(A11, B11, AB11);
double ip12 = InnerProduct(AB11, D11, 3);
std::cout << "2. Inner Product = " << ip12 << std::endl;
double tp12 = TripleProduct(A11, B11, C11);
std::cout << "2. Triple Product ^ 2 = " << tp12 * tp12 << std::endl;
std::vector<double> ABC11(3), BAC11(3), CAB11(3);
VectorProduct(A11, BC11, ABC11);
VectorProduct(B11, CA11, BAC11);
VectorProduct(C11, AB11, CAB11);
double zx = ABC11[0] + BAC11[0] + CAB11[0];
double zy = ABC11[1] + BAC11[1] + CAB11[1];
double zz = ABC11[2] + BAC11[2] + CAB11[2];
std::cout << "3. Zero Vector = " << zx << ' ' << zy << ' ' << zz;
std::cout << std::endl;
}
int main()
{
Exercises_Section_1_13_1_Triple_Products();
Exercises_Section_1_14_Vector_Identities();
return 0;
}
Blog Entry © Tuesday, April 14, 2026, by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition © 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Page 48
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider
#include <iostream>
#include <vector>
static double InnerProduct(
std::vector<double> A,
std::vector<double> B,
int n)
{
double sum = 0.0;
for (int i = 0; i < n; i++)
sum += A[i] * B[i];
return sum;
}
static void VectorProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double>&C)
{
C.resize(3);
C[0] = A[1] * B[2] - A[2] * B[1];
C[1] = A[0] * B[2] - A[2] * B[0];
C[2] = A[0] * B[1] - A[1] * B[0];
}
static double TripleProduct(
std::vector<double> A,
std::vector<double> B,
std::vector<double> C)
{
double sum0 = 0.0, sum1 = 0.0;
sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
return sum0 - sum1;
}
static void Exercises_Section_1_13_1_Triple_Products()
{
// triple products (a)
std::vector<double> A1 = { 2, 0, 0 };
std::vector<double> B1 = { 0, 3, 0 };
std::vector<double> C1 = { 0, 0, 5 };
double tp_1 = TripleProduct(A1, B1, C1);
std::cout << "Exercise 1 (a)" << '\t';
std::cout << tp_1 << std::endl;
// triple products (b)
std::vector<double> A2 = { 1, 1, 1 };
std::vector<double> B2 = { 3, 1, 0 };
std::vector<double> C2 = { 0, -1, 5 };
std::cout << "Exercise 1 (b)" << '\t';
double tp_2 = TripleProduct(A2, B2, C2);
std::cout << tp_2 << std::endl;
// triple products (c)
std::vector<double> A3 = { 2, -1, 1 };
std::vector<double> B3 = { 1, 1, 1 };
std::vector<double> C3 = { 2, 0, 3 };
double tp_3 = TripleProduct(A3, B3, C3);
std::cout << "Exercise 1 (c)" << '\t';
std::cout << tp_3 << std::endl;
// triple products (d)
std::vector<double> A4 = { 0, 0, 1 };
std::vector<double> B4 = { 1, 0, 0 };
std::vector<double> C4 = { 0, 1, 0 };
double tp_4 = TripleProduct(A4, B4, C4);
std::cout << "Exercise 1 (d)" << '\t';
std::cout << tp_4 << std::endl;
// volume of a parallelpipped
std::vector<double> A5 = { 3, 4, 1 };
std::vector<double> B5 = { 2, 3, 4 };
std::vector<double> C5 = { 0, 0, 5 };
double tp_5 = TripleProduct(A5, B5, C5);
std::cout << "Exercise 2" << '\t';
std::cout << tp_5 << std::endl;
// volume of a parallelpipped
std::vector<double> A6 = { 3, 2, 1 };
std::vector<double> B6 = { 4, 2, 1 };
std::vector<double> C6 = { 0, 1, 4 };
std::vector<double> D6 = { 0, 0, 7 };
std::vector<double> AB = { -1, 0, 0 };
std::vector<double> AC = { 3, 1, -3 };
std::vector<double> AD = { 3, 2, -6 };
std::cout << "Exercise 3" << '\t';
double tp_6 = TripleProduct(AB, AC, AD);
std::cout << tp_6 << std::endl;
// volume of a tetrahedron
std::vector<double> AB1 = { 1, 1, 0 };
std::vector<double> AC1 = { 1, -1, 0 };
std::vector<double> AD1 = { 0, 0, 2 };
std::cout << "Exercise 4" << '\t';
double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
std::cout << fabs(tp_7) << std::endl;
std::vector<double> P1 = { 0, 0, 0 };
std::vector<double> P2 = { 1, 1, 0 };
std::vector<double> P3 = { 3, 4, 0 };
std::vector<double> P4 = { 4, 5, 0 };
std::vector<double> P5 = { 0, 0, 1 };
std::vector<double> Q1 = { 1, 1, 0 };
std::vector<double> Q2 = { 2, 3, 0 };
std::vector<double> Q3 = { 1, 1, 1 };
std::cout << "Exercise 5" << '\t';
double tp_8 = TripleProduct(Q1, Q2, Q3);
std::cout << fabs(tp_8) << std::endl;
std::vector<double> A10 = { 1, 1, 1 };
std::vector<double> B10 = { 2, 4, -1 };
std::vector<double> C10 = { 1, 1, 3 };
double tp_10 = TripleProduct(A10, B10, C10);
std::vector<double> D10 = { 0 };
VectorProduct(A10, B10, D10);
double magnitude = sqrt(InnerProduct(D10, D10, 3));
std::cout << "Exercise 10" << '\t';
std::cout << tp_10 / magnitude << '\t';
std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
}
int main()
{
Exercises_Section_1_13_1_Triple_Products();
return 0;
}
Blog Entry © Tuesday, April 7, 2026, by James Pate Williams, Jr., Hydrogen-like Atom Polar and Azimuthal Wavefunctions
Blog Entry © Saturday April 4, 2026, by James Pate Williams, Jr., Hydrogen-like Radial Electron Distribution Functions in CPP
#pragma once
class RadialWaveFunction
{
private:
static double Factorial(int n);
static double Laguerre(double rho, int n, int l);
public:
static double R(double r, int Z, int n, int l);
};
#include <math.h>
#include "RadialWaveFunction.h"
double RadialWaveFunction::Factorial(int n)
{
double factorial = 1.0;
for (int i = 2; i <= n; i++)
factorial *= i;
return factorial;
}
double RadialWaveFunction::Laguerre(double rho, int n, int l)
{
double sum = 0.0;
for (int k = 0; k <= n - l - 1; k++)
{
double factor1 = pow(-1, k + 2 * l + 1);
double factor2 = pow(Factorial(n + l), 2.0);
double factor3 = pow(rho, k);
double numer = factor1 * factor2 * factor3;
double factor4 = Factorial(n - l - 1 - k);
double factor5 = Factorial(2 * l + 1 + k);
double factor6 = Factorial(k);
double denom = factor4 * factor5 * factor6;
sum += numer / denom;
}
return sum;
}
double RadialWaveFunction::R(double r, int Z, int n, int l)
{
double rho = 2.0 * Z * r / n;
double numer1 = pow(2.0 * Z / n, 3.0);
double numer2 = Factorial(n - l - 1);
double denom1 = Factorial(n + n);
double denom2 = pow(Factorial(n + l), 3.0);
double numer3 = -sqrt(numer1 * numer2 /
(denom1 * denom2));
double exp2 = exp(-0.5 * rho);
double rhol = pow(rho, l);
return numer3 * exp2 * rhol * Laguerre(rho, n, l);
}
Blog Entry © Thursday, April 2, 2026, by James Pate Williams, Jr., Matrix Inverses that Exist and Their Characteristic Polynomials
#pragma once
#include <vector>
class Inverse
{
public:
static bool Compute(
int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X);
static void CharacteristicPolynomial(
int n, std::vector<double>& a,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& Madj);
};
// Reference: 5.2: The Characteristic Polynomial - Mathematics LibreTexts
#include "Inverse.h"
#include <vector>
bool Inverse::Compute(
int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X) {
double d = 0.0;
int i = 0, j = 0, k = 0, l = 0;
std::vector<double> C(n);
std::vector<std::vector<double>> B;
B.resize(n);
X.resize(n);
for (i = 0; i < n; i++)
{
B[i].resize(n);
X[i].resize(n);
for (j = 0; j < n; j++)
{
B[i][j] = 0.0;
X[i][j] = 0.0;
}
B[i][i] = 1.0;
}
j = -1;
Step2:
j++;
if (j == n)
goto Step6;
for (i = j; i < n; i++) {
if (M[i][j] != 0)
break;
}
if (i == n)
return false;
if (i > j) {
for (l = j; l < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
}
for (l = 0; l < n; l++) {
double temp = B[i][l];
B[i][l] = B[j][l];
B[j][l] = temp;
}
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
for (k = j + 1; k < n; k++) {
for (l = 0; l < n; l++) {
B[k][l] -= C[k] * B[j][l];
}
}
goto Step2;
Step6:
for (i = n - 1; i >= 0; i--) {
std::vector<double> sum(n, 0);
for (j = i + 1; j < n; j++) {
for (k = 0; k < n; k++) {
sum[k] += M[i][j] * X[j][k];
}
}
for (j = 0; j < n; j++)
X[i][j] = (B[i][j] - sum[j]) / M[i][i];
}
return true;
}
double delta(int i, int j) {
return i == j ? 1.0 : 0.0;
}
void Inverse::CharacteristicPolynomial(
int n, std::vector<double>& a,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& Madj) {
int i = 0, j = 0, k = 0, l = 0;
std::vector<double> X(n);
std::vector<std::vector<double>> MC;
a.resize(n + 1LL);
C.resize(n);
M.resize(n);
MC.resize(n);
Madj.resize(n);
for (i = 0; i < n; i++) {
C[i].resize(n);
M[i].resize(n);
MC[i].resize(n);
Madj[i].resize(n);
for (j = 0; j < n; j++) {
C[i][j] = delta(i, j);
}
}
a[0] = 1.0;
i = 0;
Step2:
i++;
if (i == n) {
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
double sum = 0.0;
for (l = 0; l < n; l++) {
sum += M[j][l] * C[l][k];
}
MC[j][k] = sum;
}
}
double trace = 0.0;
for (j = 0; j < n; j++) {
trace += MC[j][j];
}
a[n] = -trace / n;
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
Madj[j][k] = pow(-1, n - 1) * C[j][k];
}
}
return;
}
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
double sum = 0.0;
for (l = 0; l < n; l++) {
sum += M[j][l] * C[l][k];
}
MC[j][k] = sum;
}
}
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
C[j][k] = MC[j][k];
}
}
double tr = 0.0;
for (j = 0; j < n; j++) {
tr += C[j][j];
}
a[i] = -tr / i;
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
C[j][k] += a[i] * delta(j, k);
}
}
goto Step2;
}
Blog Entry © April 1, 2026, by James Pate Williams, Jr. More Romberg Integration Results
/*
Author: Pate Willimas c 1995
The following program is a translation of the FORTRAN77
subroutine found in "Elementary Numerical Analysis
An Algorithmic Approach Third Edition" (c) 1980 by S. D.
Conte and Carl de Boor pages 343-344. The program
uses Romberg extrapolation to find the integral of a
function. Also see: https://dlmf.nist.gov/3.5#E10
*/
#include <math.h>
#include <stdio.h>
#include <vector>
#include "Bessel.h"
typedef double real;
long fe;
real t[10][10];
real f(real x)
{
fe++;
return(expl(-x * x));
}
real g(real x)
{
fe++;
std::vector<real> j(4);
Bessel::bessj(x, 0, j);
return exp(-x) * j[0];
}
real romberg(
real a, real b,
int start, int row,
real (*fx)(real))
{
int i, /*j, */k, m;
real h, ratio, sum;
m = start;
h = (b - a) / m;
sum = 0.5 * (fx(a) + fx(b));
if (m > 1)
for (i = 1; i < m; i++)
sum += fx(a + i * h);
t[1][1] = h * sum;
//printf_s("Romberg T-Table\n%9.7lf\n", t[1][1]);
if (row < 2) return(t[1][1]);
for (k = 2; k <= row; k++)
{
h = 0.5 * h;
m *= 2;
sum = 0.0;
for (i = 1; i <= m; i += 2)
sum += fx(a + h * i);
t[k][1] = 0.5 * t[k - 1][1] + sum * h;
for (i = 1; i < k; i++)
{
t[k - 1][i] = t[k][i] - t[k - 1][i];
t[k][i + 1] = t[k][i] - t[k - 1][i] /
(powl(4.0, (real)i) - 1.0);
//for (j = 1; j < k; j++)
//printf_s("%9.7lf ", t[k][j]);
}
//printf_s("\n");
}
if (row < 3) return(t[2][2]);
//printf_s("Table of ratios\n");
for (k = 1; k <= row - 2; k++) {
for (i = 1; i <= k; i++)
{
if (t[k + 1][i] == 0.0)
ratio = 0.0;
else
ratio = t[k][i] / t[k + 1][i];
t[k][i] = ratio;
//printf_s("%5.2lf ", ratio);
}
//printf_s("\n");
}
return(t[row][row - 1]);
}
real CompositeTrapezoidalRule(
real a, real b, int n,
real(*fx)(real)) {
real endPts = 0.5 * (fx(a) + fx(b));
real sum = 0, xk = 0.0;
real h = (b - a) / n;
for (int k = 1; k <= n - 1; k++)
{
xk = a + k * h;
sum += fx(xk);
}
return h * (endPts + sum);
}
double CompositeSimpsonsRule(
double a, double b, int n,
real(*fx)(real))
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < n; i += 2)
{
s += fx(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < n; i += 2)
{
t += fx(x);
x += h2;
}
return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}
void Integrate(int trial, real(*fx)(real))
{
real a = 0.0, b = 0.0;
int row, start;
if (trial == 0)
{
a = 0.0;
b = 1.0;
}
else
{
a = 0.0;
b = 30.0;
}
if (trial == 0)
printf("Romberg integration of f(x) = exp(- x * x)\n");
else
printf("Romberg integration of g(x) = exp(-x) * J0(x)\n");
printf("number of trapezoidal intervals = ");
scanf_s("%d", &start);
printf("number of rows in table (<= 8) = ");
scanf_s("%d", &row);
fe = 0;
printf_s("integral = %14.11lf\tevals = %d\n",
romberg(a, b, start, row, fx), fe);
printf("number of trapezoidal intervals = ");
scanf_s("%d", &row);
fe = 0;
printf_s("integral = %14.11lf\tevals = %d\n",
CompositeTrapezoidalRule(a, b, row, fx), fe);
fe = 0;
printf("number of Simpson's intervals = ");
scanf_s("%d", &row);
printf_s("integral = %14.11lf\tevals = %d\n",
CompositeSimpsonsRule(a, b, row, fx), fe);
if (trial == 0)
printf_s("integral = %14.11lf\n", 0.74682413279);
else
{
real pi = 4.0 * atan(1.0);
real si = sin(pi * 45.0 / 180.0);
printf_s("integral = %14.11lf\n", si);
}
}
int main(void)
{
Integrate(0, f);
Integrate(1, g);
return(0);
}
Blog Entry © Sunday, March 29, 2026, James Pate Williams, Jr. Properties of Determinants
// Reference: "A Course in Computational
// Algebraic Number Theory by Henri Cohen
// Chapter 2 Algorithm 2.2.1 and Algorithm
// 2.2.3
#pragma once
#include <vector>
class Determinants
{
public:
static bool GaussianElimnation(
int n,
std::vector<std::vector<double>>& M,
std::vector<double> B,
std::vector<double> X);
static bool Determinant(
double& det, int n,
std::vector<std::vector<double>>& M);
};
#include "Determinants.h"
bool Determinants::GaussianElimnation(
int n,
std::vector<std::vector<double>>& M,
std::vector<double> B,
std::vector<double> X) {
double d = 0.0;
int i = 0, j = -1, k = 0, l = 0;
std::vector<double> C(n);
Step2:
j++;
if (j == n)
goto Step6;
for (i = j; i < n; i++) {
if (M[i][j] > 0)
break;
}
if (i == n)
return false;
if (i > j) {
for (l = j; j < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
temp = B[i];
B[i] = B[j];
B[j] = temp;
}
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
for (k = j + 1; k < n; k++) {
B[k] -= C[k] * B[j];
}
goto Step2;
Step6:
for (i = n - 1; i >= 0; i--) {
double sum = 0.0;
for (j = i + 1; j < n; j++) {
sum += M[i][j] * X[j];
}
X[i] = (B[i] - sum) / M[i][i];
}
return true;
}
bool Determinants::Determinant(
double& det, int n,
std::vector<std::vector<double>>& M) {
double d = 0.0, x = 1.0;
int i = 0, j = -1, k = 0, l = 0;
std::vector<double> C(n, 0);
Step2:
j++;
if (j == n) {
det = x;
return true;
}
for (i = j; i < n; i++) {
if (M[i][j] != 0.0)
break;
}
if (i == n) {
det = 0.0;
return false;
}
if (i > j) {
for (l = j; j < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
}
x = -x;
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
x *= M[j][j];
goto Step2;
}
// ComputeRealDeterminant.cpp : Defines the entry point for the application.
//
#include "framework.h"
#include "ComputeRealDeterminant.h"
#include "Determinants.h"
#include <string>
#include <vector>
#define MAX_LOADSTRING 100
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
WCHAR line[128], numWstr[128];
WCHAR inpStr[16384], outStr[16384];
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK DataInputDialog(HWND, UINT, WPARAM, LPARAM);
int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
_In_opt_ HINSTANCE hPrevInstance,
_In_ LPWSTR lpCmdLine,
_In_ int nCmdShow)
{
UNREFERENCED_PARAMETER(hPrevInstance);
UNREFERENCED_PARAMETER(lpCmdLine);
// TODO: Place code here.
// Initialize global strings
LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
LoadStringW(hInstance, IDC_COMPUTEREALDETERMINANT, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_COMPUTEREALDETERMINANT));
MSG msg;
// Main message loop:
while (GetMessage(&msg, nullptr, 0, 0))
{
if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
{
TranslateMessage(&msg);
DispatchMessage(&msg);
}
}
return (int) msg.wParam;
}
//
// FUNCTION: MyRegisterClass()
//
// PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW wcex = { 0 };
wcex.cbSize = sizeof(WNDCLASSEX);
wcex.style = CS_HREDRAW | CS_VREDRAW;
wcex.lpfnWndProc = WndProc;
wcex.cbClsExtra = 0;
wcex.cbWndExtra = 0;
wcex.hInstance = hInstance;
wcex.hIcon = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_COMPUTEREALDETERMINANT));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_COMPUTEREALDETERMINANT);
wcex.lpszClassName = szWindowClass;
wcex.hIconSm = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));
return RegisterClassExW(&wcex);
}
//
// FUNCTION: InitInstance(HINSTANCE, int)
//
// PURPOSE: Saves instance handle and creates main window
//
// COMMENTS:
//
// In this function, we save the instance handle in a global variable and
// create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
hInst = hInstance; // Store instance handle in our global variable
HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);
if (!hWnd)
{
return FALSE;
}
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
return TRUE;
}
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
switch (message)
{
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDM_START:
DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG), hWnd,
DataInputDialog);
break;
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
// TODO: Add any drawing code that uses hdc here...
EndPaint(hWnd, &ps);
}
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
INT_PTR CALLBACK DataInputDialog(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam) {
static int n = 0;
static std::wstring wstr;
static std::vector < std::vector<double>> A;
static HFONT hFont = nullptr;
static HWND hEditMultiline1 = nullptr;
static HWND hEditMultiline2 = nullptr;
UNREFERENCED_PARAMETER(lParam);
switch (message) {
case WM_INITDIALOG:
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_SETCURSEL, 0, 0);
hFont = CreateFont(
-MulDiv(10, GetDeviceCaps(GetDC(hDlg), LOGPIXELSY), 72),
0, 0, 0, FW_BOLD, FALSE, FALSE, FALSE,
DEFAULT_CHARSET, OUT_DEFAULT_PRECIS,
CLIP_DEFAULT_PRECIS, DEFAULT_QUALITY,
FIXED_PITCH | FF_MODERN,
TEXT("Courier New")
);
hEditMultiline1 = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | WS_VSCROLL |
ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL,
10, 100, 250, 250, // Position and size
hDlg, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE1, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
hEditMultiline2 = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | ES_LEFT | ES_READONLY,
120, 360, 120, 30, // Position and size
hDlg, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE2, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
break;
case WM_COMMAND:
if (LOWORD(wParam) == IDC_BUTTON_COMPUTE) {
GetDlgItemText(hDlg, IDC_COMBO_N, inpStr, 16384);
WCHAR* endptr = nullptr;
n = (int)wcstol(inpStr, &endptr, 10);
A.resize(n);
GetDlgItemText(hDlg, IDC_EDIT_MULTILINE1, inpStr, 16384);
int k = 0;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
A[i].resize(n);
int l = 0;
while (inpStr[k] != L' ') {
numWstr[l++] = inpStr[k++];
}
k++;
numWstr[l++] = L'\0';
A[i][j] = (double)wcstof(numWstr, &endptr);
}
}
double det = 0.0;
bool result = Determinants::Determinant(det, n, A);
swprintf_s(line, 128, L"%lf", det);
SetDlgItemText(hDlg, IDC_EDIT_MULTILINE2, line);
break;
}
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL) {
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
//{{NO_DEPENDENCIES}}
// Microsoft Visual C++ generated include file.
// Used by ComputeRealDeterminant.rc
#define IDS_APP_TITLE 103
#define IDR_MAINFRAME 128
#define IDD_COMPUTEREALDETERMINANT_DIALOG 102
#define IDD_ABOUTBOX 103
#define IDM_ABOUT 104
#define IDM_EXIT 105
#define IDI_COMPUTEREALDETERMINANT 107
#define IDI_SMALL 108
#define IDC_COMPUTEREALDETERMINANT 109
#define IDC_MYICON 2
#ifndef IDC_STATIC
#define IDC_STATIC -1
#endif
// Next default values for new objects
//
#ifdef APSTUDIO_INVOKED
#ifndef APSTUDIO_READONLY_SYMBOLS
#define _APS_NO_MFC 130
#define _APS_NEXT_RESOURCE_VALUE 129
#define _APS_NEXT_COMMAND_VALUE 32771
#define _APS_NEXT_CONTROL_VALUE 1000
#define _APS_NEXT_SYMED_VALUE 110
#endif
#endif
#define IDD_DATA_INPUT_DIALOG 1000
#define IDC_COMBO_N 1010
#define IDC_EDIT_MULTILINE1 1020
#define IDC_EDIT_MULTILINE2 1030
#define IDC_BUTTON_COMPUTE 1040
#define IDM_START 32771
//Microsoft Visual C++ generated resource script.
//
#include "resource.h"
#define APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
#ifndef APSTUDIO_INVOKED
#include "targetver.h"
#endif
#define APSTUDIO_HIDDEN_SYMBOLS
#include "windows.h"
#undef APSTUDIO_HIDDEN_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
#undef APSTUDIO_READONLY_SYMBOLS
#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU)
LANGUAGE 9, 1
/////////////////////////////////////////////////////////////////////////////
//
// Icon
//
// Icon with lowest ID value placed first to ensure application icon
// remains consistent on all systems.
IDI_COMPUTEREALDETERMINANT ICON "ComputeRealDeterminant.ico"
IDI_SMALL ICON "small.ico"
/////////////////////////////////////////////////////////////////////////////
//
// Menu
//
IDC_COMPUTEREALDETERMINANT MENU
BEGIN
POPUP "&Begin"
BEGIN
MENUITEM "&Start", IDM_START
MENUITEM SEPARATOR
MENUITEM "E&xit", IDM_EXIT
END
POPUP "&Help"
BEGIN
MENUITEM "&About ...", IDM_ABOUT
END
END
/////////////////////////////////////////////////////////////////////////////
//
// Accelerator
//
IDC_COMPUTEREALDETERMINANT ACCELERATORS
BEGIN
"?", IDM_ABOUT, ASCII, ALT
"/", IDM_ABOUT, ASCII, ALT
END
/////////////////////////////////////////////////////////////////////////////
//
// Dialog
//
IDD_ABOUTBOX DIALOGEX 0, 0, 170, 62
STYLE DS_SETFONT | DS_MODALFRAME | DS_FIXEDSYS | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "About ComputeRealDeterminant"
FONT 8, "MS Shell Dlg"
BEGIN
ICON IDI_COMPUTEREALDETERMINANT,IDC_STATIC,14,14,21,20
LTEXT "ComputeRealDeterminant, Version 1.0",IDC_STATIC,42,14,114,8,SS_NOPREFIX
LTEXT "Copyright (c) 2026",IDC_STATIC,42,26,114,8
DEFPUSHBUTTON "OK",IDOK,113,41,50,14,WS_GROUP
END
IDD_DATA_INPUT_DIALOG DIALOGEX 0, 0, 280, 260
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog"
FONT 10, "Courier New", 700
BEGIN
LTEXT "n:", IDC_STATIC, 10, 10, 30, 10
COMBOBOX IDC_COMBO_N, 70, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "input matrix:", IDC_STATIC, 10, 30, 200, 10
LTEXT "determinant:", IDC_STATIC, 10, 185, 200, 10
PUSHBUTTON "Compute", IDC_BUTTON_COMPUTE, 10, 220, 50, 16
PUSHBUTTON "Cancel", IDCANCEL, 220, 220, 50, 16
END
/////////////////////////////////////////////////////////////////////////////
//
// DESIGNINFO
//
#ifdef APSTUDIO_INVOKED
GUIDELINES DESIGNINFO
BEGIN
IDD_ABOUTBOX, DIALOG
BEGIN
LEFTMARGIN, 7
RIGHTMARGIN, 163
TOPMARGIN, 7
BOTTOMMARGIN, 55
END
END
#endif // APSTUDIO_INVOKED
#ifdef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// TEXTINCLUDE
//
1 TEXTINCLUDE
BEGIN
"resource.h\0"
END
2 TEXTINCLUDE
BEGIN
"#ifndef APSTUDIO_INVOKED\r\n"
"#include ""targetver.h""\r\n"
"#endif\r\n"
"#define APSTUDIO_HIDDEN_SYMBOLS\r\n"
"#include ""windows.h""\r\n"
"#undef APSTUDIO_HIDDEN_SYMBOLS\r\n"
"\0"
END
3 TEXTINCLUDE
BEGIN
"\r\n"
"\0"
END
#endif // APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// String Table
//
STRINGTABLE
BEGIN
IDC_COMPUTEREALDETERMINANT "COMPUTEREALDETERMINANT"
IDS_APP_TITLE "ComputeRealDeterminant"
END
#endif
/////////////////////////////////////////////////////////////////////////////
#ifndef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
/////////////////////////////////////////////////////////////////////////////
#endif // not APSTUDIO_INVOKED