Category: Uncategorized

Blog Entry Tuesday, June 16, 2026, by James Pate Williams, Jr., Linear Systems of Equations Solutions

Blog Entry Wednesday June 17 2026Download
#pragma once
#include <vector>

class DirectMethods
{
public:
    static void Substitute(
        int n,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x,
        std::vector<int>& pivot);
    static bool GaussianElimimation(
        int n,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x);
    static void LUDecomposition(
        int n,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x);
};


#include "DirectMethods.h"
#include <vector>

void DirectMethods::Substitute(
    int n,
    std::vector<std::vector<double>>& w,
    std::vector<double>& b,
    std::vector<double>& x,
    std::vector<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 = 1; i < n; i++)
    {
        for (j = 0, sum = 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; j < n; j++)
            sum += w[i][j] * x[j];
        x[i] = (x[i] - sum) / w[i][i];
    }
}

bool DirectMethods::GaussianElimimation(
    int n,
    std::vector<std::vector<double>>& w,
    std::vector<double>& b,
    std::vector<double>& x)
    // returns false if matrix is singular
{
    double awikod, col_max, ratio, row_max, temp;
    std::vector<double> d(n);
    int flag = 1, i, i_star, j, k;
    std::vector<int> pivot(n);

    for (i = 0; i < n; i++)
    {
        pivot[i] = i;
        row_max = 0;
        for (j = 0; j < n; j++)
            row_max = fmax(row_max, fabs(w[i][j]));

        if (row_max == 0)
        {
            flag = 0;
            row_max = 1;
        }

        d[i] = row_max;
    }

    if (n <= 1) return flag != 0;

    // factorization
    
    for (k = 0; k < n - 1; k++)
    {
        // determine pivot row the row i_star
        col_max = fabs(w[k][k]) / d[k];
        i_star = k;
        
        for (i = k + 1; i < n; i++)
        {
            awikod = fabs(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;

    if (flag == 0)
        return false;

    Substitute(n, w, b, x, pivot);
    return true;
}

void DirectMethods::LUDecomposition(
    int n,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x)
{
    std::vector<double> y(n);
    std::vector<std::vector<double>> l(n);
    std::vector<std::vector<double>> u(n);

    for (int k = 0; k < n; k++)
    {
        l[k].resize(n);
        u[k].resize(n);
    }

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

            for (int p = 0; p <= k - 1; p++)
                sum += l[k][p] * u[p][j];

            u[k][j] = a[k][j] - sum;
        }

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

            for (int p = 0; p < k - 1; p++)
                sum += l[i][p] * u[p][k];

            l[i][k] = (a[i][k] - sum) / u[k][k];
        }
    }

    // forward substitution

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

        for (int j = 0; j < k; j++)
            sum += l[k][j] * y[j];

        y[k] = b[k] - sum;
    }

    // backward substitution

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

        for (int j = k + 1; j < n; j++)
            sum += u[k][j] * x[j];

        x[k] = (y[k] - sum) / u[k][k];
    }
}


#pragma once
#include <vector>

class ClassicalIterativeMethods
{
public:
    static double Jacobi(
        int n,
        int maxIterations,
        double tolerance,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x0);
    static double GaussSeidel(
        int n,
        int maxIterations,
        double tolerance,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x0);
    static double SOR(
        int n,
        int maxIterations,
        double omega,
        double tolerance,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x0);
    static double InnerProduct(
        int n,
        std::vector<double>& a,
        std::vector<double>& b);
    static double GradientMethod(
        int n,
        int maxIterations,
        double tolerance,
        std::vector<std::vector<double>>& a,
        std::vector<double>& b,
        std::vector<double>& x0);
};


#include "ClassicalIterativeMethods.h"
#include <vector>

double ClassicalIterativeMethods::Jacobi(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

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

            for (int j = 0; j != i && j < n; j++)
                sum += a[i][j] * x0[j];

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

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::GaussSeidel(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

        for (int i = 0; i < n; i++)
        {
            double sum0 = 0;
            double sum1 = 0;

            for (int j = 0; j <= i - 1; j++)
                sum0 += a[i][j] * x[j];

            for (int j = i + 1; j < n; j++)
                if (i != n)
                    sum1 += a[i][j] * x0[j];

            x[i] = (b[i] - sum0 - sum1) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::SOR(
    int n,
    int maxIterations,
    double omega,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

        for (int i = 0; i < n; i++)
        {
            double sum0 = 0;
            double sum1 = 0;

            for (int j = 0; j <= i - 1; j++)
                sum0 += a[i][j] * x[j];

            for (int j = i + 1; j < n; j++)
                if (i != n)
                    sum1 += a[i][j] * x0[j];

            x[i] = (1.0 - omega) * x0[i] +
                omega * (b[i] - sum0 - sum1) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::InnerProduct(
    int n, 
    std::vector<double>& a,
    std::vector<double>& b)
{
    double sum = 0.0;

    for (int i = 0; i < n; i++)
        sum += a[i] * b[i];

    return sum;
}

double ClassicalIterativeMethods::GradientMethod(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> r(n);
    std::vector<double> x(n);
    std::vector<double> y(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

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

            for (int j = 0; j < n; j++)
                sum += a[i][j] * x0[j];

            r[i] = b[i] - sum;
        }

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

            for (int j = 0; j < n; j++)
                sum += a[i][j] * r[j];

            y[i] = sum;
        }

        double alpha = InnerProduct(n, r, r);

        alpha /= InnerProduct(n, r, y);

        for (int i = 0; i < n; i++)
            x[i] = x0[i] + alpha * r[i];

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}


#include "ClassicalIterativeMethods.h"
#include <vector>

double ClassicalIterativeMethods::Jacobi(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

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

            for (int j = 0; j != i && j < n; j++)
                sum += a[i][j] * x0[j];

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

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::GaussSeidel(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

        for (int i = 0; i < n; i++)
        {
            double sum0 = 0;
            double sum1 = 0;

            for (int j = 0; j <= i - 1; j++)
                sum0 += a[i][j] * x[j];

            for (int j = i + 1; j < n; j++)
                if (i != n)
                    sum1 += a[i][j] * x0[j];

            x[i] = (b[i] - sum0 - sum1) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::SOR(
    int n,
    int maxIterations,
    double omega,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> x(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

        for (int i = 0; i < n; i++)
        {
            double sum0 = 0;
            double sum1 = 0;

            for (int j = 0; j <= i - 1; j++)
                sum0 += a[i][j] * x[j];

            for (int j = i + 1; j < n; j++)
                if (i != n)
                    sum1 += a[i][j] * x0[j];

            x[i] = (1.0 - omega) * x0[i] +
                omega * (b[i] - sum0 - sum1) / a[i][i];
        }

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}

double ClassicalIterativeMethods::InnerProduct(
    int n, 
    std::vector<double>& a,
    std::vector<double>& b)
{
    double sum = 0.0;

    for (int i = 0; i < n; i++)
        sum += a[i] * b[i];

    return sum;
}

double ClassicalIterativeMethods::GradientMethod(
    int n,
    int maxIterations,
    double tolerance,
    std::vector<std::vector<double>>& a,
    std::vector<double>& b,
    std::vector<double>& x0)
{
    double t = 0.0;
    std::vector<double> r(n);
    std::vector<double> x(n);
    std::vector<double> y(n);
    int its = 0;

    for (int i = 0; i < n; i++)
        x[i] = x0[i];

    while (its < maxIterations)
    {
        its++;

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

            for (int j = 0; j < n; j++)
                sum += a[i][j] * x0[j];

            r[i] = b[i] - sum;
        }

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

            for (int j = 0; j < n; j++)
                sum += a[i][j] * r[j];

            y[i] = sum;
        }

        double alpha = InnerProduct(n, r, r);

        alpha /= InnerProduct(n, r, y);

        for (int i = 0; i < n; i++)
            x[i] = x0[i] + alpha * r[i];

        t = 0.0;

        for (int i = 0; i < n; i++)
            t += pow(x[i] - x0[i], 2.0);

        t = sqrt(t);

        if (t < tolerance)
            break;

        for (int i = 0; i < n; i++)
            x0[i] = x[i];
    }

    return t;
}


#include "Utility.h"
#include <ctime>
#include <cstring>
#include <windows.h>

wchar_t* Utility::ConvertToWide(
    char* input,
    UINT codePage = CP_UTF8)
{
    if (!input) 
    {
        return nullptr;
    }

    // Determine required wide-character size
    int requiredSize = MultiByteToWideChar(
        codePage,
        MB_ERR_INVALID_CHARS,
        input,
        -1,       // Input is null-terminated
        nullptr,
        0
    );

    if (requiredSize == 0) {
        return nullptr;
    }

    // Allocate buffer (requiredSize includes null terminator)
    wchar_t* output = new (std::nothrow) wchar_t[requiredSize];

    if (!output) {
        return nullptr;
    }

    // Perform conversion
    int result = MultiByteToWideChar(
        codePage,
        MB_ERR_INVALID_CHARS,
        input,
        -1,
        output,
        requiredSize
    );

    if (result == 0) {
        delete[] output;
        output = nullptr;
        return output;
    }

    return output;
}

double Utility::Factorial(int n)
{
    if (n <= 1)
        return 1.0;
    else
        return n * Factorial(n - 1);
}

double Utility::Binomial(int m, int n)
{
    return Factorial(m) / Factorial(m - n) / Factorial(n);
}

std::wstring Utility::Format(double x)
{
    wchar_t line[128] = L"";
    std::wstring result = L"";

    if (x > 0.0)
        result += L"+";
    else if (x == 0.0)
        result += L" ";
    else
        result += L"-";

    x = fabs(x);
    swprintf_s(line, 32, L"%14.8e", x);
    return result + line;
}

bool Utility::FormatTime(
    char buffer[], size_t bufferSize, std::time_t t)
{
    bool result = false;
    const char format[] = "%Y-%m-%d %H:%M:%S";
    
    if (buffer == nullptr || format == nullptr || bufferSize == 0)
    {
        return false;
    }

    std::tm local0 = {};
    
    localtime_s(&local0, &t);

    size_t written = std::strftime(buffer, bufferSize, format, &local0);
    result = written > 0;

    if (!result)
        return false;

    return true;
}

void Utility::Matrix(
    int matrix, int n,
    std::vector<std::vector<double>>& a,
    std::vector<std::vector<double>>& c,
    std::vector<double>& b)
{
    if (matrix == 1)
    {
        // Pascal matrix

        for (int i = 0; i < n; i++)
            b[i] = pow(2, i + 1);

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                a[i][j] = Binomial(i + j, j);
    }

    else if (matrix == 2)
    {
        b[0] = -0.5;

        for (int i = 1; i < n - 1; i++)
            b[i] = -1.5;

        b[n - 1] = +0.5;

        for (int i = 0; i < n; i++)
        {
            a[i][i] = -2.0;

            if (i < n - 1)
                a[i][i + 1LL] = 1.0;

            if (i > 0)
                a[i][i - 1LL] = 1.0;
        }
    }

    else if (matrix == 3)
    {
        // Cauchy matrix

        for (int i = 0; i < n; i++)
        {
            double xi = 2 * i;

            for (int j = 0; j < n; j++)
            {
                double yj = 2 * j + 1;

                a[i][j] = 1.0 / (xi + yj);
            }

            b[i] = i + 1;
        }
    }

    else if (matrix == 4)
    {
        // Lehmer matrix

        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
                a[i][j] = fmin(i + 1, j + 1) + fmax(i + 1, j + 1);

            a[i][i] = i + 1;
            b[i] = i;
        }
    }

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            c[i][j] = a[i][j];
}


// LinearSystems.cpp : Defines the entry point for the application.
//

#include "framework.h"
#include "LinearSystems.h"
#include "DirectMethods.h"
#include "ClassicalIterativeMethods.h"
#include "Utility.h"

#define MAX_LOADSTRING 100

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
WCHAR line[256], text[16384];                   // wide character buffers

// Forward declarations of functions included in this code module:
ATOM                MyRegisterClass(HINSTANCE hInstance);
BOOL                InitInstance(HINSTANCE, int);
LRESULT CALLBACK    WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    About(HWND, UINT, WPARAM, LPARAM);

int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    UNREFERENCED_PARAMETER(hPrevInstance);
    UNREFERENCED_PARAMETER(lpCmdLine);

    // TODO: Place code here.

    // Initialize global strings
    LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
    LoadStringW(hInstance, IDC_LINEARSYSTEMS, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_LINEARSYSTEMS));

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, nullptr, 0, 0))
    {
        if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
        {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }

    return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
    WNDCLASSEXW wcex = { 0 };

    wcex.cbSize = sizeof(WNDCLASSEX);

    wcex.style          = CS_HREDRAW | CS_VREDRAW;
    wcex.lpfnWndProc    = WndProc;
    wcex.cbClsExtra     = 0;
    wcex.cbWndExtra     = 0;
    wcex.hInstance      = hInstance;
    wcex.hIcon          = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_LINEARSYSTEMS));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_LINEARSYSTEMS);
    wcex.lpszClassName  = szWindowClass;
    wcex.hIconSm        = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

    return RegisterClassExW(&wcex);
}

//
//   FUNCTION: InitInstance(HINSTANCE, int)
//
//   PURPOSE: Saves instance handle and creates main window
//
//   COMMENTS:
//
//        In this function, we save the instance handle in a global variable and
//        create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   hInst = hInstance; // Store instance handle in our global variable

   HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    static char ascLine[256] = "";
    static wchar_t* unicodePtr = nullptr;
    static int matrix = 0, method = 0;
    static HWND hEditN = NULL;
    static HWND hEditMultiline = NULL;
    static HFONT hFont = NULL;
    static WCHAR matrixStr[32] = L"";
    static std::time_t t0 = {}, t1 = {};

    switch (message)
    {
    case WM_CREATE:
        CreateWindowEx(0, L"STATIC", L"n", WS_CHILD | WS_VISIBLE,
            10, 10, 60, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
        hEditN = CreateWindowEx(
            WS_EX_CLIENTEDGE,                       // Extended style for sunken border
            TEXT("EDIT"),                           // Class name
            TEXT("2"),                              // Initial text (can be blank)
            WS_CHILD | WS_VISIBLE | ES_LEFT | ES_AUTOHSCROLL,
            100, 10, 60, 20,                         // Position and size
            hWnd,                                   // Parent window handle
            (HMENU)IDC_EDIT_N,                      // Unique control ID
            hInst,                                  // Application instance
            NULL                                    // Extra parameter
        );

        hFont = CreateFont(
            16,                // Height of font
            0,                 // Width of font (0 = default)
            0,                 // Escapement angle
            0,                 // Orientation angle
            FW_BOLD,           // Weight (FW_NORMAL, FW_BOLD, etc.)
            FALSE,             // Italic
            FALSE,             // Underline
            FALSE,             // Strikeout
            UNICODE,           // Character set
            OUT_DEFAULT_PRECIS,// Output precision
            CLIP_DEFAULT_PRECIS,// Clipping precision
            DEFAULT_QUALITY,   // Output quality
            FIXED_PITCH | FF_MODERN, // Pitch and family
            TEXT("Courier New")// Typeface name
        );

        hEditMultiline = CreateWindowEx(
            WS_EX_CLIENTEDGE,                       // Extended style for sunken border
            TEXT("EDIT"),                           // Class name
            TEXT(""),                               // Initial text (can be blank)
            WS_CHILD | WS_VISIBLE | WS_VSCROLL | ES_AUTOHSCROLL |
            ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL | WS_HSCROLL | WS_VSCROLL,
            10, 50, 990, 400,                      // Position and size
            hWnd,                                   // Parent window handle
            (HMENU)IDC_EDIT_MULTILINE,              // Unique control ID
            hInst,                                  // Application instance
            NULL                                    // Extra parameter
        );

        CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            200, 10, 80, 30, hWnd, (HMENU)IDC_COMPUTE_BUTTON, hInst, NULL);
        CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            300, 10, 80, 30, hWnd, (HMENU)IDC_CANCEL_BUTTON, hInst, NULL);
        CreateWindowEx(0, L"BUTTON", L"Clear", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            400, 10, 80, 30, hWnd, (HMENU)IDC_CLEAR_BUTTON, hInst, NULL);

        SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, 0);
        return 0;
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);
            // Parse the menu selections:
            switch (wmId)
            {
            case IDM_ABOUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
                break;
            case IDM_CAUCHY:
                matrix = 1;
                break;
            case IDM_LEHMER:
                matrix = 2;
                break;
            case IDM_OTHER:
                matrix = 3;
                break;
            case IDM_PASCAL:
                matrix = 4;
                break;
            case IDM_GAUSSIAN:
                method = 1;
                break;
            case IDM_GAUSS_SEIDEL:
                method = 2;
                break;
            case IDM_GRADIENT:
                method = 3;
                break;
            case IDM_JACOBI:
                method = 4;
                break;
            case IDM_LUD:
                method = 5;
                break;
            case IDM_SOR:
                method = 6;
                break;
            case IDC_COMPUTE_BUTTON:
            {
                t0 = time(NULL);
                int n = GetDlgItemInt(hWnd, IDC_EDIT_N, FALSE, FALSE);
                int maxIterations = 20000 * n;
                double omega = 1.5;
                double tolerance = n <= 5 ? 1.0e-12 : 1.0e-8;
                std::vector<double> b(n);
                std::vector<double> x(n);
                std::vector<std::vector<double>> a(n);
                std::vector<std::vector<double>> c(n);

                for (int i = 0; i < n; i++)
                {
                    a[i].resize(n);
                    c[i].resize(n);
                }
                
                if (matrix == 4)
                {
                    // Pascal matrix

                    for (int i = 0; i < n; i++)
                        b[i] = pow(2, i + 1);

                    for (int i = 0; i < n; i++)
                        for (int j = 0; j < n; j++)
                            a[i][j] = Utility::Binomial(i + j, j);

                    wcscpy_s(matrixStr, L"Matrix Type Pascal");
                }

                else if (matrix == 3)
                {
                    wcscpy_s(matrixStr, L"Matrix Type Other");

                    b[0] = -0.5;

                    for (int i = 1; i < n - 1; i++)
                        b[i] = -1.5;

                    b[n - 1LL] = +0.5;

                    for (int i = 0; i < n; i++)
                    {
                        a[i][i] = -2.0;

                        if (i < n - 1)
                            a[i][i + 1LL] = 1.0;

                        if (i > 0)
                            a[i][i - 1LL] = 1.0;
                    }
                }

                else if (matrix == 1)
                {
                    // Cauchy matrix

                    for (int i = 0; i < n; i++)
                    {
                        double xi = 2 * i;

                        for (int j = 0; j < n; j++)
                        {
                            double yj = 2 * j + 1;

                            a[i][j] = 1.0 / (xi + yj);
                        }

                        b[i] = i + 1;
                    }

                    wcscpy_s(matrixStr, L"Matrix Type Cauchy");
                }

                else if (matrix == 2)
                {
                    // Lehmer matrix

                    for (int i = 0; i < n; i++)
                    {
                        for (int j = 0; j < n; j++)
                            a[i][j] = fmin(i + 1, j + 1) + fmax(i + 1, j + 1);

                        a[i][i] = i + 1;
                        b[i] = i;
                    }

                    wcscpy_s(matrixStr, L"Matrix Type Lehmer");
                }

                for (int i = 0; i < n; i++)
                    for (int j = 0; j < n; j++)
                        c[i][j] = a[i][j];

                text[0] = L'\0';
                wcscpy_s(text, matrixStr);
                wcscat_s(text, L"\r\n");

                switch (method)
                {
                case 1:
                    wcscat_s(text, L"Gaussian Elimination:\r\n");
                    DirectMethods::GaussianElimimation(n, a, b, x);
                    break;
                case 2:
                    wcscat_s(text, L"Gauss-Seidel Method:\r\n");
                    ClassicalIterativeMethods::GaussSeidel(
                        n, maxIterations, tolerance, a, b, x);
                    break;
                case 3:
                    wcscat_s(text, L"Gradient Method:\r\n");
                    ClassicalIterativeMethods::GradientMethod(
                        n, maxIterations, tolerance, a, b, x);
                    break;
                case 4:
                    wcscat_s(text, L"Jacobi Method:\r\n");
                    ClassicalIterativeMethods::Jacobi(
                        n, maxIterations, tolerance, a, b, x);
                    break;
                case 5:
                    wcscat_s(text, L"LU Decomposition:\r\n");
                    DirectMethods::LUDecomposition(n, a, b, x);
                    break;
                case 6:
                    wcscat_s(text, L"Successive Over Relaxation:\r\n");
                    ClassicalIterativeMethods::SOR(
                        n, maxIterations, omega, tolerance, a, b, x);
                    break;
                }

                for (int i = 0; i < n; i++)
                {
                    std::wstring format = Utility::Format(x[i]) + L"\r\n";
                    wcscat_s(text, format.c_str());
                }

                double error = 0.0;

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

                    for (int j = 0; j < n; j++)
                        sum += c[i][j] * x[j];

                    error += pow(b[i] - sum, 2.0);
                }

                error = sqrt(error / n);
                swprintf_s(line, 256, L"RMS Error = %14.8e\r\n", error);
                wcscat_s(text, line);

                t1 = time(NULL);
                Utility::FormatTime(ascLine, 256, t0);
                unicodePtr = Utility::ConvertToWide(ascLine, CP_UTF8);
                wcscat_s(text, unicodePtr);
                wcscat_s(text, L"\r\n");
                Utility::FormatTime(ascLine, 256, t1);
                unicodePtr = Utility::ConvertToWide(ascLine, CP_UTF8);
                wcscat_s(text, unicodePtr);
                wcscat_s(text, L"\r\n");
                SetWindowText(hEditMultiline, text);
                return 1;
            }
            case IDC_CANCEL_BUTTON:
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            case IDC_CLEAR_BUTTON:
                text[0] = L'\0';
                SetWindowText(hEditMultiline, text);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
        {
            PAINTSTRUCT ps;
            HDC hdc = BeginPaint(hWnd, &ps);
            // TODO: Add any drawing code that uses hdc here...
            EndPaint(hWnd, &ps);
        }
        break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

// Message handler for about box.

INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

Blog Entry © Wednesday, December 24, 2025, by James Pate Williams, Jr. ID3 Decision Tree Metadata Parser

Blog Entry Wednesday December 24 2025Download
// ID3MetadataParser.cpp (c) December 2025
// by James Pate Williams, Jr.

#include "pch.h"

#define FILE_EOF			0
#define NO_ERROR			1
#define EMPTY_FILE			2
#define INVALID_LINE		3
#define MISSING_NAME		4
#define INVALID_NAME		5
#define INVALID_DESCRIPTION 6
#define MISSING_DESCRIPTION 7
#define INVALID_TYPE		8
#define MISSING_TYPE		9
#define INVALID_RANGE		10
#define INVALID_CATEGORICAL	11
#define INVALID_DOUBLE		12
#define INVALID_FLOAT		13
#define INVALID_INTEGER		14
#define INVALID_ROLE		15
#define MISSING_ROLE		16

enum AttributeType {
	categorical, integer, doubleReal, FloatReal
};

typedef struct tagCategoricalAttribute {
	std::string name = "";
	std::string description = "";
	std::vector<char> category;
} CategoricalAttribute, * PCategoricalAttribute;

typedef struct tagIntegerAttribute {
	std::string name = "";
	std::string description = "";
	int loValue= -1, hiValue = -1;
} IntegerAttribute, * PIntegerAttribute;

typedef struct tagDoubleAttribute {
	std::string name = "";
	std::string description = "";
	double loValue = -1.0, hiValue = -1.0;
} DoubleAttribute, * PDoubleAttribute;

typedef struct tagFloatAttribute {
	std::string name = "";
	std::string description = "";
	float loValue = -1.0f, hiValue = -1.0f;
} FloatAttribute, * PFloatAttribute;

static bool parseName(
	char line[],
	int length,
	int& errorCode,
	int& index,
	bool& feature,
	std::string& name)
{
	char* cptr1 = std::strstr(line, "#name: feature ");
	char* cptr2 = std::strstr(line, "#name: target ");

	if (cptr1 == nullptr && cptr2 == nullptr) {
		errorCode = MISSING_NAME;
		return false;
	}

	if (cptr1) {
		feature = true;
		index = static_cast<int>(strlen("#name: feature "));
	}

	else if (cptr2) {
		feature = false;
		index = static_cast<int>(strlen("#name: target "));
	}

	else {
		errorCode = INVALID_NAME;
		return false;
	}

	if (index >= static_cast<int>(strlen(line))) {
		errorCode = INVALID_NAME;
		return false;
	}

	if (line[index] >= L'A' && line[index] <= 'Z' ||
		line[index] >= L'a' && line[index] <= 'z') {
		bool first = true;

		name = "";

		while (index < strlen(line)) {
			if (line[index] >= 'A' && line[index] <= 'Z' ||
				line[index] >= 'a' && line[index] <= 'z' ||
				line[index] == ' ') {
				if (first)
					name += line[index++];
				else if (first &&
					line[index] >= '0' &&
					line[index] <= '9') {
					first = false;
					name += line[index++];
				}

				if (!first)
					name += line[index++];
			}

			else if (!first) {
				errorCode = INVALID_NAME;
				return false;
			}
		}
	}

	errorCode = 0;
	index = length;
	return true;
}

static bool parseDescription(
	char line[],
	int length,
	int& errorCode,
	int& index,
	std::string& description) {
	
	char* cptr = std::strstr(line, "#description: ");

	if (cptr == nullptr) {
		errorCode = MISSING_DESCRIPTION;
		return false;
	}

	int lengthDesc = static_cast<int>(
		strlen("#description: "));

	if (lengthDesc == length) {
		errorCode = INVALID_DESCRIPTION;
		return false;
	}

	index = lengthDesc;

	while (index < length)
		description += line[index++];
	
	errorCode = 0;
	return true;
}

static bool parseCategorical(
	char line[],
	int length,
	int& errorCode,
	int& index,
	std::vector<char>& category) {
	int delta = static_cast<int>(strlen("#type: categorical: {"));
	char* cptr = line + delta - 1;
	char ch = *cptr++;

	while (ch != '}' && index < length) {
		while (ch != ',' && index < length) {
			
			if (ch == '}') {
				if (index == length - 1)
					break;
				
				else {
					errorCode = INVALID_TYPE;
					return false;
				}

			}
			
			category.push_back(ch);
			index++;
			break;
		}

		cptr++;
		ch = *cptr;
	}

	if (category.size() != 0 && ch == '}') {
		errorCode = 0;
		return true;
	}

	else {
		errorCode = INVALID_CATEGORICAL;
		return false;
	}
}

static bool parseDoubleRange(
	char line[],
	int length,
	int& errorCode,
	int& index,
	double& hiDouble,
	double& loDouble)
{
	index = static_cast<int>(strlen("#type: doubleReal ["));
	char ch = line[index++];
	std::string doubleStr;

	while (ch != ',' &&
		index < static_cast<int>(strlen(line))) {
		doubleStr.push_back(ch);
		ch = line[index++];
	}

	if (doubleStr.size() == 0) {
		errorCode = INVALID_DOUBLE;
		return false;
	}

	try {
		loDouble = std::stod(doubleStr);
		doubleStr = "";
		ch = line[index++];

		while (ch != ']' && index < strlen(line)) {
			doubleStr.push_back(ch);
			ch = line[index++];
		}

		if (doubleStr.size() == 0) {
			errorCode = INVALID_DOUBLE;
			return false;
		}

		hiDouble = std::stod(doubleStr);
		errorCode = 0;
		return true;
	}
	catch (std::exception ex) {
		errorCode = INVALID_DOUBLE;
		return false;
	}

	errorCode = INVALID_RANGE;
	return false;
}

static bool parseFloatRange(
	char line[],
	int length,
	int& errorCode,
	int& index,
	float& hiFloat,
	float& loFloat)
{
	char ch = '\0';
	std::string floatStr;
	ch = line[index++];

	while (ch != ',' && index < strlen(line)) {
		floatStr.push_back(ch);
		ch = line[index++];
	}

	if (floatStr.size() == 0) {
		errorCode = INVALID_INTEGER;
		return false;
	}

	else {
		try {
			loFloat = std::stof(floatStr);
			floatStr = "";
			ch = line[++index];

			while (ch != ']' && index < strlen(line)) {
				floatStr.push_back(ch);
				ch = line[index++];
			}

			if (floatStr.size() == 0) {
				errorCode = INVALID_FLOAT;
				return false;
			}

			hiFloat = std::stof(floatStr);
			errorCode = 0;
			return true;
		}
		catch (std::exception ex) {
			errorCode = INVALID_FLOAT;
			return false;
		}
	}

	errorCode = INVALID_RANGE;
	return false;
}

static bool parseIntegerRange(
	char line[],
	int length,
	int& errorCode,
	int& index,
	int& hiInteger,
	int& loInteger) {
	char ch = '\0';
	int i = 0;
	std::string integerStr;
	
	ch = line[i++];

	if (ch < '0' || ch > '9') {
		errorCode = INVALID_INTEGER;
		return false;
	}

	while (ch != ',' &&	index < length) {
		integerStr.push_back(ch);
		ch = line[i++];
		index++;
	}

	if (integerStr.size() == 0) {
		errorCode = INVALID_INTEGER;
		return false;
	}

	else {
		try {
			loInteger = std::stoi(integerStr);
			integerStr = "";
			i = 0;
			ch = line[i++];
			ch = line[i++];
			ch = line[i++];
			index += 3;

			while (
				ch >= '0' && ch <= '9' &&
				ch != ']' && index < length) {
				integerStr.push_back(ch);
				ch = line[i++];
				index++;
			}

			if (integerStr.size() == 0) {
				errorCode = INVALID_INTEGER;
				return false;
			}

			hiInteger = std::stoi(integerStr);
			errorCode = 0;
			return true;
		}
		catch (std::exception ex) {
			errorCode = INVALID_INTEGER;
			return false;
		}
	}

	errorCode = INVALID_RANGE;
	return false;
}

static bool parseType(
	char line[],
	int length,
	double& hiDouble,
	double& loDouble,
	float& hiFloat,
	float& loFloat,
	int& errorCode,
	int& index,
	int& hiInteger,
	int& loInteger,
	std::string& type,
	std::vector<char>& alphabet) {

	char* cptr = std::strstr(line, "#type: ");

	if (cptr == nullptr) {
		errorCode = MISSING_TYPE;
		return false;
	}

	int lengthType = static_cast<int>(strlen("#type: "));

	if (lengthType >= length) {
		errorCode = INVALID_TYPE;
		return false;
	}

	index = lengthType;
	cptr = line + index;

	if (std::strstr(cptr, "categorical {") != nullptr) {
		if (parseCategorical(line, length, errorCode,
			index, alphabet)) {
			errorCode = 0;
			type = "categorical";
			return true;
		}

		else {
			errorCode = INVALID_CATEGORICAL;
			return false;
		}
	}

	if (std::strstr(cptr, "integer [") != nullptr) {
		bool pir = parseIntegerRange(
			line + index + strlen("integer ["),
			length,
			errorCode,
			index,
			hiInteger,
			loInteger);
		if (pir) {
			type = "integer";
			return true;
		}

		else
			return false;
	}

	if (std::strstr(cptr, "doubleReal [") != nullptr) {
		bool pdr = parseDoubleRange(
			line,
			length,
			errorCode,
			index,
			hiDouble,
			loDouble);

		if (pdr) {
			type = "doubleReal";
			return true;
		}

		else
			return false;
	}

	if (std::strstr(cptr, "floatReal [") != nullptr) {
		bool pfr = parseFloatRange(
			line,
			length,
			errorCode,
			index,
			hiFloat,
			loFloat);

		if (pfr) {
			type = "floatReal";
			return true;
		}
	}

	errorCode = INVALID_TYPE;
	return false;
}

static bool readMetaDataLine(
	std::ifstream& file1,
	char line[],
	int& errorCode,
	int& index) {
	file1.getline(line, 256);

	if (strlen(line) == 0 && index == -1) {
		errorCode = EMPTY_FILE;
		return false;
	}

	if (file1.eof()) {
		errorCode = 0;
		index = FILE_EOF;
		return false;
	}

	if (strlen(line) > 0 &&
		std::strstr(line, "#endheader") != nullptr)
		return false;

	if (strlen(line) > 0)
		return true;
	else
		return false;
}

double dbl_max[8] = { 0 };
double dbl_min[8] = { 0 };
int int_max = 0;
int int_min = 0;

static void readDatasetFile(
	std::ifstream& file2) {
	char line[256] = "";

	for (int i = 0; i < 8; i++) {
		dbl_min[i] = DBL_MAX;
		dbl_max[i] = DBL_MIN;
	}

	int_min = INT_MAX;
	int_max = INT_MIN;

	while (!file2.eof()) {
		file2.getline(line, 256);
		int count = 0, index = 0;

		while (
			count <= 9 &&
			index < static_cast<int>(strlen(line))) {
			char ch = line[index++], subline[256] = "";
			int cp = 0;

			while (ch != ',' && cp < static_cast<int>(strlen(line))) {
				subline[cp++] = ch;
				ch = line[index++];
			}

			count++;

			if (strlen(subline) >= 1)
				subline[cp] = '\0';

			if (count >= 1 && count <= 8 && cp > 1) {
				std::string substr(subline);
				double x = std::stod(subline);

				if (x > dbl_max[count - 1])
					dbl_max[count - 1] = x;
				if (x < dbl_min[count - 1])
					dbl_min[count - 1] = x;
			}

			else if (count == 9 && !(
				strstr(subline, "F") ||
				strstr(subline, "I") ||
				strstr(subline, "M"))) {
				std::string substr(subline);
				int x = std::stoi(substr);

				if (x > int_max)
					int_max = x;
				if (x < int_min)
					int_min = x;
			}
		}
	}

	file2.close();
}

int main()
{
	bool feature = false;
	char filename1[256] = "C:\\Users\\James\\OneDrive\\Desktop\\ID3MetadataParser\\x64\\Debug\\ID3MetadataParserDataFile.txt";
	char filename2[256] = "C:\\Users\\James\\OneDrive\\Desktop\\ID3MetadataParser\\x64\\Debug\\abalone.data.txt";
	char line[256] = "";
	int errorCode = -1, index = -1, role = -1;
	std::ifstream file1(filename1);
	std::ifstream file2(filename2);

	// file1 format
	std::string name, description, type;
	std::vector<char> category;
	
	// file2 format
	std::string cat, length, diameter, height, whole;
	std::string shucked, viscera, shell, rings;
	
	std::vector<CategoricalAttribute> categoricalAttributes;
	std::vector<IntegerAttribute> integerAttributes;
	std::vector<DoubleAttribute> doubleAttributes;
	std::vector<FloatAttribute> floatAttributes;

	std::vector<std::string> names;
	std::vector<std::string> descriptions;
	std::vector<std::string> types;

	while (!file1.eof()) {
		index = -1;

		bool result = readMetaDataLine(
			file1,
			line,
			errorCode,
			index);

		if (!result)
			break;

		index = 0;
		int length = static_cast<int>(strlen(line));

		if (length == 0)
			break;
		
		name = "";

		bool pn = parseName(
			line,
			length,
			errorCode,
			index,
			feature,
			name);

		if (pn) {
			bool result = readMetaDataLine(
				file1,
				line,
				errorCode,
				index);

			if (!result)
				break;
		
			length = static_cast<int>(strlen(line));
			index = 0;
			description = "";

			bool pd = parseDescription(
				line,
				length,
				errorCode,
				index,
				description);

			if (pd) {
				bool result = readMetaDataLine(
					file1,
					line,
					errorCode,
					index);

				if (!result)
					break;

				length = static_cast<int>(strlen(line));
				index = 0;
				type = "";

				double hiDouble = DBL_MIN;
				double loDouble = DBL_MAX;
				float hiFloat = FLT_MIN;
				float loFloat = FLT_MAX;
				int hiInteger = INT_MIN;
				int loInteger = INT_MAX;

				bool pt = parseType(
					line,
					length,
					hiDouble,
					loDouble,
					hiFloat,
					loFloat,
					errorCode,
					index,
					hiInteger,
					loInteger,
					type,
					category);
					length = static_cast<int>(strlen(line));

				if (pt) {
					if (type == "categorical") {
						CategoricalAttribute ca;
						ca.category = category;
						ca.description = description;
						ca.name = name;
						categoricalAttributes.push_back(ca);
					}

					else if (type == "integer") {
						IntegerAttribute ia;
						ia.loValue = loInteger;
						ia.hiValue = hiInteger;
						ia.description = description;
						ia.name = name;
						integerAttributes.push_back(ia);
					}

					else if (type == "doubleReal") {
						DoubleAttribute da;
						da.loValue = loDouble;
						da.hiValue = hiDouble;
						da.description = description;
						da.name = name;
						doubleAttributes.push_back(da);
					}

					else if (type == "floatReal") {
						FloatAttribute fa;
						fa.loValue = loFloat;
						fa.hiValue = hiFloat;
						fa.description = description;
						fa.name = name;
						floatAttributes.push_back(fa);
					}

					else {
						errorCode = INVALID_TYPE;
						break;
					}
				}

				else {
					errorCode = MISSING_TYPE;
					break;
				}
			}

			else {
				errorCode = INVALID_DESCRIPTION;
				break;
			}
		}

		else {
			errorCode = INVALID_NAME;
			return false;
		}
	}

	readDatasetFile(file2);

	for (int i = 1; i <= 7; i++) {
		std::cout << i << '\t' << dbl_min[i];
		std::cout << '\t' << dbl_max[i];
		std::cout << std::endl;
	}

	std::cout << "8\t" << int_min << '\t' << int_max;
	std::cout << std::endl;
	std::cout << std::endl;

	for (int i = 0; i < static_cast<int>(categoricalAttributes.size()); i++) {
		std::cout << categoricalAttributes[i].name << ' ';
		std::cout << categoricalAttributes[i].description << ' ';
		std::cout << std::endl;
	}

	for (int i = 0; i < static_cast<int>(doubleAttributes.size()); i++) {
		std::cout << doubleAttributes[i].name << ' ';
		std::cout << doubleAttributes[i].description << ' ';
		std::cout << doubleAttributes[i].loValue << ' ';
		std::cout << doubleAttributes[i].hiValue;
		std::cout << std::endl;
	}

	for (int i = 0; i < static_cast<int>(floatAttributes.size()); i++) {
		std::cout << floatAttributes[i].name << ' ';
		std::cout << floatAttributes[i].description << ' ';
		std::cout << floatAttributes[i].loValue << ' ';
		std::cout << floatAttributes[i].hiValue;
		std::cout << std::endl;
	}

	for (int i = 0; i < static_cast<int>(integerAttributes.size()); i++) {
		std::cout << integerAttributes[i].name << ' ';
		std::cout << integerAttributes[i].description << ' ';
		std::cout << integerAttributes[i].loValue << ' ';
		std::cout << integerAttributes[i].hiValue;
		std::cout << std::endl;
	}

	file1.close();
	return 0;
}

Blog Entry (c) Wednesday, August 13, 2025, by James Pate Williams, Jr. Exercises from an Online Textbook

Exercises Szabo OstlundDownload
#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 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 std::complex<double> CmpDeterminant(
        bool& failure, int n,
        std::vector<std::vector<std::complex<double>>>& A);
    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);
};

#include <vector>

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(
        bool& failure, int n,
        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);
};

// 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.

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

int main()
{
	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 } };
	int m = 3, n = 3, p = 3;

	std::vector<double> br(3);
	std::vector<size_t> pivot(3);
	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));
	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));
	
	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];
		}
	}

	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];
		}
	}
	
	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;
	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
	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 = " << 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;
	std::cout << "Ac" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
	std::cout << std::endl;
	inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
	std::cout << "Ac inverse = " << 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);
}

Blog Entry Wednesday, July 10, 2024, © James Pate Williams, Jr. My Dual Interests in Cryptography and Number Theory

I became fascinated with secret key cryptography as a child. Later, as an adult, in around 1979, I started creating crude symmetric cryptographic algorithms. I became further enthralled with cryptography and number theory in 1996 upon reading Applied Cryptography, Second Edition: Protocols, Algorithms, and Source Code in C by Bruce Schneier and later the Handbook of Applied Cryptography by Alfred J. Menezes, Paul C. van Oorschot, and Scott A. Vanstone. After implementing many of the algorithms in both tomes, I communicated my results to two of the authors namely Bruce Schneier and Professor Alfred J. Menezes. In 1997 I developed a website devoted to constraint satisfaction problems and their solutions, cryptography, and number theory. I posted legal C and C++ source code. Professor Menezes advertised my website along with his treatise. See the following blurb:

In the spirit of my twin scientific infatuations, I offer yet another C integer factoring implementation utilizing the Free Large Integer Package (known more widely as lip) which was created by Arjen K. Lenstra (now a Professor Emeritus). This implementation includes Henri Cohen’s Trial Division algorithm, the Brent-Cohen-Pollard rho method, the Cohen-Pollard p – 1 stage 1 method, and the Lenstra lip Elliptic Curve Method. If I can get the proper authorization, I will later post the source code.

total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
1
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.014000 seconds:
2
5 ^ 2
41
397
2113
enter number below:
0

total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
2
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 1.531000 seconds:
2
5 ^ 2
41
397
2113
415878438361
3630105520141
enter number below:
0

total time required for initialization: 0.055000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
3
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.066000 seconds:
2
5 ^ 2
41
838861
415878438361
3630105520141
enter number below:
0

total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
4
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.013000 seconds:
2
5
205
838861
415878438361
3630105520141
enter number below:
0

Blog Entry Sunday, June 23, 2024 (c) James Pate Williams, Jr.

The object of this C Win32 application is to find a multiple of 9 with its digits summing to a multiple of 9 also. The first column below is a multiple of 9 whose digits sum to 9 also. The second column is the sum of digits found in the column one number. The last column is the first column divided by 9.

Enter PRNG seed:
1
Enter number of bits (4 to 16):
4
    9       9       1
Enter number of bits (4 to 16):
5
   27       9       3
Enter number of bits (4 to 16):
6
   45       9       5
Enter number of bits (4 to 16):
7
  117       9      13
Enter number of bits (4 to 16):
8
  252       9      28
Enter number of bits (4 to 16):
0

C:\Users\james\source\repos\CProductOf9Console\Debug\CProductOf9Console.exe (process 23280) exited with code 0.
Press any key to close this window . . .
Enter PRNG seed:
1
Enter number of bits (4 to 16):
9
  369      18      41
Enter number of bits (4 to 16):
10
  846      18      94
Enter number of bits (4 to 16):
11
 1080       9     120
Enter number of bits (4 to 16):
12
 3015       9     335
Enter number of bits (4 to 16):
13
 5040       9     560
Enter number of bits (4 to 16):
14
10350       9    1150
Enter number of bits (4 to 16):
15
30870      18    3430
Enter number of bits (4 to 16):
16
57798      36    6422
Enter number of bits (4 to 16):
0
// CProductOf9Console.c (c) Sunday, June 23, 2024
// by James Pate Williams, Jr., BA, BS, MSwE, PhD

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

char nextStr[256], numbStr[256];

void ConvertToString(int number, int radix)
{
	int i = 0;

	while (number > 0)
	{
		nextStr[i++] = (char)(number % radix + '0');
		number /= radix;
	}

	nextStr[i++] = '\0';
	_strrev(nextStr);
}

int Sum(int next)
{
	long sum = 0;

	ConvertToString(next, 10);

	for (int i = 0; i < (int)strlen(nextStr); i++)
		sum += (long)nextStr[i] - '0';

	if (sum % 9 == 0 && sum != 0)
		return sum;

	return -1;
}

long GetNext(int numBits, int* next)
{
	long hi = 0, lo = 0, nine = 0;

	nextStr[0] = '\0';
	numbStr[0] = '\0';

	if (numBits == 4)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 16);

			if (*next != 0 && *next >= 8 && *next < 16)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 5)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 32);

			if (*next >= 16 && *next < 32)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 6)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 64);

			if (*next >= 32 && *next < 64)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 7)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 128);

			if (*next >= 64 && *next < 128)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 8)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 256);

			if (*next >= 128 && *next < 256)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 9)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 512);

			if (*next >= 256 && *next < 512)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 10)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 1024);

			if (*next >= 512 && *next < 1024)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 11)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 2048);

			if (*next >= 1024 && *next < 2048)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 12)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 4096);

			if (*next >= 2048 && *next < 4096)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 13)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 8192);

			if (*next >= 4096 && *next < 8192)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 14)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 16384);

			if (*next >= 8192 && *next < 16384)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 15)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 32768);

			if (*next >= 16384 && *next < 32768)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	else if (numBits == 16)
	{
		while (1)
		{
			*next = 9 * (long)(rand() % 65536);

			if (*next >= 32768 && *next < 65536)
			{
				nine = Sum(*next);

				if (nine % 9 == 0)
					return nine;
			}
		}
	}

	return -1;
}

int main()
{
	char buffer[256] = { '\0' };
	long seed = 0;

	printf_s("Enter PRNG seed:\n");
	scanf_s("%s", buffer, sizeof(buffer));
	seed = atol(buffer);
	srand((unsigned int)seed);

	while (1)
	{
		int next = 0, nine = 0, numberBits = 0;

		printf_s("Enter number of bits (4 to 16):\n");
		scanf_s("%s", buffer, sizeof(buffer));
		numberBits = atol(buffer);

		if (numberBits == 0)
			break;

		if (numberBits < 4 || numberBits > 16)
		{
			printf_s("illegal number of bits must >= 4 and <= 16\n");
			continue;
		}

		nine = GetNext(numberBits, &next);

		if (nine == -1)
		{
			printf_s("illegal result, try again\n");
			continue;
		}

		printf_s("%5ld\t%5ld\t%5ld\n", next, nine, next / 9);
	}

	return 0;
}

Blog Entry (c) Friday, June 21, 2024, by James Pate Williams, Jr. Comparison of Two Prime Number Sieves

First the C++ results:

Limit = 1000000
Number of primes <= 1000000 78498
Milliseconds taken by Sieve of Atkin: 12
Number of primes <= 1000000 78498
Milliseconds taken by Sieve of Eratosthenes: 14
Limit = 10000000
Number of primes <= 10000000 664579
Milliseconds taken by Sieve of Atkin: 159
Number of primes <= 10000000 664579
Milliseconds taken by Sieve of Eratosthenes: 204
Limit = 100000000
Number of primes <= 100000000 5761455
Milliseconds taken by Sieve of Atkin: 1949
Number of primes <= 100000000 5761455
Milliseconds taken by Sieve of Eratosthenes: 2343
Limit = 0

Next, we have the Java results:

C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 1000000 0
number of primes less than equal 1000000 = 78498
total computation time in seconds = 0.008

C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 10000000 0
number of primes less than equal 10000000 = 664579
total computation time in seconds = 0.098

C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 1000000 0
number of primes less than equal 1000000 = 78498
total computation time in seconds = 0.011

C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 10000000 0
number of primes less than equal 10000000 = 664579
total computation time in seconds = 0.151

C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 100000000 0
number of primes less than equal 100000000 = 5761455
total computation time in seconds = 1.511

C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 100000000 0
number of primes less than equal 100000000 = 5761455
total computation time in seconds = 1.995

Notice that the Java application outperforms the C++ application.

// PrimeSieveComparison.cpp (c) Friday, June 21, 2024
// by James Pate Williams, Jr.
//
//  SieveOfAtkin.java
//  SieveOfAtkin
//
//  Created by James Pate Williams, Jr. on 9/29/07.
//  Copyright (c) 2007 James Pate Williams, Jr. All rights reserved.
//
//  SieveOfEratosthenes.java
//  SieveOfEratosthenes
//
//  Created by James Pate Williams, Jr. on 9/29/07.
//  Copyright (c) 2007 James Pate Williams, Jr. All rights reserved.
//

#include <math.h>
#include <iostream>
#include <chrono>
using namespace std::chrono;
using namespace std;

const int Maximum = 100000000;
bool sieve[Maximum + 1];

void SieveOfAtkin(int limit)
{
	auto start = high_resolution_clock::now();
	int e, k, n, p, x, xx3, xx4, y, yy;
	int primeCount = 2, sqrtLimit = (int)sqrt(limit);

	for (n = 5; n <= limit; n++)
		sieve[n] = false;

	for (x = 1; x <= sqrtLimit; x++) {
		xx3 = 3 * x * x;
		xx4 = 4 * x * x;
		for (y = 1; y <= sqrtLimit; y++) {
			yy = y * y;
			n = xx4 + yy;
			if (n <= limit && (n % 12 == 1 || n % 12 == 5))
				sieve[n] = !sieve[n];
			n = xx3 + yy;
			if (n <= limit && n % 12 == 7)
				sieve[n] = !sieve[n];
			n = xx3 - yy;
			if (x > y && n <= limit && n % 12 == 11)
				sieve[n] = !sieve[n];
		}
	}

	for (n = 5; n <= sqrtLimit; n++) {
		if (sieve[n]) {
			e = 1;
			p = n * n;
			while (true) {
				k = e * p;
				if (k > limit)
					break;
				sieve[k] = false;
				e++;
			}
		}
	}
	
	for (n = 5; n <= limit; n++)
		if (sieve[n])
			primeCount++;

	auto stop = high_resolution_clock::now();
	auto duration = duration_cast<milliseconds>(stop - start);

	std::cout << "Number of primes <= " << limit << ' ';
	std::cout << primeCount << endl;
	std::cout << "Milliseconds taken by Sieve of Atkin: "
		<< duration.count() << endl;
}

void SieveOfEratosthenes(int limit)
{
	auto start = high_resolution_clock::now();
	int i = 0, k = 0, n = 0, nn = 0;
	int primeCount = 0, sqrtLimit = (int)sqrt(limit);

	// initialize the prime number sieve

	for (n = 2; n <= limit; n++)
		sieve[n] = true;

	// eliminate the multiples of n

	for (n = 2; n <= sqrtLimit; n++)
		for (i = 2; i <= n - 1; i++)
			sieve[i * n] = false;

	// eliminate squares

	for (n = 2; n <= sqrtLimit; n++) {
		if (sieve[n]) {
			k = 0;
			nn = n * n;
			i = nn + k * n;
			while (i <= limit) {
				sieve[i] = false;
				i = nn + k * n;
				k++;
			}
		}
	}

	primeCount = 0;

	for (n = 2; n <= limit; n++)
		if (sieve[n])
			primeCount++;

	auto stop = high_resolution_clock::now();
	auto duration = duration_cast<milliseconds>(stop - start);

	std::cout << "Number of primes <= " << limit << ' ';
	std::cout << primeCount << endl;
	std::cout << "Milliseconds taken by Sieve of Eratosthenes: "
		<< duration.count() << endl;
}

int main()
{
	while (true)
	{
		int limit = 0;
		std::cout << "Limit = ";
		cin >> limit;

		if (limit == 0)
			break;

		SieveOfAtkin(limit);
		SieveOfEratosthenes(limit);
	}

	return 0;
}

Factorizations of Some Fibonacci Sequence Numbers, Lucas Sequence Numbers and Some Other Numbers Using Arjen K. Lenstra’s Free Large Integer Package and the Elliptic Curve Method (c) January 28, 2024, by James Pate Williams, Jr.

All of the following computations were performed on a late 2015 Dell XPS 8900 personal computer with a 64-bit Intel Core I7 processor @ 4.0GHz with 16GB of DDR2 RAM.

Factorization of Six Fibonacci Sequence Numbers:

 

Fibonacci 500

# digits 105

5 ^ 2 p # digits 1

15 c # digits 2

101 p # digits 3

401 p # digits 3

1661 c # digits 4

3001 p # digits 4

10291 c # digits 5

570601 p # digits 6

112128001 p # digits 9

1353439001 p # digits 10

28143378001 p # digits 11

5465167948001 p # digits 13

84817574770589638001 p # digits 20

158414167964045700001 p # digits 21

Runtime (s) = 1.206000



Fibonacci 505

# digits 106

5 p # digits 1

743519377 p # digits 9

44614641121 p # digits 11

770857978613 p # digits 12

960700389041 p # digits 12

12588421794766514566269164716286291055826556238643852856601641 p # digits 62

Runtime (s) = 1.959000



Fibonacci 510

# digits 107

2 ^ 3 p # digits 1

11 p # digits 2

61 p # digits 2

1021 p # digits 4

1597 p # digits 4

3469 p # digits 4

3571 p # digits 4

9521 p # digits 4

53551 p # digits 5

95881 p # digits 5

142445 c # digits 6

1158551 p # digits 7

3415914041 p # digits 10

20778644396941 p # digits 14

20862774425341 p # digits 14

81358225616651 c # digits 14

162716451241291 p # digits 15

Runtime (s) = 2.682000



Fibonacci 515

# digits 108

5 p # digits 1

519121 p # digits 6

5644193 p # digits 7

512119709 p # digits 9

84388938382141 p # digits 14

300367026458796424297447559250634818495937628065437243817852436228914621 p # digits 72

Runtime (s) = 7.861000



Fibonacci 520

# digits 109

131 p # digits 3

451 c # digits 3

521 p # digits 3

2081 p # digits 4

2161 p # digits 4

3121 p # digits 4

24571 p # digits 5

90481 p # digits 5

2519895 c # digits 7

21183761 p # digits 8

57089761 p # digits 8

102193207 p # digits 9

1932300241 p # digits 10

14736206161 p # digits 11

5836312049326721 p # digits 16

42426476041450801 p # digits 17

Runtime (s) = 5.155000



Fibonacci 525

# digits 110

2 p # digits 1

5 p # digits 1

421 p # digits 3

701 p # digits 3

3001 p # digits 4

3965 c # digits 4

4201 p # digits 4

141961 p # digits 6

2553601 p # digits 7

230686501 p # digits 9

8288823481 p # digits 10

82061511001 p # digits 11

19072991752501 c # digits 14

8481116649425701 p # digits 16

17231203730201189308301 p # digits 23

Runtime (s) = 2.026000

Factorization of Six Lucas Sequence Numbers



Lucas 340

113709744839525149336680459091826532688903186653162057995534262332121127

# digits 72

7 p # digits 1

2161 p # digits 4

5441 p # digits 4

897601 p # digits 6

23230657239121 p # digits 14

17276792316211992881 p # digits 20

3834936832404134644974961 p # digits 25

Runtime (s) = 109.103000



Lucas 345

# digits 73

2 ^ 2 p # digits 1

31 p # digits 2

461 p # digits 3

1151 p # digits 4

1529 c # digits 4

324301 p # digits 6

686551 p # digits 6

1485571 p # digits 7

4641631 p # digits 7

19965899801 c # digits 11

117169733521 p # digits 12

3490125311294161 p # digits 16

Runtime (s) = 0.032000



Lucas 350

13985374084677485786380981408251904922622980674054858121032362563653278123

# digits 74

3 p # digits 1

401 p # digits 3

2801 p # digits 4

11521 c # digits 5

28001 p # digits 5

570601 p # digits 6

12317523121 p # digits 11

248773766357061401 p # digits 18

7358192362316341243805801 p # digits 25

Runtime (s) = 21.047000



Lucas 355

69362907070206748494476200566565775354902428015845969798000696945226974645

# digits 74

5 p # digits 1

4261 p # digits 4

6673 p # digits 4

75309701 p # digits 8

309273161 p # digits 9

46165371073 p # digits 11

9207609261398081 p # digits 16

49279722643391864192801 p # digits 23

Runtime (s) = 40.726000



Lucas 360

769246427201094785080787978422393713094534885688979999504447628313150135520

# digits 75

2 ^ 5 p # digits 1

3 ^ 2 p # digits 1

23 p # digits 2

41 p # digits 2

105 c # digits 3

107 p # digits 3

241 p # digits 3

2161 p # digits 4

2521 p # digits 4

3439 c # digits 4

8641 p # digits 4

20641 p # digits 5

103681 p # digits 6

109441 p # digits 6

191306797 c # digits 9

10783342081 p # digits 11

13373763765986881 p # digits 17

Runtime (s) = 0.032000



Lucas 365

19076060504701386559675231910437330047906343529583769121365013189782992678011

# digits 77

11 p # digits 2

151549 p # digits 6

514651 p # digits 6

7015301 p # digits 7

8942501 p # digits 7

9157663121 p # digits 10

11899937029 p # digits 11

3252336525249736694804553589211 p # digits 31

The following two numbers were first factorized by J. M. Pollard on an 8-bit Phillips P2012 personal computer with 64 KB RAM and two 640 KB disc drives. The times required by Pollard were 41 and 47 hours.

2^144-3
22300745198530623141535718272648361505980413

# digits 44

492729991333 p # digits 12

45259565260477899162010980272761 p # digits 32

Runtime (s) = 0.086000

2^153+3
11417981541647679048466287755595961091061972995

# digits 47

5 p # digits 1

11 p # digits 2

600696432006490087537 p # digits 21

345598297796034189382757 p # digits 24

Runtime (s) = 0.676000

Partial factorization of the Twelfth Fermat Number 2^4096+1
# digits 1234

114689 p # digits 6

26017793 p # digits 8

63766529 p # digits 8

190274191361 p # digits 12

Runtime (s) = 1532.878000