Blog Entry © Tuesday, August 19, 2025, by James Pate Williams, Jr., Continuation of Answers to the Exercises in Chapter 1 of Modern Quantum Chemistry an Introduction to Advanced Electronic Structure Theory by Attila Szabo and Neil S. Ostlund

Note: Later on, Tuesday, August 19, 2025, I added five C++ source code files.

Blog Entry Tuesday August 19 2025Download
#include <vector>
#include <random>

class DblLinearAlgebra
{
public:
    static void DblPrintMatrix(
        int m, int n, std::vector<std::vector<double>>& A);
    static void DblAddition(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblAnticommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblCommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static double DblDeterminant(
        int n, int row, int col,
        std::vector<std::vector<double>>& A);
    static bool DblGaussianElimination(
        int m, int n, std::vector<std::vector<double>>& A,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblGaussianFactor(
        int n, std::vector<std::vector<double>>& M,
        std::vector<size_t>& pivot);
    static bool DblGaussianSolution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblSubstitution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblInverse(
        int n, std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& A);
    static void DblCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<double>>& C,
        std::vector<std::vector<double>>& I,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& adjoint,
        std::vector<double>& a);
    static void DblMatrixKernel(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& X,
        size_t& r);
    static void DblMatrixImage(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& N,
        std::vector<std::vector<double>>& X,
        int rank);
    static void DblGenerateNonSingular(
        double scale, double& determinant,
        int n, unsigned int seed,
        std::vector<std::vector<double>>& Mr);
};

#include "DblLinearAlgebra.h"
#include <iomanip>
#include <iostream>

void DblLinearAlgebra::DblPrintMatrix(
    int m, int n, std::vector<std::vector<double>>& A)
{
    for (int i = 0; i < m; i++)
    {
        for (int j = 0; j < n; j++)
        {
            std::cout << std::setprecision(6) << std::setw(9);

            if (fabs(A[i][j]) > 1.0e-12)
            {
                std::cout << A[i][j] << ' ';
            }

            else
            {
                std::cout << 0 << ' ';
            }
        }

        std::cout << std::endl;
    }
}

void DblLinearAlgebra::DblAddition(
    size_t m, size_t n,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            C[i][j] = A[i][j] + B[i][j];
        }
    }
}

void DblLinearAlgebra::DblSubtraction(
    size_t m, size_t n,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            C[i][j] = A[i][j] - B[i][j];
        }
    }
}

void DblLinearAlgebra::DblMultiply(
    size_t m, size_t n, size_t p,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            double sum = { 0 };

            for (size_t k = 0; k < p; k++)
            {
                sum += A[i][k] * B[k][j];
            }

            C[i][j] = sum;
        }
    }
}

void DblLinearAlgebra::DblAnticommutator(
    size_t n,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    std::vector<std::vector<double>> D(n,
        std::vector<double>(n));
    std::vector<std::vector<double>> E(n,
        std::vector<double>(n));

    DblMultiply(n, n, n, A, B, D);
    DblMultiply(n, n, n, B, A, E);
    DblAddition(n, n, D, E, C);
}

void DblLinearAlgebra::DblCommutator(
    size_t n,
    std::vector<std::vector<double>>& A,
    std::vector<std::vector<double>>& B,
    std::vector<std::vector<double>>& C)
{
    std::vector<std::vector<double>> D(n,
        std::vector<double>(n));
    std::vector<std::vector<double>> E(n,
        std::vector<double>(n));

    DblMultiply(n, n, n, A, B, D);
    DblMultiply(n, n, n, B, A, E);
    DblSubtraction(n, n, D, E, C);
}

// https://www.geeksforgeeks.org/dsa/determinant-of-a-matrix/

double getDet(std::vector<std::vector<double>>& mat, int n) {

    // Base case: if the matrix is 1x1
    if (n == 1) {
        return mat[0][0];
    }

    // Base case for 2x2 matrix
    if (n == 2) {
        return mat[0][0] * mat[1][1] -
            mat[0][1] * mat[1][0];
    }

    // Recursive case for larger matrices
    double  res = 0;
    for (int col = 0; col < n; ++col) {

        // Create a submatrix by removing the first 
        // row and the current column
        std::vector<std::vector<double>> sub(n - 1,
            std::vector<double>(n - 1));
        for (int i = 1; i < n; ++i) {
            int subcol = 0;
            for (int j = 0; j < n; ++j) {

                // Skip the current column
                if (j == col) continue;

                // Fill the submatrix
                sub[i - 1LL][subcol++] = mat[i][j];
            }
        }

        // Cofactor expansion
        int sign = (col % 2 == 0) ? 1 : -1;
        res += sign * mat[0][col] * getDet(sub, n - 1);
    }

    return res;
}

double DblLinearAlgebra::DblDeterminant(
    int n, int row, int col,
    std::vector<std::vector<double>>& A)
{
    return getDet(A, A.size());
}

bool DblLinearAlgebra::DblGaussianElimination(
    int m, int n, std::vector<std::vector<double>>& M,
    std::vector<double>& b, std::vector<double>& x, 
    std::vector<size_t>& pivot)
{
    bool failure = false;
    std::vector<double> c(m);

    b.resize(n);
    x.resize(n);

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

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

        while (i < n && !found)
        {
            if (M[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[i][l];
                M[i][l] = M[j][l];
                M[j][l] = t;
                t = b[i];
                b[i] = b[j];
                b[j] = t;
            }
        }

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

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

        for (size_t k = j + 1; k < n; k++)
        {
            for (size_t l = j + 1; l < n; l++)
                M[k][l] = M[k][l] - c[k] * M[j][l];

            b[k] = b[k] - c[k] * b[j];
        }
    }

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

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

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

    return failure;
}

bool DblLinearAlgebra::DblGaussianFactor(
    int n, std::vector<std::vector<double>>& M,
    std::vector<size_t>& pivot)
{
    // returns false if matrix is singular
    std::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[i] = i;
        row_max = 0;
        
        for (size_t j = 0; j < n; j++)
            row_max = fmax(row_max, fabs(M[i][j]));
        
        if (row_max == 0)
        {
            flag = 0;
            row_max = 1.0;
        }

        d[i] = row_max;
    }

    if (n <= 1) return flag != 0;
    
    // factorization
    
    for (size_t k = 0; k < (size_t)(n - 1LL); k++)
    {
        // determine pivot row the row i_star

        col_max = fabs(M[k][k]) / d[k];
        i_star = k;

        for (size_t i = k + 1; i < n; i++)
        {
            awikod = fabs(M[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;
                itemp = pivot[i_star];
                pivot[i_star] = pivot[k];
                pivot[k] = itemp;
                temp = d[i_star];
                d[i_star] = d[k];
                d[k] = temp;

                for (size_t j = 0; j < n; j++)
                {
                    temp = M[i_star][j];
                    M[i_star][j] = M[k][j];
                    M[k][j] = temp;
                }
            }

            // eliminate x[k]

            for (size_t i = k + 1; i < n; i++)
            {
                M[i][k] /= M[k][k];
                ratio = M[i][k];

                for (size_t j = k + 1; j < n; j++)
                    M[i][j] -= ratio * M[k][j];
            }
        }

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

    if (flag == 0)
        return false;

    return true;
}

bool DblLinearAlgebra::DblGaussianSolution(
    int n, std::vector<std::vector<double>>& M,
    std::vector<double>& b, std::vector<double>& x,
    std::vector<size_t>& pivot)
{
    if (!DblGaussianFactor(n, M, pivot))
        return false;

    return DblSubstitution(n, M, b, x, pivot);
}

bool DblLinearAlgebra::DblSubstitution(
    int n, std::vector<std::vector<double>>& M,
    std::vector<double>& b, std::vector<double>& x,
    std::vector<size_t>& pivot)
{
    double sum = 0.0;
    size_t n1 = n - 1LL;

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

    // forward substitution

    x[0] = b[pivot[0]];

    for (size_t i = 1; i < n; i++)
    {
        double sum = 0.0;

        for (size_t j = 0; j < i; j++)
            sum += M[i][j] * x[j];

        x[i] = b[pivot[i]] - sum;
    }

    // backward substitution

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

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

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

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

    return true;
}

bool DblLinearAlgebra::DblInverse(
    int n, std::vector<std::vector<double>>& M,
    std::vector<std::vector<double>>& Mi)
{
    std::vector<double> b(n);
    std::vector<double> x(n);
    std::vector<size_t> pivot(n);
    std::vector<std::vector<double>> Mc(n,
        std::vector<double>(n));

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            Mc[i][j] = M[i][j];
        }
    }

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

    for (size_t i = 0; i < n; i++)
    {
        b[i] = 0;
    }

    for (size_t i = 0; i < n; i++)
    {
        b[i] = 1;

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

        b[i] = 0;

        for (size_t j = 0; j < n; j++)
            Mi[j][i] = x[j];
    }

    return true;
}

void DblLinearAlgebra::DblCharPolyAndAdjoint(
    int n,
    std::vector<std::vector<double>>& C,
    std::vector<std::vector<double>>& I,
    std::vector<std::vector<double>>& M,
    std::vector<std::vector<double>>& adjoint,
    std::vector<double>& a)
{
    C.resize(n, std::vector<double>(n));
    I.resize(n, std::vector<double>(n));

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            C[i][j] = I[i][j] = 0;
    }

    for (size_t i = 0; i < n; i++)
        C[i][i] = I[i][i] = 1;

    a[0] = 1;

    for (size_t i = 1; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
            {
                double sum = 0.0;

                for (size_t l = 0; l < n; l++)
                    sum += M[j][l] * C[l][k];

                C[j][k] = sum;
            }
        }

        double tr = 0.0;

        for (size_t j = 0; j < n; j++)
            tr += C[j][j];

        a[i] = -tr / i;

        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
                C[j][k] += a[i] * I[j][k];
        }
    }

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

            for (size_t k = 0; k < n; k++)
                sum += M[i][k] * C[k][j];

            C[i][j] = sum;
        }
    }

    double trace = 0.0;

    for (size_t i = 0; i < n; i++)
        trace += C[i][i];

    trace /= n;
    a[n - 1LL] = -trace;

    double factor = 1.0;

    if ((n - 1) % 2 != 0)
        factor = -1.0;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            adjoint[i][j] = factor * C[i][j];
    }
}

void DblLinearAlgebra::DblMatrixKernel(
    int m, int n,
    std::vector<std::vector<double>>& M,
    std::vector<std::vector<double>>& X,
    size_t& r)
{
    double D = 0.0;
    std::vector <int> c(m);
    std::vector <int> d(n);

    r = 0;

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    size_t j, k = 1;

Step2:
    
    for (j = 0; j < m; j++)
    {
        if (M[j][k] != 0 && c[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d[k] = 0;
        goto Step4;
    }

    D = -1.0 / M[j][k];

    M[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        M[j][s] = D * M[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = M[i][k];
                M[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            M[i][s] += D * M[j][s];
        }
    }

    c[j] = (int)k;
    d[k] = (int)j;

Step4:
    
    if (k < n - 1)
    {
        k++;
        goto Step2;
    }

    X.resize(n, std::vector<double>(n));

    if (r != 0)
    {
        for (k = 0; k < n; k++)
        {
            if (d[k] == 0)
            {
                for (size_t i = 0; i < n; i++)
                {
                    if (d[i] > 0)
                        X[k][i] = M[d[i]][k];
                    else if (i == k)
                        X[k][i] = 1;
                    else
                        X[k][i] = 0;
                }
            }
        }
    }
}

void DblLinearAlgebra::DblMatrixImage(
    int m, int n,
    std::vector<std::vector<double>>& M,
    std::vector<std::vector<double>>& N,
    std::vector<std::vector<double>>& X,
    int rank)
{
    double D = 0.0;
    size_t r = 0;
    std::vector<std::vector<double>> copyM(
        m, std::vector<double>(n));
    std::vector <int> c(m);
    std::vector <int> d(n);

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    size_t j = 0, k = 1;

    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            N[i][j] = copyM[i][j] = M[i][j];
        }
    }

Step2:
    
    for (size_t j = 0; j < m; j++)
    {
        if (copyM[j][k] != 0 && c[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d[k] = 0;
        goto Step4;
    }

    D = -1.0 / copyM[j][k];

    copyM[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        copyM[j][s] = D * copyM[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = copyM[i][k];
                copyM[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            copyM[i][s] += D * copyM[j][s];
        }
    }

    c[j] = (int)k;
    d[k] = (int)j;

Step4:
    
    if (k < (size_t)(n - 1LL))
    {
        k++;
        goto Step2;
    }

    rank = (int)(n - r) ;

    for (j = 0; j < m; j++)
    {
        if (c[j] != 0)
        {
            for (size_t i = 0; i < m; i++)
            {
                N[i][c[j]] = M[i][c[j]];
            }
        }
    }
}

void DblLinearAlgebra::DblGenerateNonSingular(
    double scale, double& determinant,
    int n, unsigned int seed,
    std::vector<std::vector<double>>& Mr)
{
    bool failure = false;
    std::mt19937 rng(seed);
    std::uniform_real_distribution<double> dist(0.0, 1.0);

    while (true)
    {
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                Mr[i][j] = scale * dist(rng);
            }
        }

        determinant = DblDeterminant(n, 0, 0, Mr);
        failure = determinant == 0;

        if (!failure)
            return;
    }
}

#include <complex>
#include <vector>

class CmpLinearAlgebra
{
public:
    static void CmpPrintMatrix(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& Ac);
    static void CmpAddition(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpAnticommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpCommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static std::complex<double> CmpDeterminant(
        int n,
        std::vector<std::vector<std::complex<double>>>& Ac);
    static void CmpAdjoint(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& Ac,
        std::vector<std::vector<std::complex<double>>>& Ad);
    static bool CmpGaussianElimination(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpGaussianFactor(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<size_t>& pivot);
    static bool CmpGaussianSolution(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpSubstitution(
        int m, int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpInverse(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& Mi);
    static void CmpCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<std::complex<double>>>& C,
        std::vector<std::vector<std::complex<double>>>& I,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& adjoint,
        std::vector<std::complex<double>>& a);
    static void CmpMatrixKernel(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& X,
        size_t& r);
    static void CmpMatrixImage(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& N,
        std::vector<std::vector<std::complex<double>>>& X,
        int rank);
    static void CmpGenerateNonSingular(
        double scale, std::complex<double>& determinant,
        int n, unsigned int seed,
        std::vector<std::vector<std::complex<double>>>& Mc);
};

#include "CmpLinearAlgebra.h"
#include <iomanip>
#include <iostream>
#include <random>

void CmpLinearAlgebra::CmpPrintMatrix(
    int m, int n,
    std::vector<std::vector<std::complex<double>>>& Ac)
{
    for (int i = 0; i < m; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (Ac[i][j]._Val[0] >= 0)
                std::cout << '+';
            else
                std::cout << '-';

            std::cout << std::setprecision(6) << std::setw(9);

            if (fabs(Ac[i][j]._Val[0]) > 1.0e-12)
            {
                std::cout << fabs(Ac[i][j]._Val[0]) << ' ';
            }

            else
            {
                std::cout << 0 << ' ';
            }

            if (Ac[i][j]._Val[1] >= 0)
                std::cout << '+';
            else
                std::cout << '-';

            std::cout << std::setprecision(6) << std::setw(9);

            if (fabs(Ac[i][j]._Val[1]) > 1.0e-12)
            {
                std::cout << fabs(Ac[i][j]._Val[1]) << "i\t";
            }

            else
            {
                std::cout << 0 << "i\t";
            }
        }

        std::cout << std::endl;
    }
}

void CmpLinearAlgebra::CmpAddition(
    size_t m, size_t n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            C[i][j] = A[i][j] + B[i][j];
        }
    }
}

void CmpLinearAlgebra::CmpSubtraction(
    size_t m, size_t n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            C[i][j] = A[i][j] - B[i][j];
        }
    }
}

void CmpLinearAlgebra::CmpMultiply(
    size_t m, size_t n, size_t p,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            std::complex<double> sum = 0;

            for (size_t k = 0; k < p; k++)
            {
                sum += A[i][k] * B[k][j];
            }

            C[i][j] = sum;
        }
    }
}

void CmpLinearAlgebra::CmpAnticommutator(
    size_t n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    std::vector<std::vector<std::complex<double>>> D(n,
        std::vector<std::complex<double>>(n));
    std::vector<std::vector<std::complex<double>>> E(n,
        std::vector<std::complex<double>>(n));

    CmpMultiply(n, n, n, A, B, D);
    CmpMultiply(n, n, n, B, A, E);
    CmpAddition(n, n, D, E, C);
}

void CmpLinearAlgebra::CmpCommutator(
    size_t n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::vector<std::complex<double>>>& B,
    std::vector<std::vector<std::complex<double>>>& C)
{
    std::vector<std::vector<std::complex<double>>> D(n,
        std::vector<std::complex<double>>(n));
    std::vector<std::vector<std::complex<double>>> E(n,
        std::vector<std::complex<double>>(n));

    CmpMultiply(n, n, n, A, B, D);
    CmpMultiply(n, n, n, B, A, E);
    CmpSubtraction(n, n, D, E, C);
}

void CmpLinearAlgebra::CmpAdjoint(
    size_t m, size_t n,
    std::vector<std::vector<std::complex<double>>>& Ac,
    std::vector<std::vector<std::complex<double>>>& Ad)
{
    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            Ad[j][i] = std::conj(Ac[i][j]);
        }
    }
}

// https://www.geeksforgeeks.org/dsa/determinant-of-a-matrix/

std::complex<double> getDet(
    std::vector<std::vector<std::complex<double>>>& mat, int n) {

    // Base case: if the matrix is 1x1
    if (n == 1) {
        return mat[0][0];
    }

    // Base case for 2x2 matrix
    if (n == 2) {
        return mat[0][0] * mat[1][1] -
            mat[0][1] * mat[1][0];
    }

    // Recursive case for larger matrices
    std::complex<double> res = 0;
    for (int col = 0; col < n; ++col) {

        // Create a submatrix by removing the first 
        // row and the current column
        std::vector<std::vector<std::complex<double>>> sub(n - 1,
            std::vector<std::complex<double>>(n - 1));
        for (int i = 1; i < n; ++i) {
            int subcol = 0;
            for (int j = 0; j < n; ++j) {

                // Skip the current column
                if (j == col) continue;

                // Fill the submatrix
                sub[i - 1LL][subcol++] = mat[i][j];
            }
        }

        // Cofactor expansion
        int sign = (col % 2 == 0) ? 1 : -1;
        std::complex<double> csign(sign, 0.0);
        res = res + csign * mat[0][col] * getDet(sub, n - 1);
    }

    return res;
}

std::complex<double> CmpLinearAlgebra::CmpDeterminant(
    int n, std::vector<std::vector<std::complex<double>>>& A)
{
    return getDet(A, A.size());
}

bool CmpLinearAlgebra::CmpGaussianElimination(
    int m, int n,
    std::vector<std::vector<std::complex<double>>>& A,
    std::vector<std::complex<double>>& b,
    std::vector<std::complex<double>>& x,
    std::vector<size_t>& pivot)
{
    bool failure = false;
    std::vector<std::complex<double>> c(m);

    b.resize(n);
    x.resize(n);

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

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

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

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

        if (i > j)
        {
            for (size_t l = j; l < n; l++)
            {
                std::complex<double> t = A[i][l];
                A[i][l] = A[j][l];
                A[j][l] = t;
                t = b[i];
                b[i] = b[j];
                b[j] = t;
            }
        }

        std::complex<double> d = 1.0 / A[j][j];

        for (size_t k = j + 1; k < n; k++)
            c[k] = d * A[k][j];

        for (size_t k = j + 1; k < n; k++)
        {
            for (size_t l = j + 1; l < n; l++)
                A[k][l] = A[k][l] - c[k] * A[j][l];

            b[k] = b[k] - c[k] * b[j];
        }
    }

    for (long long i = (long long)n - 1; i >= 0; i--)
    {
        std::complex<double> sum = 0;

        for (size_t j = i + 1; j < n; j++)
            sum += A[i][j] * x[j];

        x[i] = (b[i] - sum) / A[i][i];
    }

    return failure;
}

bool CmpLinearAlgebra::CmpSubstitution(
    int m, int n, std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::complex<double>>& b,
    std::vector<std::complex<double>>& x,
    std::vector<size_t>& pivot)
{
    std::complex<double> sum = 0;
    size_t n1 = n - 1LL;

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

    // forward substitution

    x[0] = b[pivot[0]];

    for (size_t i = 1; i < n; i++)
    {
        std::complex<double> sum = 0;

        for (size_t j = 0; j < i; j++)
            sum += M[i][j] * x[j];

        x[i] = b[pivot[i]] - sum;
    }

    // backward substitution

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

    for (long long i = n - 2LL; i >= 0; i--)
    {
        std::complex<double> sum = 0;

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

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

    return true;
}

static std::complex<double> complex_max(
    std::complex<double> a, std::complex<double> b) {
    return (std::abs(a) > std::abs(b)) ? a : b;
}

bool CmpLinearAlgebra::CmpGaussianFactor(
    int n, std::vector<std::vector<std::complex<double>>>& M,
    std::vector<size_t>& pivot)
{
    // returns false if matrix is singular
    std::vector<std::complex<double>> d(n);
    std::complex<double> awikod = 0, col_max = 0, ratio = 0, row_max = 0, temp = 0;
    int flag = 1;
    size_t i_star, itemp;

    for (size_t i = 0; i < n; i++)
    {
        pivot[i] = i;
        row_max = 0;

        for (size_t j = 0; j < n; j++)
            row_max = complex_max(row_max, abs(M[i][j]));
        
        if (abs(row_max) == 0)
        {
            flag = 0;
            row_max = 1;
        }

        d[i] = row_max;
    }
    if (n <= 1) return flag != 0;
    
    // factorization
    
    for (size_t k = 0; k < (size_t)n - 1LL; k++)
    {
        // determine pivot row the row i_star

        col_max = abs(M[k][k]) / d[k];
        i_star = k;

        for (size_t i = k + 1; i < n; i++)
        {
            awikod = abs(M[i][k]) / d[i];

            if (abs(awikod) > abs(col_max))
            {
                col_max = awikod;
                i_star = i;
            }
        }
        
        if (abs(col_max) == 0)
            flag = 0;
        
        else
        {
            if (i_star > k)
            {
                // make k the pivot row by
                // interchanging with i_star
                flag *= -1;
                itemp = pivot[i_star];
                pivot[i_star] = pivot[k];
                pivot[k] = itemp;
                temp = d[i_star];
                d[i_star] = d[k];
                d[k] = temp;

                for (size_t j = 0; j < n; j++)
                {
                    temp = M[i_star][j];
                    M[i_star][j] = M[k][j];
                    M[k][j] = temp;
                }
            }

            // eliminate x[k]
            
            for (size_t i = k + 1; i < n; i++)
            {
                M[i][k] /= M[k][k];
                ratio = M[i][k];
                
                for (size_t j = k + 1; j < n; j++)
                    M[i][j] -= ratio * M[k][j];
            }
        }

        if (abs(M[n - 1LL][n - 1LL]) == 0) flag = 0;
    }

    if (flag == 0)
        return false;

    return true;
}

bool CmpLinearAlgebra::CmpGaussianSolution(
    int n, std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::complex<double>>& b,
    std::vector<std::complex<double>>& x,
    std::vector<size_t>& pivot)
{
    if (!CmpGaussianFactor(n, M, pivot))
        return false;

    return CmpSubstitution(n, n, M, b, x, pivot);
}

bool CmpLinearAlgebra::CmpInverse(
    int n, std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::vector<std::complex<double>>>& Mi)
{
    std::vector<std::complex<double>> b(n);
    std::vector<std::complex<double>> x(n);
    std::vector<size_t> pivot(n);
    std::vector<std::vector<std::complex<double>>> Mc(n,
        std::vector<std::complex<double>>(n));

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            Mc[i][j] = M[i][j];
        }
    }

    if (!CmpGaussianFactor(n, Mc, pivot))
        return false;

    for (size_t i = 0; i < n; i++)
    {
        b[i] = 0;
    }

    for (size_t i = 0; i < n; i++)
    {
        b[i] = 1;

        if (!CmpSubstitution(n, n, Mc, b, x, pivot))
            return false;

        b[i] = 0;

        for (size_t j = 0; j < n; j++)
            Mi[j][i] = x[j];
    }

    return true;
}

void CmpLinearAlgebra::CmpCharPolyAndAdjoint(
    int n,
    std::vector<std::vector<std::complex<double>>>& C,
    std::vector<std::vector<std::complex<double>>>& I,
    std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::vector<std::complex<double>>>& adjoint,
    std::vector<std::complex<double>>& a)
{
    C.resize(n, std::vector<std::complex<double>>(n));
    I.resize(n, std::vector<std::complex<double>>(n));

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            C[i][j] = I[i][j] = 0;
    }

    for (size_t i = 0; i < n; i++)
        C[i][i] = I[i][i] = 1;

    a[0] = 1;

    for (size_t i = 1; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
            {
                std::complex<double> sum = 0.0;

                for (size_t l = 0; l < n; l++)
                    sum += M[j][l] * C[l][k];

                C[j][k] = sum;
            }
        }

        std::complex<double> tr = 0.0;

        for (size_t j = 0; j < n; j++)
            tr += C[j][j];

        std::complex<double> ci = 0;
        ci._Val[0] = (double)i;

        a[i] = -tr / ci;

        for (size_t j = 0; j < n; j++)
        {
            for (size_t k = 0; k < n; k++)
                C[j][k] += a[i] * I[j][k];
        }
    }

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            std::complex<double> sum = 0.0;

            for (size_t k = 0; k < n; k++)
                sum += M[i][k] * C[k][j];

            C[i][j] = sum;
        }
    }

    std::complex<double> trace = 0.0;

    for (size_t i = 0; i < n; i++)
        trace += C[i][i];

    trace /= n;
    a[n - 1LL] = -trace;

    std::complex<double> factor = 1.0;

    if ((n - 1) % 2 != 0)
        factor = -1.0;

    for (size_t i = 0; i < n; i++)
    {
        for (size_t j = 0; j < n; j++)
            adjoint[i][j] = factor * C[i][j];
    }
}

void CmpLinearAlgebra::CmpMatrixKernel(
    int m, int n,
    std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::vector<std::complex<double>>>& X,
    size_t& r)
{
    std::complex<double> D = 0;
    std::vector<int> c(m);
    std::vector<int> d(n);

    r = 0;

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    size_t j = 0, k = 1;

Step2:

    for (j = 0; j < m; j++)
    {
        if (abs(M[j][k]) != 0 && c[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d[k] = 0;
        goto Step4;
    }

    D = -1.0 / M[j][k];

    M[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        M[j][s] = D * M[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = M[i][k];
                M[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            M[i][s] += D * M[j][s];
        }
    }

    c[j] = (int)k;
    d[k] = (int)j;

Step4:

    if (k < n - 1)
    {
        k++;
        goto Step2;
    }

    X.resize(n, std::vector<std::complex<double>>(n));

    if (r != 0)
    {
        for (k = 0; k < n; k++)
        {
            if (d[k] == 0)
            {
                for (size_t i = 0; i < n; i++)
                {
                    if (d[i] > 0)
                        X[k][i] = M[d[i]][k];
                    else if (i == k)
                        X[k][i] = 1;
                    else
                        X[k][i] = 0;
                }
            }
        }
    }
}

void CmpLinearAlgebra::CmpMatrixImage(
    int m, int n,
    std::vector<std::vector<std::complex<double>>>& M,
    std::vector<std::vector<std::complex<double>>>& N,
    std::vector<std::vector<std::complex<double>>>& X,
    int rank)
{
    std::complex<double> D = 0.0;
    size_t r = 0;
    std::vector<std::vector<std::complex<double>>> copyM(
        m, std::vector<std::complex<double>>(n));
    std::vector<int> c(m);
    std::vector<int> d(n);

    for (size_t i = 0; i < m; i++)
        c[i] = -1;

    size_t j = 0, k = 1;

    for (size_t i = 0; i < m; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            N[i][j] = copyM[i][j] = M[i][j];
        }
    }

Step2:

    for (size_t j = 0; j < m; j++)
    {
        if (abs(copyM[j][k]) != 0 && c[j] == -1)
            break;
    }

    if (j == m)
    {
        r++;
        d[k] = 0;
        goto Step4;
    }

    D = -1.0 / copyM[j][k];

    copyM[j][k] = -1;

    for (size_t s = k + 1; s < n; s++)
    {
        copyM[j][s] = D * copyM[j][s];

        for (size_t i = 0; i < m; i++)
        {
            if (i != j)
            {
                D = copyM[i][k];
                copyM[i][k] = 0;
            }
        }
    }

    for (size_t s = k + 1; s < n; s++)
    {
        for (size_t i = 0; i < m; i++)
        {
            copyM[i][s] += D * copyM[j][s];
        }
    }

    c[j] = (int)k;
    d[k] = (int)j;

Step4:
    
    if (k < (size_t)(n - 1LL))
    {
        k++;
        goto Step2;
    }

    rank = (int)(n - r);

    for (j = 0; j < m; j++)
    {
        if (c[j] != 0)
        {
            for (size_t i = 0; i < m; i++)
            {
                N[i][c[j]] = M[i][c[j]];
            }
        }
    }
}

void CmpLinearAlgebra::CmpGenerateNonSingular(
    double scale, std::complex<double>& cDeterminant,
    int n, unsigned int seed,
    std::vector<std::vector<std::complex<double>>>& Mc)
{
    bool failure = false;
    std::mt19937 rng(seed);
    std::uniform_real_distribution<double> dist(0.0, 1.0);

    while (true)
    {
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                Mc[i][j]._Val[0] = scale * dist(rng);
                Mc[i][j]._Val[1] = scale * dist(rng);
            }
        }

       cDeterminant = CmpDeterminant(n, Mc);

       if (cDeterminant._Val[0] != 0 || cDeterminant._Val[1] != 0)
           break;
    }
}

// Exercises from "Modern Quantum Chemistry an Introduction to Advanced
// Electronic Structure Theory" by Attila Szabo and Neil S. Ostlund
// https://chemistlibrary.wordpress.com/wp-content/uploads/2015/02/modern-quantum-chemistry.pdf
// Program (c) Tuesday, August 12, 2025 by James Pate Williams, Jr.
// Program (c) Tuesday, August 19, 2025 by James Pate Williams, Jr.

#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"

int main()
{
	// static data matrices
	double AArcb[3][3] = { { 2, 3, -1 }, { 4, 4, -3 }, { -2, 3, -1 } };
	double AArso[3][3] = { { 1, 1, 0 }, { 1, 2, 2 }, { 0, 2, -1 } };
	double BBrso[3][3] = { { 1, -1, 1 }, { -1, 0, 0 }, { 1, 0, 1} };
	double BBr[3][3] = { { 1, -1, 1 }, { -1 , 0, 0 }, { 1, 0, 1 } };
	double AAcr[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	double AAci[3][3] = { { 1, 1, 2 }, { 3, 0, 1 }, { 0, 2, 4 } };
	double BBcr[3][3] = { { 1, 0, 1 }, { 1 , 1, 0 }, { 0, 1, 1 } };
	double BBci[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	// some array dimensions
	int m = 3, n = 3, p = 3;
	// a couple of 3x1 vectors
	std::vector<double> br(3);
	std::vector<size_t> pivot(3);
	// 3x3 real matrices
	std::vector<std::vector<double>> Arcb(3, std::vector<double>(3));
	std::vector<std::vector<double>> Arso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Brso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Br(3, std::vector<double>(3));
	std::vector<std::vector<double>> Cr(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ai(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ari(3, std::vector<double>(3));
	// a 4x4 real matrix
	std::vector<std::vector<double>> Mr(4, std::vector<double>(4));
	// 3x3 complex matrices
	std::vector<std::vector<std::complex<double>>> Ac(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Bc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Cc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Dc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Ec(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Fc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Gc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Mc(4,
		std::vector<std::complex<double>>(4));
	// copy static real matrices to dynamic matrices
	for (int i = 0; i < m; i++)
	{
		for (int j = 0; j < p; j++)
		{
			Arcb[i][j] = AArcb[i][j];
			Arso[i][j] = AArso[i][j];
			Brso[i][j] = BBrso[i][j];
			Ac[i][j]._Val[0] = AAcr[i][j];
			Ac[i][j]._Val[1] = AAci[i][j];
		}
	}
	// copy static complex matrices to dynamic matrices
	for (int i = 0; i < p; i++)
	{
		for (int j = 0; j < n; j++)
		{
			Br[i][j] = BBr[i][j];
			Bc[i][j]._Val[0] = BBcr[i][j];
			Bc[i][j]._Val[1] = BBci[i][j];
		}
	}
	// See "Elementary Numerical Analysis an
	// Algorithmic Approach" (c) 1980 by S. D. Conte
	// and Carl de Boor
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Br, Cr);
	std::cout << "Ar * Br = Cr Conte & de Boor" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	// complex matrix multiplication
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << "Ac * Bc = Cc" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
	std::cout << std::endl;
	// Exercise 1.2 from Szabo and Ostlund
	std::cout << "Exercise 1.2 page 5 Commutator" << std::endl;
	DblLinearAlgebra::DblCommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	std::cout << "Exercise 1.2 page 5 Anticommutator" << std::endl;
	DblLinearAlgebra::DblAnticommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Cc, Dc);
	std::cout << "Exercise 1.3 page 6 Cc adjoint" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Dc);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Ac, Ec);
	CmpLinearAlgebra::CmpAdjoint(3, 3, Bc, Fc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Fc, Ec, Gc);
	std::cout << "Exercise 1.3 page 6 Bc adjoint * Ac adjoint"
		<< std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Gc);
	std::cout << std::endl;
	std::cout << "Ar matrix" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Arcb);
	bool inv = DblLinearAlgebra::DblInverse(n, Arcb, Ai);
	std::cout << std::endl;
	std::cout << "Ar Conte & de Boor inverse flag = "
		<< inv << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ai);
	std::cout << std::endl;
	std::cout << "Ar * Ar inverse" << std::endl;
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Ai, Ari);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ari);
	std::cout << std::endl;
	double rDeterminant = 0;
	DblLinearAlgebra::DblGenerateNonSingular(
		2.0, rDeterminant, 4, 1, Mr);
	std::cout << "rDeterminant = ";
	std::cout << rDeterminant << std::endl;
	std::cout << std::endl;
	std::cout << "Ac" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
	std::cout << std::endl;
	inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
	std::cout << "Ac inverse flag = " << inv << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Bc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << std::endl;
	std::cout << "Ac * Ac inverse" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
	std::complex<double> cDeterminant = 0;
	CmpLinearAlgebra::CmpGenerateNonSingular(
		2.0, cDeterminant, 4, 1, Mc);
	std::cout << std::endl;
	std::cout << "complex determinant = ";
	std::cout << cDeterminant << std::endl;
	double rDeterminantA = 0;
	std::vector<std::vector<double>> A44r(4,
		std::vector<double>(4));
	DblLinearAlgebra::DblGenerateNonSingular(
		2.0, rDeterminantA, 4, 2, A44r);
	double rDeterminantB = 0;
	std::vector<std::vector<double>> B44r(4,
		std::vector<double>(4));
	DblLinearAlgebra::DblGenerateNonSingular(
		2.0, rDeterminantB, 4, 3, B44r);
	std::cout << std::endl;
	std::vector<std::vector<double>> C44r(4,
		std::vector<double>(4));
	DblLinearAlgebra::DblMultiply(4, 4, 4, A44r, B44r, C44r);
	std::cout << "|A| = " << rDeterminantA << std::endl;
	std::cout << "|B| = " << rDeterminantB << std::endl;
	bool failure = false;
	double rDeterminantC =
		DblLinearAlgebra::DblDeterminant(4, 0, 0, C44r);
	std::cout << "|AB| = " << rDeterminantC << std::endl;
	std::cout << "|A||B| = " << rDeterminantA *
		rDeterminantB << std::endl;
	// Exercise 1.6 with 4x4 complex determinants
	std::vector<std::vector<std::complex<double>>> A44c(4,
		std::vector<std::complex<double>>(4));
	std::cout << std::endl;
	std::complex<double> cDeterminantA = 0;
	CmpLinearAlgebra::CmpGenerateNonSingular(
		2.0, cDeterminantA, 4, 2, A44c);
	std::vector<std::vector<std::complex<double>>> B44c(4,
		std::vector<std::complex<double>>(4));
	std::complex<double> cDeterminantB = 0;
	CmpLinearAlgebra::CmpGenerateNonSingular(
		2.0, cDeterminantB, 4, 3, B44c);
	std::vector<std::vector<std::complex<double>>> C44c(4,
		std::vector<std::complex<double>>(4));
	CmpLinearAlgebra::CmpMultiply(4, 4, 4, A44c, B44c, C44c);
	std::cout << "|A| = " << cDeterminantA << std::endl;
	std::cout << "|B| = " << cDeterminantB << std::endl;
	failure = false;
	std::complex<double> cDeterminantC =
		CmpLinearAlgebra::CmpDeterminant(4, C44c);
	std::cout << "|AB| = " << cDeterminantC << std::endl;
	std::cout << "|A||B| = " << cDeterminantA *
		cDeterminantB << std::endl;
	std::cout << "\nEnter any key to halt: ";
	char line[128] = "";
	std::cin.getline(line, 128);
}
Unknown's avatar

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.

Leave a comment