Category: Legendre Functions
Blog Entry © Tuesday, March 24, 2026, by James Pate Williams, Jr. Cubature on a Disk
// CubatureDisk.cpp (c) Tuesday, March 24, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <iomanip>
#include <iostream>
#include <vector>
#include "Mueller.h"
static void SelectionSort(
std::vector<double>& x,
std::vector<double>& w,
int n) {
for (int i = 0; i < n - 1; i++) {
for (int j = i + 1; j < n; j++) {
if (x[i] > x[j]) {
double t = x[i];
x[i] = x[j];
x[j] = t;
t = w[i];
w[i] = w[j];
w[j] = t;
}
}
}
}
static double fr(double rho, int i, int j) {
return fabs(rho) * pow(rho, i + j);
}
static double ft(double t, int i, int j) {
double t1 = 1.0 - t * t;
return (1.0 / sqrt(t1)) * pow(t1, i / 2.0) * pow(t, j);
}
static double GaussLegendre(
std::vector<double> x, std::vector<double> w,
double (*fx)(double, int, int), int i, int j, int n) {
double integral = 0.0;
for (int k = 0; k < n; k++) {
integral += fx(x[k], i, j) * w[k];
}
return integral;
}
static double Integrate(
double (*fx)(double, int, int),
int i, int j, int n,
unsigned int seed) {
std::string msg;
std::vector<double> x(n);
std::vector<double> weights(n);
Mueller::OrthogonalPolynomialRoots(
"Legendre", "Zeros", n, n, seed, msg, x);
for (int i = 0; i < n; i++) {
double xi = x[i];
double dpn = Mueller::g("Legendre", "Extrema",
xi, n, 0, x);
weights[i] = 2.0 / ((1.0 - xi * xi) * dpn * dpn);
}
SelectionSort(x, weights, n);
return GaussLegendre(x, weights, fx, i, j, n);
}
static double SimpsonsRule(
int i, int j, int n, double a, double b,
double (*fx)(double, int, int))
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int k = 1; k < n; k += 2)
{
s += fx(x, i, j);
x += h2;
}
x = a + h2;
for (int k = 2; k < n; k += 2)
{
t += fx(x, i, j);
x += h2;
}
return h * (fx(a, i, j) + 4 * s + 2 * t + fx(b, i, j)) / 3.0;
}
double MonteCarlo(int i, int j, int N) {
double integral = 0.0;
for (int k = 0; k < N; k++)
{
double x = (double)rand() / RAND_MAX;
double y = (double)rand() / RAND_MAX;
if (x * x + y * y <= 1)
integral += pow(x, i) * pow(y, j);
}
return 4.0 * integral / N;
}
static double Integrand(std::vector<double> x, double n) {
double power = 1.0;
for (int k = 0; k < n; k++) {
power *= sqrt(fabs(x[k] - 1.0 / 3.0));
}
return power;
}
static double MonteCarloPower(int N, int n) {
double integration = 0.0;
std::vector<double> x(n);
for (int i = 0; i < N; i++) {
for (int k = 0; k < n; k++) {
x[k] = (double)rand() / RAND_MAX;
}
integration += Integrand(x, n);
}
return integration / N;
}
static double MonteCarloPi(int N) {
double integral = 0.0;
for (int i = 0; i < N; i++) {
double x = (double)rand() / RAND_MAX;
double y = (double)rand() / RAND_MAX;
if (sqrt(x * x + y * y) <= 1)
integral++;
}
return 4.0 * integral / N;
}
int main() {
int i = 4, j = 4, n = 0, N = 10000000;
unsigned int seed = 1;
std::cout << "Integral # 1" << std::endl;
for (int n = 5; n <= 25; n += 5)
{
double integral1 = Integrate(fr, i, j, n, seed++);
double integral2 = Integrate(ft, i, j, n, seed++);
double integral = integral1 * integral2;
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << n << "\t " << integral << std::endl;
}
for (int n = 50000000; n <= 90000000; n += 10000000)
{
double integral1 = SimpsonsRule(i, j, n, -1 + 1.0e-15, 1 - 1.0e-15, fr);
double integral2 = SimpsonsRule(i, j, n, -1 + 1.0e-15, 1 - 1.0e-15, ft);
double integral = integral1 * integral2;
std::cout << "Simpson ";
if (n < 10000000)
std::cout << n << "\t\t";
else
std::cout << n << '\t';
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral << std::endl;
}
double integral3 = MonteCarlo(i, j, N);
std::cout << "Monte Carlo ";
std::cout << N << "\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral3 << std::endl;
N = 10000000;
n = 6;
double integral4 = MonteCarloPower(N, n);
double actual4 = pow(2.0 / 3.0, n) * pow(pow(2.0 / 3.0, 1.5) +
pow(1.0 / 3.0, 1.5), n);
std::cout << "Integral # 2" << std::endl;
std::cout << "Monte Carlo ";
std::cout << N << "\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral4 << std::endl;
std::cout << "actual\t\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << actual4 << std::endl;
double perror = 100.0 * fabs(integral4 - actual4) /
fabs(actual4);
std::cout << "% error\t\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << perror << std::endl;
double actual5 = 4.0 * atan(1.0);
double integral5 = MonteCarloPi(N);
std::cout << "Integral # 3" << std::endl;
std::cout << "Monte Carlo ";
std::cout << N << "\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral5 << std::endl;
std::cout << "actual\t\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << actual5 << std::endl;
perror = 100.0 * fabs(integral5 - actual5) /
fabs(actual5);
std::cout << "% error\t\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << perror << std::endl;
return 0;
}
#pragma once
/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Numerical Analysis: An Algorithmic
* Approach" by S. D. Conte and Carl de Boor
* Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1
* Translated to better C++ on Saturday, March 21,
* 2026
*/
#include <string>
#include <vector>
class Mueller
{
private:
static double epsilon1, epsilon2;
public:
static double f(std::string name,
double x, int n, int m,
std::vector<double>& zeros);
static double g(
std::string name,
std::string graph,
double x, int n, int m,
std::vector<double>& zeros);
static double OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros);
static void OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros);
};
#include "Mueller.h"
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
double Mueller::epsilon1 = 1.0e-12;
double Mueller::epsilon2 = 1.0e-15;
double Mueller::f(std::string name,
double x, int n, int m,
std::vector<double>& zeros)
{
double a1n = 1, a2n = 1, a3n = 1, a4n = 1;
double f0 = 1, f1 = x, f2 = 1;
if (name == "Laguerre")
f1 = -x + 1;
for (int i = 2; i <= n; i++)
{
if (name == "Chebyshev")
{
a1n = 1.0;
a2n = 0.0;
a3n = 2.0;
a4n = 1.0;
}
else if (name == "Laguerre")
{
a1n = i - 1 + 1;
a2n = 2 * (i - 1) + 1;
a3n = -1.0;
a4n = i - 1;
}
else if (name == "Legendre")
{
a1n = i + 1 - 1;
a2n = 0.0;
a3n = 2.0 * (i - 1) + 1;
a4n = i - 1;
}
f2 = ((a2n + a3n * x) * f1 - a4n * f0) / a1n;
f0 = f1;
f1 = f2;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) < epsilon2)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
}
return f2;
}
double Mueller::g(std::string name,
std::string graph, double x,
int n, int m, std::vector<double>& zeros)
{
double f0 = 1, f1 = x, f2 = 1;
if (graph != "Extrema")
return f(name, x, n, m, zeros);
else
{
if (name == "Laguerre") {
if (x == 1.0)
x -= 0.0000001;
else if (x == -1.0)
x += 0.0000001;
}
double g0x = 1.0, g1x = 1.0, g2x = 1.0;
if (name == "Chebyshev") {
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
else if (name == "Laguerre") {
g0x = -n;
g1x = n;
g2x = x;
}
else if (name == "Legendre")
{
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
if (fabs(g2x) < epsilon2)
{
std::cout << "g2x is too small!";
return 0.0;
}
f0 = f(name, x, n - 1, 0, zeros);
f1 = f(name, x, n, 0, zeros);
f2 = (g1x * f1 + g0x * f0) / g2x;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) == 0.0)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
return f2;
}
}
double Mueller::OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros)
{
double x3 = 0.0;
int i = 2;
while (true)
{
double h1 = x1 - x0;
if (fabs(h1) < epsilon2)
{
msg = "h1 is too small!";
return 0.0;
}
double h2 = x2 - x1;
if (fabs(h2) < epsilon2)
{
msg = "h2 is too small!";
return 0.0;
}
double f0 = g(name, graph, x0, n, m, zeros);
double f1 = g(name, graph, x1, n, m, zeros);
double f2 = g(name, graph, x2, n, m, zeros);
double g1 = (f1 - f0) / h1;
double g2 = (f2 - f1) / h2;
if (fabs(f2) < epsilon2)
{
x3 = x2;
break;
}
if (fabs(h1 + h2) < epsilon2)
{
msg = "h1 + h2 is too small!";
return 0.0;
}
double k2 = (g2 - g1) / (h1 + h2);
double c1 = g2 + h2 * k2;
double d = c1 * c1 - 4.0 * f2 * k2;
if (d < 0.0)
d = 0.0;
double s1 = sqrt(d);
double m1 = c1 - s1;
double p1 = c1 + s1;
double denom;
if (fabs(m1) > fabs(p1))
denom = m1;
else
denom = p1;
if (fabs(denom) < epsilon2)
{
msg = "denom is too small!";
return 0.0;
}
double h3 = -2.0 * f2 / denom;
x3 = x2 + h3;
double f3 = g(name, graph, x3, n, m, zeros);
if (fabs(f3) < epsilon2)
break;
if (fabs(h3) < epsilon2)
{
msg = "h3 is too small!";
return 0.0;
}
double g3 = (f3 - f2) / h3;
double delta = fabs(x3 - x2);
if (delta < epsilon1)
break;
i++;
if (i > maxIt)
break;
x0 = x1;
x1 = x2;
x2 = x3;
if (fabs(x1 - x0) < epsilon1)
break;
if (fabs(x2 - x1) < epsilon1)
break;
}
return x3;
}
void Mueller::OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros) {
double epsilon1 = 1.0e-10;
epsilon2 = 1.0e-20;
double x0 = -0.5, x1 = 0.0, x2 = 0.5;
double z;
int maxIt = 128;
zeros.resize(nr);
srand(seed);
for (int m = 0; m < nr; m++) {
x0 = (double)rand() / RAND_MAX;
x1 = (double)rand() / RAND_MAX;
x2 = (double)rand() / RAND_MAX;
z = zeros[m] = OrthogonalPolynomialRoot(
name, graph, n, epsilon1, epsilon2,
maxIt, m, x0, x1, x2, msg, zeros);
}
}
Blog Entry © Monday, March 23, 2026, by James Pate Williams, Jr. Gauss-Legendre Quadrature
#pragma once
/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Numerical Analysis: An Algorithmic
* Approach" by S. D. Conte and Carl de Boor
* Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1
* Translated to better C++ on Saturday, March 21,
* 2026
*/
#include <string>
#include <vector>
class Mueller
{
private:
static double epsilon1, epsilon2;
public:
static double f(std::string name,
double x, int n, int m,
std::vector<double>& zeros);
static double g(
std::string name,
std::string graph,
double x, int n, int m,
std::vector<double>& zeros);
static double OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros);
static void OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros);
};
#include "Mueller.h"
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
double Mueller::epsilon1 = 1.0e-12;
double Mueller::epsilon2 = 1.0e-15;
double Mueller::f(std::string name,
double x, int n, int m,
std::vector<double>& zeros)
{
double a1n = 1, a2n = 1, a3n = 1, a4n = 1;
double f0 = 1, f1 = x, f2 = 1;
if (name == "Laguerre")
f1 = -x + 1;
for (int i = 2; i <= n; i++)
{
if (name == "Chebyshev")
{
a1n = 1.0;
a2n = 0.0;
a3n = 2.0;
a4n = 1.0;
}
else if (name == "Laguerre")
{
a1n = i - 1 + 1;
a2n = 2 * (i - 1) + 1;
a3n = -1.0;
a4n = i - 1;
}
else if (name == "Legendre")
{
a1n = i + 1 - 1;
a2n = 0.0;
a3n = 2.0 * (i - 1) + 1;
a4n = i - 1;
}
f2 = ((a2n + a3n * x) * f1 - a4n * f0) / a1n;
f0 = f1;
f1 = f2;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) < epsilon2)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
}
return f2;
}
double Mueller::g(std::string name,
std::string graph, double x,
int n, int m, std::vector<double>& zeros)
{
double f0 = 1, f1 = x, f2 = 1;
if (graph != "Extrema")
return f(name, x, n, m, zeros);
else
{
if (name == "Laguerre") {
if (x == 1.0)
x -= 0.0000001;
else if (x == -1.0)
x += 0.0000001;
}
double g0x = 1.0, g1x = 1.0, g2x = 1.0;
if (name == "Chebyshev") {
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
else if (name == "Laguerre") {
g0x = -n;
g1x = n;
g2x = x;
}
else if (name == "Legendre")
{
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
if (fabs(g2x) < epsilon2)
{
std::cout << "g2x is too small!";
return 0.0;
}
f0 = f(name, x, n - 1, 0, zeros);
f1 = f(name, x, n, 0, zeros);
f2 = (g1x * f1 + g0x * f0) / g2x;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) == 0.0)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
return f2;
}
}
double Mueller::OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros)
{
double x3 = 0.0;
int i = 2;
while (true)
{
double h1 = x1 - x0;
if (fabs(h1) < epsilon2)
{
msg = "h1 is too small!";
return 0.0;
}
double h2 = x2 - x1;
if (fabs(h2) < epsilon2)
{
msg = "h2 is too small!";
return 0.0;
}
double f0 = g(name, graph, x0, n, m, zeros);
double f1 = g(name, graph, x1, n, m, zeros);
double f2 = g(name, graph, x2, n, m, zeros);
double g1 = (f1 - f0) / h1;
double g2 = (f2 - f1) / h2;
if (fabs(f2) < epsilon2)
{
x3 = x2;
break;
}
if (fabs(h1 + h2) < epsilon2)
{
msg = "h1 + h2 is too small!";
return 0.0;
}
double k2 = (g2 - g1) / (h1 + h2);
double c1 = g2 + h2 * k2;
double d = c1 * c1 - 4.0 * f2 * k2;
if (d < 0.0)
d = 0.0;
double s1 = sqrt(d);
double m1 = c1 - s1;
double p1 = c1 + s1;
double denom;
if (fabs(m1) > fabs(p1))
denom = m1;
else
denom = p1;
if (fabs(denom) < epsilon2)
{
msg = "denom is too small!";
return 0.0;
}
double h3 = -2.0 * f2 / denom;
x3 = x2 + h3;
double f3 = g(name, graph, x3, n, m, zeros);
if (fabs(f3) < epsilon2)
break;
if (fabs(h3) < epsilon2)
{
msg = "h3 is too small!";
return 0.0;
}
double g3 = (f3 - f2) / h3;
double delta = fabs(x3 - x2);
if (delta < epsilon1)
break;
i++;
if (i > maxIt)
break;
x0 = x1;
x1 = x2;
x2 = x3;
if (fabs(x1 - x0) < epsilon1)
break;
if (fabs(x2 - x1) < epsilon1)
break;
}
return x3;
}
void Mueller::OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros) {
double epsilon1 = 1.0e-10;
epsilon2 = 1.0e-20;
double x0 = -0.5, x1 = 0.0, x2 = 0.5;
double z;
int maxIt = 128;
zeros.resize(nr);
srand(seed);
for (int m = 0; m < nr; m++) {
x0 = (double)rand() / RAND_MAX;
x1 = (double)rand() / RAND_MAX;
x2 = (double)rand() / RAND_MAX;
z = zeros[m] = OrthogonalPolynomialRoot(
name, graph, n, epsilon1, epsilon2,
maxIt, m, x0, x1, x2, msg, zeros);
}
}
// Quadrature1d.cpp (c) Saturday, March 21, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Reference: https://dlmf.nist.gov/3.5#x
#include <iomanip>
#include <iostream>
#include <vector>
#include "Mueller.h"
static void SelectionSort(
std::vector<double>& x,
std::vector<double>& w,
int n) {
for (int i = 0; i < n - 1; i++) {
for (int j = i + 1; j < n; j++) {
if (x[i] > x[j]) {
double t = x[i];
x[i] = x[j];
x[j] = t;
t = w[i];
w[i] = w[j];
w[j] = t;
}
}
}
}
static double fx(double x) {
return x * x * x * exp(-x);
}
static double GaussLegendre(
std::vector<double> x, std::vector<double> w,
double (*fx)(double), int n) {
double integral = 0.0;
for (int i = 0; i < n; i++) {
integral += fx(x[i]) * w[i];
}
return integral;
}
static double Integrate(double (*fx)(double), int n,
unsigned int seed) {
std::string msg;
std::vector<double> x(n);
std::vector<double> weights(n);
Mueller::OrthogonalPolynomialRoots(
"Legendre", "Zeros", n, n, seed, msg, x);
for (int i = 0; i < n; i++) {
double xi = x[i];
double dpn = Mueller::g("Legendre", "Extrema",
xi, n, 0, x);
weights[i] = 2.0 / ((1.0 - xi * xi) * dpn * dpn);
}
SelectionSort(x, weights, n);
std::cout << "Abscissas (x) and Weights (w)" << std::endl;
std::cout << "for Gauss-Legendre Quadrature";
std::cout << std::endl;
std::cout << "x\t\t\tw" << std::endl;
for (int i = 0; i < n; i++) {
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << x[i] << '\t';
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << weights[i] << std::endl;
}
return GaussLegendre(x, weights, fx, n);
}
static double SimpsonsRule(
int n, double a, double b,
double (*fx)(double))
{
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;
}
int main() {
unsigned int seed = 1;
for (int n = 10; n <= 40; n += 10)
{
double integral = Integrate(fx, n, seed++);
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << n << '\t' << integral << std::endl;
}
for (int n = 50000; n <= 100000; n += 10000)
{
double integral = SimpsonsRule(n, -1, 1, fx);
std::cout << "Simpson ";
std::cout << n << "\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral << std::endl;
}
return 0;
}
Blog Entry © Monday, November 10, 2025, by James Pate Williams, Jr. More COMP 640 Advanced Computer Graphics Results
Blog Entry © Sunday, November 9, 2025, by James Pate Williams, Jr. Hydrogenic Wavefunctions, Radial Probability Functions, Distribution Functions, and First Moment Integrals
Four Methods of Numerical Double Integration: Sequential Simpson’s Rule, Multitasking Simpson’s Rule, Sequential Monte Carlo and Multitasking Monte Carlo Methods © Wednesday April 16 – 18, 2025, by James Pate Williams, Jr.
Approximation of the Ground-State Total Energy of a Beryllium Atom © Sunday, March 30 to Tuesday April 1, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Computer Science
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;
}
Numerical Explorations 