Category: Numerical Analysis

Determinant Calculations: Bareiss Algorithm and Gaussian Elimination by James Pate Williams, Jr.

Below are three screenshots of two methods of calculating the determinant of a matrix, namely the Bareiss Algorithm and Gaussian Elimination:

using System;
using System.Diagnostics;
using System.Windows.Forms;

namespace MatrixInverseComparison
{
    public partial class MainForm : Form
    {
        public MainForm()
        {
            InitializeComponent();
        }
        static private string FormatNumber(double x)
        {
            string result = string.Empty;

            if (x > 0)
                result += x.ToString("F5").PadLeft(10);
            else
                result += x.ToString("F5").PadLeft(10);

            return result;
        }

        private void MultiplyPrintMatricies(
            double[,] A, double[,] B, int n)
        {
            double[,] I = new double[n, n];

            textBox1.Text += "Matrix Product:\r\n";

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    double sum = 0.0;

                    for (int k = 0; k < n; k++)
                        sum += A[i, k] * B[k, j];

                    I[i, j] = sum;
                }
            }

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    textBox1.Text += FormatNumber(I[i, j]) + " ";
                }

                textBox1.Text += "\r\n";
            }
        }

        private void PrintMatrix(string title, double[,] A, int n)
        {
            textBox1.Text += title + "\r\n";

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    textBox1.Text += FormatNumber(A[i, j]) + " ";
                }

                textBox1.Text += "\r\n";
            }
        }

        private void Compute(double[,] MI, int n)
        {
            double determinantGE = 1;
            double[] b = new double[n];
            double[] x = new double[n];
            double[,] MB1 = new double[n, n];
            double[,] MB2 = new double[n, n];
            double[,] MG = new double[n, n];
            double[,] MS = new double[n, n];
            double[,] IG = new double[n, n];
            double[,] IS = new double[n, n];
            int[] pivot = new int[n];
            Stopwatch sw = new();
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    MB1[i, j] = MB2[i, j] = MG[i, j] = MS[i, j] = MI[i, j];
            }
            PrintMatrix("Initial Matrix: ", MI, n);
            sw.Start();
            int flag = DirectMethods.Factor(MG, n, pivot);
            if (flag != 0)
            {
                for (int i = 0; i < n; i++)
                    determinantGE *= MG[i, i];
                determinantGE *= flag;
            }
            sw.Stop();
            PrintMatrix("Gaussian Elimination Final:", MG, n);
            for (int i = 0; i < n; i++)
                b[i] = 0;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    b[j] = 1.0;
                    DirectMethods.Substitute(MG, b, x, n, pivot);
                    for (int k = 0; k < n; k++)
                        IG[k, j] = x[k];
                    b[j] = 0.0;
                }
            }
            PrintMatrix("Gaussian Elimination Inverse:", IG, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinantGE) + "\r\n";
            MultiplyPrintMatricies(IG, MI, n);
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
            sw.Start();
            double determinant1 = BareissAlgorithm.Determinant1(MB1, n);
            sw.Stop();
            PrintMatrix("Bareiss Algorithm Final 1:", MB1, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinant1) + "\r\n";
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
            sw.Start();
            double determinant2 = BareissAlgorithm.Determinant2(MB2, n);
            sw.Stop();
            PrintMatrix("Bareiss Algorithm Final 2:", MB2, n);
            textBox1.Text += "Determinant: " +
                FormatNumber(determinant2) + "\r\n";
            textBox1.Text += "Runtime (MS) = " +
                sw.ElapsedMilliseconds + "\r\n";
        }
        private void button1_Click(object sender, EventArgs e)
        {
            int n = (int)numericUpDown1.Value;
            int seed = (int)numericUpDown2.Value;
            double[,] A = new double[n, n];
            Random random = new Random(seed);

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    double q = n * random.NextDouble();

                    while (q == 0.0)
                        q = n * random.NextDouble();

                    A[i, j] = q;
                }
            }
            
            Compute(A, n);
        }

        private void button2_Click(object sender, EventArgs e)
        {
            textBox1.Text = string.Empty;
        }
    }
}

using System;

namespace MatrixInverseComparison
{
    public class DirectMethods
    {
        // Substitute and Factor translated from FORTRAN 77
        // source code found in "Elementary Numerical Analysis:
        // An Algorithmic Approach" by S. D. Conte and Carl de
        // Boor. Translator: James Pate Williams, Jr. (c)
        // August 14 - 17, 2023
        static public void Substitute(
            double[,] w,
            double[] b,
            double[] x,
            int n,
            int[] pivot)
        {
            double sum;
            int i, j, n1 = n - 1;

            if (n == 1)
            {
                x[0] = b[0] / w[0, 0];
                return;
            }

            // forward substitution
            x[0] = b[pivot[0]];

            for (i = 0; i < n; i++)
            {
                for (j = 0, sum = 0.0; j < i; j++)
                    sum += w[i, j] * x[j];
                x[i] = b[pivot[i]] - sum;
            }

            // backward substitution
            x[n1] /= w[n1, n1];

            for (i = n - 2; i >= 0; i--)
            {
                for (j = i + 1, sum = 0.0; j < n; j++)
                    sum += w[i, j] * x[j];
                x[i] = (x[i] - sum) / w[i, i];
            }
        }

        // Factor returns +1 if an even number of exchanges
        // Factor returns -1 if an odd  number of exchanges
        // Factor retrurn 0  if matrix is singular
        static public int Factor(
            double[,] w, int n, int[] pivot)
        // returns 0 if matrix is singular
        {
            double awikod, col_max, ratio, row_max, temp;
            double[] d = new double[n];
            int flag = 1, i, i_star, j, k;

            for (i = 0; i < n; i++)
            {
                pivot[i] = i;
                row_max = 0;
                for (j = 0; j < n; j++)
                    row_max = Math.Max(row_max, Math.Abs(w[i, j]));
                if (row_max == 0)
                {
                    flag = 0;
                    row_max = 1;
                }
                d[i] = row_max;
            }
            if (n <= 1)
                return flag;
            // factorization
            for (k = 0; k < n - 1; k++)
            {
                // determine pivot row the row i_star
                col_max = Math.Abs(w[k, k]) / d[k];
                i_star = k;
                for (i = k + 1; i < n; i++)
                {
                    awikod = Math.Abs(w[i, k]) / d[i];
                    if (awikod > col_max)
                    {
                        col_max = awikod;
                        i_star = i;
                    }
                }
                if (col_max == 0)
                    flag = 0;
                else
                {
                    if (i_star > k)
                    {
                        // make k the pivot row by
                        // interchanging with i_star
                        flag *= -1;
                        i = pivot[i_star];
                        pivot[i_star] = pivot[k];
                        pivot[k] = i;
                        temp = d[i_star];
                        d[i_star] = d[k];
                        d[k] = temp;
                        for (j = 0; j < n; j++)
                        {
                            temp = w[i_star, j];
                            w[i_star, j] = w[k, j];
                            w[k, j] = temp;
                        }
                    }
                    // eliminate x[k]
                    for (i = k + 1; i < n; i++)
                    {
                        w[i, k] /= w[k, k];
                        ratio = w[i, k];
                        for (j = k + 1; j < n; j++)
                            w[i, j] -= ratio * w[k, j];
                    }
                }
            }

            if (w[n - 1, n - 1] == 0)
                flag = 0;

            return flag;
        }
    }
}

namespace MatrixInverseComparison
{
    // One implementation is based on https://en.wikipedia.org/wiki/Bareiss_algorithm
    // Another perhaps better implementation is found on the following webpage
    // https://cs.stackexchange.com/questions/124759/determinant-calculation-bareiss-vs-gauss-algorithm
    class BareissAlgorithm
    {
        static public double Determinant1(double[,] M, int n)
        {
            double M00;

            for (int k = 0; k < n; k++)
            {
                if (k - 1 == -1)
                    M00 = 1;
                else
                    M00 = M[k - 1, k - 1];

                for (int i = k + 1; i < n; i++)
                {
                    for (int j = k + 1; j < n; j++)
                    {
                        M[i, j] = (
                            M[i, j] * M[k, k] -
                            M[i, k] * M[k, j]) / M00;
                    }
                }
            }

            return M[n - 1, n - 1];
        }

        static public double Determinant2(double[,] A, int dim)
        {
            if (dim <= 0)
            {
                return 0;
            }

            double sign = 1;

            for (int k = 0; k < dim - 1; k++)
            {
                //Pivot - row swap needed
                if (A[k, k] == 0)
                {
                    int m;
                    for (m = k + 1; m < dim; m++)
                    {
                        if (A[m, k] != 0)
                        {
                            double[] tempRow = new double[dim];

                            for (int i = 0; i < dim; i++)
                                tempRow[i] = A[m, i];
                            for (int i = 0; i < dim; i++)
                                A[m, i] = A[k, i];
                            for (int i = 0; i < dim; i++)
                                A[k, i] = tempRow[i];
                            sign = -sign;
                            break;
                        }
                    }
                    //No entries != 0 found in column k -> det = 0
                    if (m == dim)
                    {
                        return 0;
                    }
                }
                //Apply formula
                for (int i = k + 1; i < dim; i++)
                {
                    for (int j = k + 1; j < dim; j++)
                    {
                        A[i, j] = A[k, k] * A[i, j] - A[i, k] * A[k, j];

                        if (k >= 1)
                        {
                            A[i, j] /= A[k - 1, k - 1];
                        }
                    }
                }
            }

            return sign * A[dim - 1, dim - 1];
        }
    }
}

Blast from the Past 1997 Image of a Matrix by James Pate Williams, Jr.

/*
  Author:  Pate Williams (c) 1997

  "Algorithm 2.3.2 (Image of a Matrix). Given an
  m by n matrix M = (m[i][i]) with 1 <= i <= m and
  1 <= j <= n having coefficients in the field K,
  this algorithm outputs a basis of the image of
  M, i. e. vector space spanned by the columns of
  M." -Henri Cohen- See "A Course in Computational
  Algebraic Number Theory" by Henri Cohen pages
  58-59. We specialize the code to the field Zp.
*/

#include <stdio.h>
#include <stdlib.h>

long** create_matrix(long m, long n)
{
    long i, ** matrix = (long**)calloc(m, sizeof(long*));

    if (!matrix) {
        fprintf(stderr, "fatal error\ninsufficient memory\n");
        fprintf(stderr, "from create_matrix\n");
        exit(1);
    }
    for (i = 0; i < m; i++) {
        matrix[i] = (long*)calloc(n, sizeof(long));
        if (!matrix[i]) {
            fprintf(stderr, "fatal error\ninsufficient memory\n");
            fprintf(stderr, "from create_matrix\n");
            exit(1);
        }
    }
    return matrix;
}

void delete_matrix(long m, long** matrix)
{
    long i;

    for (i = 0; i < m; i++) free(matrix[i]);
    free(matrix);
}

void Euclid_extended(long a, long b, long* u,
    long* v, long* d)
{
    long q, t1, t3, v1, v3;

    *u = 1, * d = a;
    if (b == 0) {
        *v = 0;
        return;
    }
    v1 = 0, v3 = b;
#ifdef DEBUG
    printf("----------------------------------\n");
    printf(" q    t3   *u   *d   t1   v1   v3\n");
    printf("----------------------------------\n");
#endif
    while (v3 != 0) {
        q = *d / v3;
        t3 = *d - q * v3;
        t1 = *u - q * v1;
        *u = v1, * d = v3;
#ifdef DEBUG
        printf("%4ld %4ld %4ld ", q, t3, *u);
        printf("%4ld %4ld %4ld %4ld\n", *d, t1, v1, v3);
#endif
        v1 = t1, v3 = t3;
    }
    *v = (*d - a * *u) / b;
#ifdef DEBUG
    printf("----------------------------------\n");
#endif
}

long inv(long number, long modulus)
{
    long d, u, v;

    Euclid_extended(number, modulus, &u, &v, &d);
    if (d == 1) return u;
    return 0;
}

void image(long m, long n, long p,
    long** M, long** X, long* r)
{
    int found;
    long D, i, j, k, s;
    long* c = (long*)calloc(m, sizeof(long));
    long* d = (long*)calloc(n, sizeof(long));
    long** N = create_matrix(m, n);

    if (!c || !d) {
        fprintf(stderr, "fatal error\ninsufficient memory\n");
        fprintf(stderr, "from kernel\n");
        exit(1);
    }
    for (i = 0; i < m; i++) {
        c[i] = -1;
        for (j = 0; j < n; j++) N[i][j] = M[i][j];
    }
    *r = 0;
    for (k = 0; k < n; k++) {
        found = 0, j = 0;
        while (!found && j < m) {
            found = M[j][k] != 0 && c[j] == -1;
            if (!found) j++;
        }
        if (found) {
            D = p - inv(M[j][k], p);
            M[j][k] = p - 1;
            for (s = k + 1; s < n; s++)
                M[j][s] = (D * M[j][s]) % p;
            for (i = 0; i < m; i++) {
                if (i != j) {
                    D = M[i][k];
                    M[i][k] = 0;
                    for (s = k + 1; s < n; s++) {
                        M[i][s] = (M[i][s] + D * M[j][s]) % p;
                        if (M[i][s] < 0) M[i][s] += p;
                    }
                }
            }
            c[j] = k;
            d[k] = j;
        }
        else {
            *r = *r + 1;
            d[k] = -1;
        }
    }
    for (j = 0; j < m; j++) {
        if (c[j] != -1) {
            for (i = 0; i < n; i++) {
                if (i < m) X[i][j] = N[i][c[j]];
                else X[i][j] = 0;
            }
        }
    }
    delete_matrix(m, N);
    free(c);
    free(d);
}

void print_matrix(long m, long n, long** a)
{
    long i, j;

    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++)
            printf("%2ld ", a[i][j]);
        printf("\n");
    }
}

int main(void)
{
    long i, j, m = 8, n = 8, p = 13, r;
    long a[8][8] = { {0,  0,  0,  0,  0,  0,  0,  0},
                    {2,  0,  7, 11, 10, 12,  5, 11},
                    {3,  6,  3,  3,  0,  4,  7,  2},
                    {4,  3,  6,  4,  1,  6,  2,  3},
                    {2, 11,  8,  8,  2,  1,  3, 11},
                    {6, 11,  8,  6,  2,  6, 10,  9},
                    {5, 11,  7, 10,  0, 11,  6, 12},
                    {3,  3, 12,  5,  0, 11,  9, 11} };
    long** M = create_matrix(m, n);
    long** X = create_matrix(n, n);

    for (i = 0; i < m; i++)
        for (j = 0; j < n; j++)
            M[i][j] = a[j][i];
    printf("the original matrix is as follows:\n");
    print_matrix(m, n, M);
    image(m, n, p, M, X, &r);
    printf("the image of the matrix is as follows:\n");
    print_matrix(n, n - r, X);
    printf("the rank of the matrix is: %ld\n", n - r);
    delete_matrix(m, M);
    delete_matrix(n, X);
    return 0;
}

Part of a New Linear Algebra Package in C++ Implemented by James Pate Williams, Jr.

// Algorithms from "A Course in Computational
// Algebraic Number Theory" by Henri Cohen
// Implemented by James Pate Williams, Jr.
// Copyright (c) 2023 All Rights Reserved

#pragma once
#include "pch.h"

template<class T> class Matrix
{
public:
	size_t m, n;
	T** data;

	Matrix() { m = 0; n = 0; data = NULL; };
	Matrix(size_t m, size_t n)
	{
		this->m = m;
		this->n = n;
		data = new T*[m];

		if (data == NULL)
			exit(-300);

		for (size_t i = 0; i < m; i++)
		{
			data[i] = new T[n];

			if (data[i] == NULL)
				exit(-301);
		}
	};
	void OutputMatrix(
		fstream& outs, char fill, int precision, int width)
	{
		for (size_t i = 0; i < m; i++)
		{
			for (size_t j = 0; j < n; j++)
			{
				outs << setfill(fill) << setprecision(precision);
				outs << setw(width) << data[i][j] << '\t';
			}

			outs << endl;
		}
	};
};

template<class T> class Vector
{
public:
	size_t n;
	T* data;

	Vector() { n = 0; data = NULL; };
	Vector(size_t n)
	{
		this->n = n;
		data = new T[n];
	};
	void OutputVector(
		fstream& outs, char fill, int precision, int width)
	{
		for (size_t i = 0; i < n; i++)
		{
			outs << setfill(fill) << setprecision(precision);
			outs << setw(width) << data[i] << '\t';
		}

		outs << endl;
	};
};

class LinearAlgebra
{
public:
	bool initialized;
	size_t m, n;
	Matrix<double> M;
	Vector<double> B;

	LinearAlgebra() { 
		initialized = false;
		m = 0; n = 0;
		M.data = NULL;
		B.data = NULL;
	};
	LinearAlgebra(size_t m, size_t n) {
		initialized = false;
		this->m = m;
		this->n = n;
		this->M.m = m;
		this->M.n = n;
		this->B.n = n;
		this->M.data = new double*[m];
		this->B.data = new double[n];

		if (M.data != NULL)
		{
			for (size_t i = 0; i < m; i++)
			{
				this->M.data[i] = new double[n];

				for (size_t j = 0; j < n; j++)
					this->M.data[i][j] = 0;
			}
		}

		if (B.data != NULL)
		{
			this->B.data = new double[n];

			for (size_t i = 0; i < n; i++)
				this->B.data[i] = 0;
		}

		initialized = this->B.data != NULL && this->M.data != NULL;
	};
	LinearAlgebra(
		size_t m, size_t n,
		double** M,
		double* B)
	{
		this->m = m;
		this->n = n;
		this->M.m = m;
		this->M.n = n;
		this->M.data = new double*[m];

		if (M != NULL)
		{
			for (size_t i = 0; i < m; i++)
			{
				this->M.data[i] = new double[n];

				for (size_t j = 0; j < n; j++)
					this->M.data[i][j] = M[i][j];
			}
		}

		if (B != NULL)
		{
			this->B.data = new double[n];

			for (size_t i = 0; i < m; i++)
				this->B.data[i] = B[i];
		}

		initialized = this->B.data != NULL && this->M.data != NULL;
	}
	~LinearAlgebra()
	{
		M.m = m;
		M.n = n;
		B.n = n;

		if (B.data != NULL)
			delete[] B.data;

		for (size_t i = 0; i < m; i++)
			if (M.data[i] != NULL)
				delete[] M.data[i];

		if (M.data != NULL)
			delete[] M.data;
	}
	double DblDeterminant(size_t n, bool failure);
	Vector<double> DblGaussianElimination(
		bool& failure);
	// The following three methods are from the
	// textbook "Elementary Numerical Analysis
	// An Algorithmic Approach" by S. D. Conte
	// and Carl de Boor Translated from the
	// original FORTRAN by James Pate Williams, Jr.
	// Copyright (c) 2023 All Rights Reserved
	bool DblGaussianFactor(
		size_t n,
		Vector<int>& pivot);
	bool DblGaussianSolution(
		int n,
		Vector<double>& x,
		Vector<int>& pivot);
	bool DblSubstitution(
		size_t n,
		Vector<double>& x,
		Vector<int>& pivot);
	bool DblInverse(
		size_t n,
		Matrix<double>& A,
		Vector<int>& pivot);
};

#include "pch.h"
#include "LinearAlgebra.h"

double LinearAlgebra::DblDeterminant(
    size_t n, bool failure)
{
    double deter = 1;
    Vector<int> pivot(n);

    if (!initialized || m != n)
    {
        failure = true;
        return 0.0;
    }

    if (!DblGaussianFactor(n, pivot))
    {
        failure = true;
        return 0.0;
    }

    for (size_t i = 0; i < n; i++)
        deter *= M.data[i][i];

    return deter;
}

Vector<double> LinearAlgebra::DblGaussianElimination(
    bool& failure)
{
    double* C = new double[m];
    Vector<double> X(n);

    X.data = new double[n];

    if (X.data == NULL)
        exit(-200);

    if (!initialized)
    {
        failure = true;
        delete[] C;
        return X;
    }
    
    for (size_t i = 0; i < m; i++)
        C[i] = -1;

    failure = false;

    for (size_t j = 0; j < n; j++)
    {
        bool found = false;
        size_t i = j;

        while (i < n && !found)
        {
            if (M.data[i][j] != 0)
                found = true;
            else
                i++;
        }

        if (!found)
        {
            failure = true;
            break;
        }

        if (i > j)
        {
            for (size_t l = j; l < n; l++)
            {
                double t = M.data[i][l];
                M.data[i][l] = M.data[j][l];
                M.data[j][l] = t;
                t = B.data[i];
                B.data[i] = B.data[j];
                B.data[j] = t;
            }
        }

        double d = 1.0 / M.data[j][j];

        for (size_t k = j + 1; k < n; k++)
            C[k] = d * M.data[k][j];

        for (size_t k = j + 1; k < n; k++)
        {
            for (size_t l = j + 1; l < n; l++)
                M.data[k][l] = M.data[k][l] - C[k] * M.data[j][l];
            
            B.data[k] = B.data[k] - C[k] * B.data[j];
        }
    }

    for (int i = (int)n - 1; i >= 0; i--)
    {
        double sum = 0;

        for (size_t j = i + 1; j < n; j++)
            sum += M.data[i][j] * X.data[j];

        X.data[i] = (B.data[i] - sum) / M.data[i][i];
    }

    delete[] C;
    return X;
}

bool LinearAlgebra::DblGaussianFactor(
    size_t n,
    Vector<int>& pivot)
    // returns false if matrix is singular
{
    Vector<double> d(n);
    double awikod, col_max, ratio, row_max, temp;
    int flag = 1;
    size_t i_star, itemp;

    for (size_t i = 0; i < n; i++)
    {
        pivot.data[i] = i;
        row_max = 0;
        for (size_t j = 0; j < n; j++)
            row_max = max(row_max, abs(M.data[i][j]));
        if (row_max == 0)
        {
            flag = 0;
            row_max = 1;
        }
        d.data[i] = row_max;
    }
    if (n <= 1) return flag != 0;
    // factorization
    for (size_t k = 0; k < n - 1; k++)
    {
        // determine pivot row the row i_star
        col_max = abs(M.data[k][k]) / d.data[k];
        i_star = k;
        for (size_t i = k + 1; i < n; i++)
        {
            awikod = abs(M.data[i][k]) / d.data[i];
            if (awikod > col_max)
            {
                col_max = awikod;
                i_star = i;
            }
        }
        if (col_max == 0)
            flag = 0;
        else
        {
            if (i_star > k)
            {
                // make k the pivot row by
                // interchanging with i_star
                flag *= -1;
                itemp = pivot.data[i_star];
                pivot.data[i_star] = pivot.data[k];
                pivot.data[k] = itemp;
                temp = d.data[i_star];
                d.data[i_star] = d.data[k];
                d.data[k] = temp;
                for (size_t j = 0; j < n; j++)
                {
                    temp = M.data[i_star][j];
                    M.data[i_star][j] = M.data[k][j];
                    M.data[k][j] = temp;
                }
            }
            // eliminate x[k]
            for (size_t i = k + 1; i < n; i++)
            {
                M.data[i][k] /= M.data[k][k];
                ratio = M.data[i][k];
                for (size_t j = k + 1; j < n; j++)
                    M.data[i][j] -= ratio * M.data[k][j];
            }
        }

        if (M.data[n - 1][n - 1] == 0) flag = 0;
    }

    if (flag == 0)
        return false;

    return true;
}

bool LinearAlgebra::DblGaussianSolution(
    int n,
    Vector<double>& x,
    Vector<int>& pivot)
{
    if (!DblGaussianFactor(n, pivot))
        return false;

    return DblSubstitution(n, x, pivot);
}

bool LinearAlgebra::DblSubstitution(
    size_t n, Vector<double>& x,
    Vector<int>& pivot)
{
    double sum;
    size_t j, n1 = n - 1;

    if (n == 1)
    {
        x.data[0] = B.data[0] / M.data[0][0];
        return true;
    }

    // forward substitution

    x.data[0] = B.data[pivot.data[0]];

    for (int i = 1; i < (int)n; i++)
    {
        for (j = 0, sum = 0; j < (size_t)i; j++)
            sum += M.data[i][j] * x.data[j];

        x.data[i] = B.data[pivot.data[i]] - sum;
    }

    // backward substitution

    x.data[n1] /= M.data[n1][n1];

    for (int i = n - 2; i >= 0; i--)
    {
        double sum = 0.0;

        for (j = i + 1; j < n; j++)
            sum += M.data[i][j] * x.data[j];

        x.data[i] = (x.data[i] - sum) / M.data[i][i];
    }

    return true;
}

bool LinearAlgebra::DblInverse(
    size_t n,
    Matrix<double>& A,
    Vector<int>& pivot)
{
    Vector<double> x(n);

    if (!DblGaussianFactor(n, pivot))
        return false;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            B.data[j] = 0;
    }
    
    for (size_t i = 0; i < n; i++)
    {
        B.data[i] = 1;

        if (!DblSubstitution(n, x, pivot))
            return false;

        B.data[i] = 0;

        for (size_t j = 0; j < n; j++)
           A.data[i][j] = x.data[pivot.data[j]];
    }

    return true;
}

/*
** Cohen's linear algebra test program
** Implemented by James Pate Williams, Jr.
** Copyright (c) 2023 All Rights Reserved
*/

#include "pch.h"
#include "LinearAlgebra.h"

double GetDblNumber(fstream& inps)
{
    char ch = inps.get();
    string numberStr;

    while (ch == ' ' || ch == '\t' || ch == '\r' || ch == '\n')
        ch = inps.get();

    while (ch == '+' || ch == '-' || ch == '.' ||
        ch >= '0' && ch <= '9')
    {
        numberStr += ch;
        ch = inps.get();
    }

    double x = atof(numberStr.c_str());
    return x;
}

int GetIntNumber(fstream& inps)
{
    char ch = inps.get();
    string numberStr;

    while (ch == ' ' || ch == '\t' || ch == '\r' || ch == '\n')
        ch = inps.get();

    while (ch == '+' || ch == '-' || ch >= '0' && ch <= '9')
    {
        numberStr += ch;
        ch = inps.get();
    }

    int x = atoi(numberStr.c_str());
    return x;
}

int main()
{
    fstream inps;

    inps.open("CLATestFile.txt", fstream::in);
    
    if (inps.fail())
    {
        cout << "Input file opening error!" << endl;
        return -1;
    }

    fstream outs;

    outs.open("CLAResuFile.txt", fstream::out | fstream::ate);

    if (outs.fail())
    {
        cout << "Output file opening error!" << endl;
        return -2;
    }

    size_t m, n;
    
    while (!inps.eof())
    {
        m = GetIntNumber(inps);

        if (inps.eof())
            return 0;

        if (m < 1)
        {
            cout << "The number of rows must be >= 1" << endl;
            return -100;
        }

        n = GetIntNumber(inps);

        if (n < 1)
        {
            cout << "The number of rows must be >= 1" << endl;
            return -101;
        }

        LinearAlgebra la(m, n);
        Matrix<double> copyM(m, n);
        Vector<double> copyB(n);

        for (size_t i = 0; i < m; i++)
        {
            for (size_t j = 0; j < n; j++)
            {
                double x = GetDblNumber(inps);

                la.M.data[i][j] = x;
                copyM.data[i][j] = x;
            }
        }

        for (size_t i = 0; i < n; i++)
        {
           la.B.data[i] = GetDblNumber(inps);
           copyB.data[i] = la.B.data[i];
        }

        bool failure = false;
        Vector<double> X = la.DblGaussianElimination(failure);

        if (!failure)
            X.OutputVector(outs, ' ', 5, 8);
        else
        {
            cout << "Cohen Gaussian elimination failure!" << endl;
            exit(-102);
        }


        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        Matrix<double> A(n, n);
        Vector<int> pivot(n);
        
        if (!la.DblGaussianSolution(n, X, pivot))
            exit(-103);

        X.OutputVector(outs, ' ', 5, 8);

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        double deter = la.DblDeterminant(n, failure);

        outs << deter << endl;

        for (size_t i = 0; i < m; i++)
        {
            la.B.data[i] = copyB.data[i];

            for (size_t j = 0; j < n; j++)
            {
                la.M.data[i][j] = copyM.data[i][j];
            }
        }

        if (!la.DblInverse(n, A, pivot))
        {
            cout << "Conte Gaussian inverse matrix failure!" << endl;
            exit(-104);
        }

        else
            A.OutputMatrix(outs, ' ', 5, 8);
    }

    inps.close();
    outs.close();
}

2
2
1	1
1	2
7	11
2
2
1	1
1	3
7	11
2
2
6	3
4	8
5	6
2
2
5	3
10	4
8	6
3
3
2	1	-1
-3	-1	2
-2	1	2
8	-11	-3

       3	       4	
       3	       4	
1
       2	      -1	
      -1	       1	
       5	       2	
       5	       2	
2
     1.5	    -0.5	
    -0.5	     0.5	
 0.61111	 0.44444	
 0.61111	 0.44444	
36
 0.22222	-0.11111	
-0.083333	 0.16667	
    -1.4	       5	
    -1.4	       5	
-10
    -0.4	       1	
     0.3	    -0.5	
       2	       3	      -1	
       2	       3	      -1	
1
       4	       5	      -2	
       3	       4	      -2	
      -1	      -1	       1

Multiple Integration Using Simpson’s Rule – Calculation of the Ground State Energies of the First Five Non-Relativistic Multiple Electron Atoms by James Pate Williams, Jr.

We calculate the ground state energies for Helium (Z = 2), Lithium (Z = 3), Beryllium (Z = 4), Boron (Z = 5), and Carbon (Z = 6). Currently, only three of the five atoms are implemented.

heliumatom-printoutDownload
lithiumatom-printoutDownload
berylliumatom-printoutDownload
integration-printoutDownload
ground-states-of-first-five-atomsDownload

Multiple Integration Using Simpson’s Rule – Calculation of the Ground State of the Non-Relativistic Helium Atom by James Pate Williams, Jr.

The six-dimensional Cartesian coordinate wavefunction is calculated by the C# method:

public double Psi(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double exp1 = Math.Exp(-alpha[0] * r1);
            double exp2 = Math.Exp(-alpha[1] * r2);
            double exp3 = Math.Exp(-alpha[2] * r12);

            return exp1 * exp2 * exp3;
        }

        public double Psi2(double[] x, double[] alpha)
        {
            double psi = Psi(x, alpha);

            return psi * psi;
        }

The wavefunction normalization method is:

public double Normalize(double[] alpha, int nSteps)
        {
            double[] lower = new double[6];
            double[] upper = new double[6];
            int[] steps = new int[6];

            lower[0] = 0.001;
            lower[1] = 0.001;
            lower[2] = 0.001;
            lower[3] = 0.001;
            lower[4] = 0.001;
            lower[5] = 0.001;

            upper[0] = 10.0;
            upper[1] = 10.0;
            upper[2] = 10.0;
            upper[3] = 10.0;
            upper[4] = 10.0;
            upper[5] = 10.0;

            for (int i = 0; i < 6; i++)
                steps[i] = nSteps;

            double norm = Math.Sqrt(integ.Integrate(
                lower, upper, alpha, Psi2, 6, steps));

            return norm;
        }

The two kinetic energy integrands are encapsulated in the following two methods:

private double Integrand1(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double term = -2.0 * alpha[0] / r1 + alpha[0] * alpha[0];
            double mul1 = Math.Exp(-2.0 * alpha[0] * r1);
            double mul2 = Math.Exp(-2.0 * alpha[1] * r2);
            double mul3 = Math.Exp(-2.0 * alpha[2] * r12);

            return N * N * term * mul1 * mul2 * mul3;
        }

        private double Integrand2(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double term = -2.0 * alpha[1] / r1 + alpha[1] * alpha[1];
            double mul1 = Math.Exp(-2.0 * alpha[0] * r1);
            double mul2 = Math.Exp(-2.0 * alpha[1] * r2);
            double mul3 = Math.Exp(-2.0 * alpha[2] * r12);

            return N * N * term * mul1 * mul2 * mul3;
        }

The two potential energy terms use the following two integrands:

private double Integrand3(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double term = 1.0 / r1;
            double mul =
                Math.Exp(-2.0 * alpha[0] * r1) *
                Math.Exp(-2.0 * alpha[1] * r2) *
                Math.Exp(-2.0 * alpha[2] * r12);

            return -N * N * Z * term * mul;
        }

        private double Integrand4(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double term = 1.0 / r2;
            double mul =
                Math.Exp(-2.0 * alpha[0] * r1) *
                Math.Exp(-2.0 * alpha[1] * r2) *
                Math.Exp(-2.0 * alpha[2] * r12);

            return -N * N * Z * term * mul;
        }

The electron-electron repulsion integrand is given by:

private double Integrand5(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));

            if (r12 == 0)
                r12 = 0.01;

            double term = 1.0 / r12;
            double mul =
                Math.Exp(-2.0 * alpha[0] * r1) *
                Math.Exp(-2.0 * alpha[1] * r2) *
                Math.Exp(-2.0 * alpha[2] * r12);

            return N * N * term * mul;
        }

The ground state non-relativistic energy is computed by the method:

public double Energy(double[] alpha, int nSteps, int Z)
        {
            double[] lower = new double[6];
            double[] upper = new double[6];
            int[] steps = new int[6];

            lower[0] = lower[1] = lower[2] = lower[3] = lower[4] = lower[5] = 0.001;
            upper[0] = upper[1] = upper[2] = upper[3] = upper[4] = upper[5] = 10.0;
            steps[0] = steps[1] = steps[2] = steps[3] = steps[4] = steps[5] = nSteps;

            N = Normalize(alpha, nSteps);

            this.Z = Z;

            double integ1 = integ.Integrate(lower, upper, alpha, Integrand1, 6, steps);
            double integ2 = integ.Integrate(lower, upper, alpha, Integrand2, 6, steps);
            double integ3 = integ.Integrate(lower, upper, alpha, Integrand3, 6, steps);
            double integ4 = integ.Integrate(lower, upper, alpha, Integrand4, 6, steps);
            double integ5 = integ.Integrate(lower, upper, alpha, Integrand5, 6, steps);
            
            return integ1 + integ2 - integ3 - integ4 + integ5;
        }

Using trial and error we calculate Alpha as: 0.535139999999

public double Integrate(double[] lower, double[] upper, double[] alpha,
            Func<double[], double[], double> f, int n, int[] steps)
        {
            double p = 1;
            double[] h = new double[n];
            double[] h2 = new double[n];
            double[] s = new double[n];
            double[] t = new double[n];
            double[] x = new double[n];
            double[] w = new double[n];

            for (int i = 0; i < n; i++)
            {
                h[i] = (upper[i] - lower[i]) / steps[i];
                h2[i] = 2.0 * h[i];
            }

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    x[j] = lower[j];
                
                x[i] = lower[i] + h[i];

                for (int j = 1; j < steps[i]; j += 2)
                {
                    s[i] += f(x, alpha);
                    x[i] += h2[i];
                }

                for (int j = 0; j < n; j++)
                    x[j] = lower[j];

                x[i] = lower[i] + h2[i];

                for (int j = 2; j < steps[i]; j += 2)
                {
                    t[i] += f(x, alpha);
                    x[i] += h2[i];
                }

                w[i] = h[i] * (f(lower, alpha) + 4 * s[i] + 2 * t[i] + f(upper, alpha)) / 3.0;
            }

            for (int i = 0; i < n; i++)
                p *= w[i];

            return p;
        }

Romberg Integration of the Logarithmic Integral by James Pate Williams, Jr.

/*
Author:  Pate Williams (c) January 22, 1995
The following program is a translation of the FORTRAN
subprogram found in Elementary Numerical Analysis by
S. D. Conte and Carl de Boor pages 343-344. The program
uses Romberg extrapolation to find the integral of a
function.
*/

#include "stdafx.h"
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>

typedef long double real;

real t[10][10];

real f(real x)
{
	return(expl(-x * x));
}

real Li(real x)
{
	return(1.0 / logl(x));
}

real romberg(
	real a, real b, int start, int row,
	real(*f)(real))
{
	int i, k, m;
	real h, ratio, sum;

	m = start;
	h = (b - a) / m;
	sum = 0.5 * (f(a) + f(b));
	if (m > 1)
		for (i = 1; i < m; i++)
			sum += f(a + i * h);
	t[1][1] = h * sum;
	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 += f(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);
		}
	}
	if (row < 3) return(t[2][2]);
	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;
		}
	return(t[row][row - 1]);
}

int main(void)
{
	long double experimental = 0.7468241328124271;
	int row = 8;

	printf("first test\n");
	printf("Romberg integration of f(x) = exp(- x * x)\n");
	printf("Internet value = 0.7468241328124271\n");
	printf("from website = https://www.integral-calculator.com/\n");
	printf("Approximation\t\tPercent Difference\tSteps\tRuntime (s)\n");
	for (long steps = 100000; steps <= 900000; steps += 100000)
	{
		clock_t start = clock();
		long double trial = romberg(0.0, 1.0, steps, row, f);
		clock_t endTm = clock();
		long double runtime = ((long double)endTm - start) /
			CLOCKS_PER_SEC;
		long double pd = 100.0 * fabsl(experimental - trial) /
			(0.5 *(experimental + trial));
		printf("%16.15lf\t%16.15le\t%6d\t%4.3lf\n", 
			trial, pd, steps, runtime);
	}
	printf("second test\n");
	printf("Romberg integration of Li(x) = 1 / ln x between 2 and some n\n");
	printf("approximates the number of prime numbers between 2 and n\n");
	printf("n\t\tLi(n)\t\tSteps\tTime (s)\n");
	for (long n = 10; n <= 1000000000; n *= 10)
	{
		long steps = 1000000;

		clock_t start = clock();
		long double trial = romberg(2.0, n, steps, row, Li);
		clock_t endTm = clock();
		long double runtime = ((long double)endTm - start) /
			CLOCKS_PER_SEC;
		printf("%12ld\t%12.0lf\t%7ld\t%4.3lf\n",
			n, trial, steps, runtime);
	}
	return(0);
}

Set Implementation of Sieve of Eratosthenes by James Pate Williams, Jr.

/*
Author:  Pate Williams (c) January 20, 1995
The following is a translation of the Pascal program
sieve found in Pascalgorithms by Edwin D. Reilly and
Francis D. Federighi page 652. This program uses sets
to represent the sieve (see C Programming Language An
Applied Perspective by Lawrence Miller and Alec Qui-
lici pages 160 - 162).
*/
#include "stdafx.h"
#include <math.h>
#include <stdio.h>
#include <time.h>

#define _WORD_SIZE 32
#define _VECT_SIZE 31250000
#define SET_MIN    0
#define SET_MAX    1000000000

typedef long LONG;
typedef long SET[_VECT_SIZE];
typedef LONG ELEMENT;

SET set;

static LONG get_bit_pos(LONG *long_ptr, LONG *bit_ptr,
	ELEMENT element)
{
	*long_ptr = element / _WORD_SIZE;
	*bit_ptr = element % _WORD_SIZE;
	return(element >= SET_MIN && element <= SET_MAX);
}

static void set_bit(ELEMENT element, LONG inset)
{
	LONG bit, word;

	if (get_bit_pos(&word, &bit, element))
		inset ? set[word] |= (01 << bit) :
		set[word] &= ~(01 << bit);
}

static int get_bit(ELEMENT element)
{
	LONG bit, word;

	return(get_bit_pos(&word, &bit, element) ?
		(set[word] >> bit) & 01 : 0);
}

void set_Add(ELEMENT element)
{
	set_bit(element, 1);
}

void set_Del(ELEMENT element)
{
	set_bit(element, 0);
}

int set_Mem(ELEMENT element)
{
	return get_bit(element);
}

void primes(LONG n)
{
	LONG c, i, inc, k;
	double x;

	clock_t now = clock();
	set_Add(2);
	for (i = 3; i <= n; i++)
		if ((i + 1) % 2 == 0)
			set_Add(i);
		else
			set_Del(i);
	c = 3;
	do
	{
		i = c * c;
		inc = c + c;
		while (i <= n)
		{
			set_Del(i);
			i = i + inc;
		}
		c += 2;
		while (set_Mem(c) == 0) c += 1;
	} while (c * c <= n);
	k = 0;
	for (i = 2; i <= n; i++)
		if (set_Mem(i) == 1) k++;
	x = n / log(n) - 5.0;
	x = x + exp(1.0 + 0.15 * log(n) * sqrt(log(n)));
	clock_t later = clock();
	double runtime = (double)(later - now) / CLOCKS_PER_SEC;
	printf("%10ld\t%10ld\t%10.0lf\t%6.4lf\n",
		n, k, x, runtime);
}

int main(void)
{
	LONG n = 10L;

	printf("--------------------------------------------------------\n");
	printf("n\t\tprimes\t\ttheory\t\ttime (s)\n");
	printf("--------------------------------------------------------\n");
	do
	{
		primes(n);
		clock_t later = clock();
		n = 10L * n;
	} while (n < 1000000000);
	printf("--------------------------------------------------------\n");
	return(0);
}

Romberg Extrapolation by James Pate Williams, Jr.

/*
Author:  Pate Williams (c) January 22, 1995
The following program is a translation of the FORTRAN
subprogram found in Elementary Numerical Analysis by
S. D. Conte and Carl de Boor pages 343-344. The program
uses Romberg extrapolation to find the integral of a
function.
*/

#include "stdafx.h"
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>

typedef long double real;

real t[10][10];

real f(real x)
{
	return(expl(-x * x));
}

real romberg(real a, real b, int start, int row)
{
	int i, k, m;
	real h, ratio, sum;

	m = start;
	h = (b - a) / m;
	sum = 0.5 * (f(a) + f(b));
	if (m > 1)
		for (i = 1; i < m; i++)
			sum += f(a + i * h);
	t[1][1] = h * sum;
	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 += f(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);
		}
	}
	if (row < 3) return(t[2][2]);
	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;
		}
	return(t[row][row - 1]);
}

int main(void)
{
	long double experimental = 0.7468241328124271;
	int row = 8;

	printf("Romberg integration of f(x) = exp(- x * x)\n");
	printf("Internet value = 0.7468241328124271\n");
	printf("from website = https://www.integral-calculator.com/\n");
	printf("Approximation\t\tPercent Difference\tSteps\tRuntime (s)\n");
	for (long steps = 100000; steps <= 900000; steps += 100000)
	{
		clock_t start = clock();
		long double trial = romberg(0.0, 1.0, steps, row);
		clock_t endTm = clock();
		long double runtime = ((long double)endTm - start) /
			CLOCKS_PER_SEC;
		long double pd = 100.0 * fabsl(experimental - trial) /
			(0.5 *(experimental + trial));
		printf("%16.15lf\t%16.15le\t%6d\t%4.3lf\n", 
			trial, pd, steps, runtime);
	}
	return(0);
}

Function Optimization by James Pate Williams, Jr.

We use two numerical analysis algorithms and two artificial intelligent methods. The first two techniques are from the excellent tome “A Numerical Library in C for Scientists and Engineers” by H.T. Lau, PhD. One does not utilize any derivatives and the other requires first partial derivatives. My homegrown algorithms are my implementation of the particle swarm optimization and an evolutionary hill-climber. The output from my C# application is below:

The two numerical algorithms PRAXIS and FLEMIN perform fairly well and do not require that many function evaluations which is better metric than wall-clock time.

See the excellent webpage: https://en.wikipedia.org/wiki/Test_functions_for_optimization

Gibb’s Overshoot Phenomenon by James Pate Williams, Jr.

The Gibb’s overshoot phenomenon is very well defined by the Fourier series expansion of the Heaviside step function. The following images were produced by a new C/C++ Win32 application that I created in the period February 10, 2023, to February 14, 2023.

The code was adapted from “Fourier Series and Boundary Value Problems” by Ruel V. Churchill and James Ward Brown. The integration algorithm was translated from C source code found in the tome “A Numerical Library in C for Scientists and Engineers” by H.T. Lau, PhD pages 299-303.