Category: C++ Computer Applications

Blog Entry © Thursday, March 19, 2026, by James Pate Williams, Jr. Solution of Laplace’s Two-Dimensional Equation on a Square Boundary

LaplaceFD2dDownload
#pragma once

class Laplace
{
private:
    static int N;

public:
    Laplace(int N);
    double** Solve(
        double L,
        double epsilon,
        int& iterations,
        double (*f)(double));
    void PrintSolution(
        double L,
        double epsilon,
        double** u,
        double** w,
        int& iterations);
};


#include "Laplace.h"
#include <float.h>
#include <math.h>
#include <iomanip>
#include <iostream>

int Laplace::N = 0;

static double fxy(double x, double y)
{
    return x * y * exp(-0.5 * (x * x + y * y));
}

static double gxy(double x, double y)
{
    return x * sin(y);
}

Laplace::Laplace(int N)
{
    this->N = N;
}

double** Laplace::Solve(
    double L, double epsilon, int& iterations,
    double (*f)(double))
{
    double h = L / (N + 1.0), uv = 0, max = 0;
    double** u = new double*[N + 2];

    if (u == nullptr)
        exit(-1);

    for (int i = 0; i < N + 2; i++)
    {
        u[i] = new double[N + 2];

        if (u[i] == nullptr)
            exit(-2);

        for (int j = 0; j < N + 2; j++)
            u[i][j] = 0.0;
    }

    double** v = new double* [N + 2];

    if (v == nullptr)
        exit(-3);

    for (int i = 0; i < N + 2; i++)
    {
        v[i] = new double[N + 2];

        if (v[i] == nullptr)
            exit(-4);

        for (int j = 0; j < N + 2; j++)
            v[i][j] = 0.0;
    }

    iterations = 0;

    for (int j = 0; j <= N + 1; j++)
        u[0][j] = 0.0;

    for (int j = 0; j <= N + 1; j++)
        u[N + 1][j] = 0.0;

    for (int i = 0; i <= N + 1; i++)
        u[i][0] = f(i * h);

    for (int i = 0; i <= N + 1; i++)
        u[i][N + 1] = 0.0;

    do
    {
        iterations++;

        for (int i = 1; i <= N; i++)
            for (int j = 1; j <= N; j++)
                u[i][j] = 0.25 * (u[i - 1][j] + u[i + 1][j] + u[i][j - 1] +
                    u[i][j + 1]);

        for (int i = 0; i <= N + 1; i++)
            for (int j = 0; j <= N + 1; j++)
                v[i][j] = u[i][j];

        for (int i = 1; i <= N; i++)
            for (int j = 1; j <= N; j++)
                v[i][j] = 0.25 * (v[i - 1][j] + v[i + 1][j] + v[i][j - 1] +
                    v[i][j + 1]);

        max = DBL_MIN;

        for (int i = 1; i <= N; i++)
        {
            for (int j = 1; j <= N; j++)
            {
                uv = fabs(v[i][j] - u[i][j]);
                max = fmax(uv, max);
            }
        }

        for (int i = 0; i <= N + 1; i++)
            for (int j = 0; j <= N + 1; j++)
                u[i][j] = v[i][j];

    } while (max > epsilon);

    for (int i = 0; i < N + 2; i++)
        delete v[i];

    return u;
}

void Laplace::PrintSolution(
    double L,
    double epsilon,
    double** u,
    double** w,
    int& iterations)
{
    double h = L / (N + 1.0);
    
    std::cout << std::endl;
    std::cout << "x\t\ty\t\tu(x, y)\t\tw(x, y)\t\tpercent error\r\n\r\n";

    for (int i = 0; i <= N + 1; i++)
    {
        double x = i * h;

        for (int j = 0; j <= N + 1; j++)
        {
            double y = j * h;
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << x << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << y << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << u[i][j] << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << w[i][j] << '\t';
            double error = 0.0;
            if (w[i][j] != 0.0)
                error = 100.0 * fabs(u[i][j] - w[i][j]) /
                    fabs(w[i][j]);
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << error << std::endl;
        }
    }

    std::cout << std::endl;
    std::cout << "iterations = " << iterations;
    std::cout << std::endl;
}


// LaplaceEquationFD2d.cpp (c) Tuesday, March 17, 2026 
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// by James Pate Williams, Jr., BA, ,BS, MSwE, PhD
// Reference: Numerical Computation 2 Methods,
// Software And Analysis (c) 1995 by Chrisoph E. Ueberhuber
// Pages 392 - 393 Poisson Matrices
// https://math.libretexts.org/Bookshelves/Differential_Equations/Introduction_to_Partial_Differential_Equations_(Herman)/06%3A_Problems_in_Higher_Dimensions/6.03%3A_Laplaces_Equation_in_2D

#include <math.h>
#include <iostream>
#include "Laplace.h"

double** w;

double SimpsonsRule(double L, double a, double b,
    int Nsr, int n, double (*f)(double, double, int))
{
    double h = (b - a) / Nsr;
    double h2 = 2.0 * h;
    double s = 0.0;
    double t = 0.0;
    double x = a + h;

    for (int i = 1; i < Nsr; i += 2)
    {
        s += f(L, x, n);
        x += h2;
    }

    x = a + h2;

    for (int i = 2; i < Nsr; i += 2)
    {
        t += f(L, x, n);
        x += h2;
    }

    return
        h * (f(L, a, n) + 4 * s + 2 * t + f(L, b, n)) / 3.0;
}

double f(double x)
{
    return x * exp(-x) * cos(x);
}

double g(double L, double x, int n)
{
    double pi = 4.0 * atan(1.0);
    double fx = f(x);
    return fx * sin(n * pi * x / L);
}

double bn(double L, int Nsr, int n)
{
    double cs = 2.0 / L;
    double sr = SimpsonsRule(L, 0.0, L, Nsr, n, g);
    return cs * sr;
}

void ComputeW(double L, int N)
{
    double h = L / (N + 1.0);

    w = new double* [N + 2];

    if (w == nullptr)
        exit(-5);

    for (int i = 0; i < N + 2; i++)
    {
        w[i] = new double[N + 2];

        if (w[i] == nullptr)
            exit(-6);

        for (int j = 0; j < N + 2; j++)
            w[i][j] = 0.0;
    }

    for (int j = 0; j <= N; j++)
        w[0][j] = 0.0;

    for (int j = 0; j <= N; j++)
        w[N + 1][j] = 0.0;

    for (int i = 0; i <= N; i++)
        w[i][0] = f(i * h);

    for (int i = 0; i <= N; i++)
        w[i][N + 1] = 0.0;

    double pi = 4.0 * atan(1.0);
    double x = h, y = 0.0;

    for (int i = 1; i <= N; i++)
    {
        y = h;

        for (int j = 1; j <= N; j++)
        {
            double s = 0.0;

            for (int n = 1; n <= 64; n++)
            {
                double npi = n * pi;
                double an = bn(L, 512, n) / sinh(npi);
                s += an * sin(npi * x / L) * sinh(npi * (L - y) / L);
            }

            w[i][j] = s;
            y += h;
        }

        x += h;
    }
}

int main()
{
    while (true)
    {
        char line[128] = "";
        std::cout << "N = ";
        std::cin.getline(line, 128);
        int N = atoi(line);

        if (N < 4 || N > 128)
        {
            std::cout << "N < 4 || N > 128, try again";
            std::cout << std::endl;
        }

        std::cout << "epsilon = ";
        std::cin.getline(line, 128);

        double epsilon = atof(line);
        Laplace l(N);

        double L = 1.0;
        int iterations = 0;
        ComputeW(L, N);
        double** u = l.Solve(L, epsilon, iterations, f);
        l.PrintSolution(L, epsilon, u, w, iterations);

        for (int i = 0; i < N + 2; i++)
        {
            if (u[i] != nullptr)
                delete u[i];

            if (w[i] != nullptr)
                delete w[i];
        }

        break;
    }

    return 0;
}

Blog Entry © Tuesday, March 17, 2026, by James Pate Williams, Jr., Comparison of Power Series Arctangent and Arcsine Functions with the C++ Built-In Functions

Blog Entry Tuesday March 17 2026Download
#pragma once

//https://scipp-legacy.pbsci.ucsc.edu/~haber/ph116A/taylor11.pdf

class Functions
{
public:
	static void initialize();
	static double factorial(int n);
	static double arccosine(double x, int n);
	static double arccosecant(double x);
	static double arccotangent(double x);
	static double arcsecant(double x);
	static double arcsine(double x, int n);
	static double arctangent(double x, int n);
};


#include "Functions.h"
#include <math.h> 

double pi2 = 0.0;

double Functions::factorial(int n)
{
	double nfactorial = 1.0;

	for (int i = 2; i <= n; i++)
		nfactorial *= i;

	return nfactorial;
}

void Functions::initialize()
{
	pi2 = 2.0 * atan(1.0);
}

double Functions::arccosine(double x, int n)
{
	double sum = pi2 - arcsine(x, n);
	return sum;
}

double Functions::arccosecant(double x)
{
	double sum = 0;
	return sum;
}

double Functions::arccotangent(double x)
{
	double sum = 0;
	return sum;
}

double Functions::arcsecant(double x)
{
	double sum = 0;
	return sum;
}

double Functions::arcsine(double x, int n)
{
	double sum = 0.0;

	if (fabs(x) <= 1.0)
	{
		for (int i = n; i >= 0; i--)
		{
			double ifact = factorial(i);
			double i2 = 2.0 * i, i21 = 2 * i + 1;
			double coeff = factorial(i2) /
				(pow(2, i2) * ifact * ifact * (i21));
			sum += coeff * pow(x, i21);
		}
	}

	return sum;
}

double Functions::arctangent(double x, int n)
{
	double sum = 0.0;

	if (fabs(x) <= 1.0)
	{
		double one = 0.0;

		if (n % 2 == 0)
			one = 1.0;
		else
			one = -1.0;

		for (int i = 2 * n + 1; i >= 0; i--)
		{
			one = -one;
			double i21 = 2 * i + 1.0;
			sum += one * pow(x, i21) / i21;
		}
	}

	return sum;
}


// InvTrigConsoleCPP.cpp (c) Monday, March 16, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD

#include "Functions.h"
#include <math.h>
#include <iomanip>
#include <iostream>

static void CreateTable(
	char title[],
	double a, double b, int n, int nPts,
	double(*fx)(double, int), double(*gx)(double))
{
	double x = a;
	double deltaX = (b - a) / nPts;
	std::cout << title << std::endl;
	std::cout << "# Terms = " << (2 * n + 2) << std::endl;
	std::cout << "x" << '\t' << "fx" << "\t\t";
	std::cout << "sx" << "\t\t" << "error" << std::endl;

	for (int i = 0; i < nPts; i++)
	{
		double f = fx(x, n);
		double s = gx(x);
		double e = 0.0;

		if (fabs(s) != 0)
			e = 100.0 * fabs(f - s) / fabs(s);

		std::cout << std::fixed << std::setprecision(4);
		std::cout << std::setw(4);
		std::cout << x << '\t';
		std::cout << std::fixed << std::setprecision(11);
		std::cout << std::setw(12);
		std::cout << f << '\t';
		std::cout << std::fixed << std::setprecision(11);
		std::cout << std::setw(12);
		std::cout << s << '\t';
		std::cout << e << std::endl;
		x += deltaX;
	}
}

int main()
{
	Functions::initialize();

	while (true)
	{
		char line[128] = "";
		std::cout << "== Menu == " << std::endl;
		std::cout << "1 arcsine" << std::endl;
		std::cout << "2 arctangent" << std::endl;
		std::cout << "3 exit" << std::endl;
		std::cout << "option # : ";
		std::cin.getline(line, 128);
		char option = line[0];

		if (option == '3')
			break;

		if (option < '1' || option > '2')
		{
			std::cout << "invalid option" << std::endl;
			continue;
		}

		std::cout << "# points: ";
		std::cin.getline(line, 128);
		int nPts = atoi(line);

		std::cout << "# terms: ";
		std::cin.getline(line, 128);
		int n = atoi(line);

		if (option == '1')
		{
			char title[] = "Series arcsin Versus C++ asin";
			CreateTable(title, 0.0, 1.0, n, nPts,
				Functions::arcsine, asin);
		}
		else if (option == '2')
		{
			char title[] = "Series arctan Versus C++ atan";
			CreateTable(title, 0.0, 1.0, n, nPts,
				Functions::arctangent, atan);
		}
	}

	return 0;
}

Blog Entry © Sunday March 15, 2026, by James Pate Williams, Jr. Graphing Integer Order Bessel Functions of the First Kind and the real Gamma Function

Blog Entry Sunday March 15 2026Download

Blog Entry © Saturday, March 14, 2026, by James Pate Williams, Jr. Power Series Solution of a Second Order Linear Ordinary Differential Equation

Blog Entry Saturday March 14 2026Download
// SecondOrderLinearODE.cpp (c) Friday March 13, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Solves the Bessel Equation of the First Kind for Order 0
// Using the power series found in the reference below - 
// References: https://math.libretexts.org/Bookshelves/Calculus/Calculus_(OpenStax)/17%3A_Second-Order_Differential_Equations/17.04%3A_Series_Solutions_of_Differential_Equations
// https://mathworld.wolfram.com/BesselDifferentialEquation.html

#include "framework.h"
#include <Windows.h>
#include "SecondOrderLinearODE.h"
#include <string>
#include <vector>

#define MAX_LOADSTRING 100

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
WCHAR fTitle[128];
double chi[1025], eta[1025], coeff[65];

// 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_SECONDORDERLINEARODE, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    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_SECONDORDERLINEARODE));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_SECONDORDERLINEARODE);
    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;
}

static double Factorial(int n)
{
    double factorial = 1.0;

    for (int i = 2; i <= n; i++)
        factorial *= i;

    return factorial;
}

static void ComputeCoefficients()
{
    double a0 = 1.0;

    for (int k = 0; k <= 64; k++)
    {
        double kFactorial = Factorial(k);

        coeff[k] = a0 * pow(-1.0, k) /
            (pow(2.0, 2.0 * k) * pow(kFactorial, 2.0));
    }
}

static void ComputePoints()
{
    double chi0 = -15.0;
    double deltaChi = 30.0 / 1024;
 
    for (int i = 0; i <= 1024; i++)
    {
        double sum = 0.0;

        for (int k = 64; k >= 0; k--)
            sum += pow(chi0, 2 * k) * coeff[k];

        chi[i] = chi0;
        eta[i] = sum;

        chi0 += deltaChi;
    }
}

static void FindMinMax(
    double x[],
    double& min,
    double& max)
{
    min = x[0];
    max = x[0];

    for (int i = 1; i <= 1024; i++)
    {
        if (x[i] < min)
            min = x[i];
        if (x[i] > max)
            max = x[i];
    }
}

static void DrawTitles(HDC hdc, RECT clientRect, const std::wstring& fTitle,
    int sx0, int sx1, int sy0, int sy1)
{
    HFONT hCustomFont = CreateFont(
        -MulDiv(10, GetDeviceCaps(hdc, LOGPIXELSY), 72), // Height in logical units
        0,                      // Width (0 = default)
        0,                      // Escapement
        0,                      // Orientation
        FW_BOLD,                // Weight (700 = bold)
        FALSE,                  // Italic
        FALSE,                  // Underline
        FALSE,                  // StrikeOut
        DEFAULT_CHARSET,        // Charset
        OUT_DEFAULT_PRECIS,     // Output precision
        CLIP_DEFAULT_PRECIS,    // Clipping precision
        DEFAULT_QUALITY,        // Quality
        FIXED_PITCH | FF_MODERN,// Pitch and family
        L"Courier New"          // Typeface name
    );

    HFONT hOldFont = (HFONT)SelectObject(hdc, hCustomFont);

    SIZE sz;
    int width = clientRect.right - clientRect.left;

    // Title: fTitle
    GetTextExtentPoint32(hdc, fTitle.c_str(), (int)fTitle.length(), &sz);
    int w = sz.cx;
    int h = sz.cy;
    TextOut(hdc, (width - w) / 2, h, fTitle.c_str(), (int)fTitle.length());

    // x-axis title: "x"
    std::wstring xTitle = L"x";
    GetTextExtentPoint32W(hdc, xTitle.c_str(), (int)xTitle.length(), &sz);
    w = sz.cx;
    TextOutW(hdc, sx0 + (sx1 - sx0 - w) / 2, sy1 + 2 * h, xTitle.c_str(), (int)xTitle.length());

    // y-axis title: "p(x)"
    std::wstring yTitle = L"J0(x)";
    GetTextExtentPoint32W(hdc, yTitle.c_str(), (int)yTitle.length(), &sz);
    w = sz.cx;
    TextOutW(hdc, sx1 + w / 5, sy0 + (sy1 - sy0) / 2 - h / 2, yTitle.c_str(), (int)yTitle.length());

    SelectObject(hdc, hOldFont); // Restore original font
}

static void DrawFormattedText(HDC hdc, char loctext[], RECT rect)
{
    // Draw the text with formatting options
    DrawTextA(hdc, loctext, (int)strlen(loctext), &rect, DT_SINGLELINE | DT_NOCLIP);
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    switch (message)
    {
    case WM_CREATE:
        ComputeCoefficients();
        ComputePoints();
        break;
    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_EXIT:
                DestroyWindow(hWnd);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
    {
        PAINTSTRUCT ps;
        HDC hdc = BeginPaint(hWnd, &ps);

        ComputePoints();
        double xMin = 0, yMin = 0;
        double xMax = 0, yMax = 0;
        FindMinMax(chi, xMin, xMax);
        FindMinMax(eta, yMin, yMax);

        double h = 0, pi = 0, plm = 0, theta = 0;
        float xSpan = (float)(xMax - xMin);
        float ySpan = (float)(yMax - yMin);
        RECT rect = { };
        GetClientRect(hWnd, &rect);
        float width = (float)(rect.right - rect.left + 1);
        float height = (float)(rect.bottom - rect.top - 32 + 1);
        float sx0 = 2.0f * width / 16.0f;
        float sx1 = 14.0f * width / 16.0f;
        float sy0 = 2.0f * height / 16.0f;
        float sy1 = 14.0f * height / 16.0f;
        float deltaX = xSpan / 8.0f;
        float deltaY = ySpan / 8.0f;
        float xSlope = (sx1 - sx0) / xSpan;
        float xInter = (float)(sx0 - xSlope * xMin);
        float ySlope = (sy0 - sy1) / ySpan;
        float yInter = (float)(sy0 - ySlope * yMax);
        float px = 0, py = 0, sx = 0, sy = 0;
        POINT wPt = { 0 };
        int i = 0;
        float x = (float)xMin;
        float y = (float)yMax;
        wcscpy_s(fTitle, 128, L"Bessel Function of the First Kind Order 0");
        DrawTitles(hdc, rect, fTitle, (int)sx0, (int)sx1, (int)sy0, (int)sy1);


        px = (float)chi[0];
        py = (float)eta[0];

        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        MoveToEx(hdc, (int)sx, (int)sy0, &wPt);

        while (i <= 8)
        {
            sx = xSlope * x + xInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
            LineTo(hdc, (int)sx, (int)sy1);
            char numberStr[128] = "";
            sprintf_s(numberStr, 128, "%5.4lf", x);
            SIZE size = { };
            GetTextExtentPoint32A(
                hdc,
                numberStr,
                (int)strlen(numberStr),
                &size);
            RECT textRect = { 0 };
            textRect.left = (long)(sx - size.cx / 2.0f);
            textRect.right = (long)(sx + size.cx / 2.0f);
            textRect.top = (long)sy1;
            textRect.bottom = (long)(sy1 + size.cy / 2.0f);
            DrawFormattedText(hdc, numberStr, textRect);
            x += deltaX;
            i++;
        }

        i = 0;
        y = (float)yMin;

        while (i <= 8)
        {
            sy = ySlope * y + yInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
            LineTo(hdc, (int)sx, (int)sy);

            if (i != 0)
            {
                char numberStr[128] = "";
                sprintf_s(numberStr, 128, "%+5.3lf", y);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    numberStr,
                    (int)strlen(numberStr),
                    &size);
                RECT textRect = { 0 };
                textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
                textRect.right = (long)(sx0 - size.cx / 2.0f);
                textRect.top = (long)(sy - size.cy / 2.0f);
                textRect.bottom = (long)(sy + size.cy / 2.0f);
                DrawFormattedText(hdc, numberStr, textRect);
            }

            y += deltaY;
            i++;
        }

        HGDIOBJ nPenNew = NULL;
        HGDIOBJ hPenOld = NULL;

        nPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
        hPenOld = SelectObject(hdc, nPenNew);

        HRGN clipRegion = CreateRectRgn(
            (int)sx0, (int)sy0,              // Left, Top
            (int)(sx1), (int)(sy1)           // Right, Bottom
        );

        // Apply clipping region

        SelectClipRgn(hdc, clipRegion);


        px = (float)chi[0];
        py = (float)eta[0];

        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        wPt.x = wPt.y = 0;
        MoveToEx(hdc, (int)sx, (int)sy, &wPt);

        for (size_t j = 1; j <= 1024; j++)
        {

            px = (float)chi[j];
            py = (float)eta[j];

            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            LineTo(hdc, (int)sx, (int)sy);
        }

        DeleteObject(nPenNew);
        nPenNew = NULL;

        SelectObject(hdc, hPenOld);
        DeleteObject(clipRegion);

        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 © Thursday and Friday, March 12 – 13, 2026, Naive Power Series Computations of the Cosine and Sine Trigonometric Functions and Polynomial Fit to the Trig Functions Using Least Squares

Blog Entry Thursday March 12 2026Download
// TrigLSPolyFit.cpp (c) Thursday, March 12, 2026
// by James Pate Williams, Jr.
// Least Squares Curve Fitting to the Trigonometric
// functions Cosine and Sine

#include <iomanip>
#include <iostream>
#include <vector>
#include "PolyLSFit.h"

static const int HornerN = 50;
double even_factors[HornerN], odd_factors[HornerN + 1];

static double factorial(int n)
{
	double nfactorial = 1.0;

	for (int i = 2; i <= n; i++)
		nfactorial *= i;
	return nfactorial;
}

static void computeEvenFactors()
{
	int one = 1;

	for (int n = 0; n < 50; n += 2)
	{
		even_factors[n] = one / factorial(n);
		one = -one;
	}
}

static void computeOddFactors()
{
	int one = 1;

	for (int n = 1; n < 51; n += 2)
	{
		odd_factors[n] = one / factorial(n);
		one = -one;
	}
}

static double fx_cos(double x, int n)
{
	double sum = even_factors[n];

	for (int i = n - 1; i >= 0; i--)
	{
		sum = sum * x + even_factors[i];
	}

	return sum;
}

static double fx_sin(double x, int n)
{
	double sum = odd_factors[n];

	for (int i = n - 1; i >= 0; i--)
	{
		sum = sum * x + odd_factors[i];
	}

	return sum;
}

static const int M = 33, N = 100, N21 = N + N + 1;
std::vector<double> cosX(N21), cosY(N21), cosCoeff;
std::vector<double> sinX(N21), sinY(N21), sinCoeff;
std::vector<double> pX(M), pY(M);

static void LSPolyFit(int HornerN, int degree)
{
	double pi = 4.0 * atan(1.0);
	double deltaX = 2.0 * pi / N;
	double X = -pi;

	for (int i = 0; i < N21; i++)
	{
		cosX[i] = X;
		cosY[i] = fx_cos(X, HornerN);
		sinX[i] = X;
		sinY[i] = fx_sin(X, HornerN + 1);
		X += deltaX;
	}

	PolyLSFit::leastSquares(cosX, cosY, cosCoeff, degree, N21);
	PolyLSFit::leastSquares(sinX, sinY, sinCoeff, degree, N21);

	std::cout << std::setw(6) << "x";
	std::cout << std::setw(17) << "fit cos";
	std::cout << std::setw(17) << "C++ cos";
	std::cout << std::setw(17) << "% Error";
	std::cout << std::endl;

	deltaX = 2.0 * pi / (M - 1);
	X = -pi;

	for (int i = 0; i < M; i++)
	{
		double Y = cosCoeff[degree];

		for (int j = degree - 1; j >= 0; j--)
		{
			Y = Y * X + cosCoeff[j];
		}

		pX[i] = X;
		pY[i] = Y;

		std::cout << std::fixed << std::setprecision(3);
		std::cout << std::setw(6);
		std::cout << X << ' ';

		if (Y >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);
		
		std::cout << fabs(Y) << ' ';
		double cx = cos(X);

		if (cx >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);
		
		std::cout << fabs(cx) << ' ';
		double pcDifference = 0.0;

		if (fabs(cx) >= 1.0e-12)
			pcDifference = 100.0 * fabs(Y - cx) / fabs(cx);
		else
			pcDifference = -1.0;

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		if (pcDifference != -1.0)
			std::cout << pcDifference << std::endl;
		else
			std::cout << "not a number" << std::endl;

		X += deltaX;
	}

	std::cout << std::endl;

	std::cout << std::setw(6) << "x";
	std::cout << std::setw(17) << "fit sin";
	std::cout << std::setw(17) << "C++ sin";
	std::cout << std::setw(17) << "% Error";
	std::cout << std::endl;

	deltaX = 2.0 * pi / (M - 1);
	X = -pi;

	for (int i = 0; i < M; i++)
	{
		double Y = sinCoeff[degree];

		for (int j = degree - 1; j >= 0; j--)
		{
			Y = Y * X + sinCoeff[j];
		}

		pX[i] = X;
		pY[i] = Y;

		std::cout << std::fixed << std::setprecision(3);
		std::cout << std::setw(6);
		std::cout << X << ' ';

		if (Y >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		std::cout << fabs(Y) << ' ';
		double sx = sin(X);

		if (sx >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		std::cout << fabs(sx) << ' ';
		double pcDifference = 0.0;

		if (fabs(sx) >= 1.0e-12)
			pcDifference = 100.0 * fabs(Y - sx) / fabs(sx);
		else
			pcDifference = -1.0;

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		if (pcDifference != -1.0)
			std::cout << pcDifference << std::endl;
		else
			std::cout << "not a number" << std::endl;

		X += deltaX;
	}

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

int main()
{
	computeEvenFactors();
	computeOddFactors();
	std::cout << "# terms 16" << std::endl;
	LSPolyFit(HornerN, 16);
	std::cout << "# terms 20" << std::endl;
	LSPolyFit(HornerN, 20);
	std::cout << "# terms 24" << std::endl;
	LSPolyFit(HornerN, 24);
	return 0;
}

Blog Entry © Thursday, February 5, 2026, by James Pate Williams, Jr., Three Toom-Cook 3-Way Multiplication Implementations

I translated some of Marco Bodrato ‘s pseudocode found on his website:

Optimal Toom-Cook Polynomial Multiplication, by Marco Bodrato

To Microsoft Win32 C++. I used Microsoft’s Visual Studio 2022 Community Version.

Balanced Toom-Cook 3-Way Multiplication

U = 7

V = 8

W = 56

U = 9

V = 9

W = 81

U = 18

V = 17

W = 306

U = 123

V = 456

W = 56088

U = 1234

V = 5678

W = 7006652

U = 12345

V = 67890

W = 838102050

U = 123456

V = 789012

W = 97408265472

D:\ToomCookBodrato\x64\Debug\ToomCookBodrato.exe (process 30392) exited with code 0 (0x0).

Unbalanced Toom-Cook 3-Way Multiplication

U = 20

V = 5

W = 100

U = 112

V = 20

W = 2240

U = 1234

V = 567

W = 699678

U = 12345

V = 6789

W = 83810205

U = 123456

V = 78901

W = 9740801856

U = 1234567

V = 890123

W = 1098916481741

U = 12345678

V = 9012345

W = 111263509394910

D:\ToomCookBodrato\x64\Debug\ToomCookBodrato.exe (process 13108) exited with code 0 (0x0).

Asymmetrical Squaring, Splitting in Five

U = 8

W = 64

U = 12

W = 144

U = 321

W = 103041

U = 1234

W = 1522756

U = 54321

W = 2950771041

U = 123456

W = 15241383936

D:\ToomCookBodrato\x64\Debug\ToomCookBodrato.exe (process 27592) exited with code 0 (0x0).

Press any key to close this window . . .

The accompanying source code is released under the GNU General Public License, version 3. The full license text is included in the downloadable archive and is also available at:
https://www.gnu.org/licenses/gpl-3.0.html

/*
   Toom–Cook Multiplication (Bodrato-inspired implementation)
   Copyright (C) 2026 James Pate Williams, Jr.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#pragma once

// Algorithm structure inspired by Marco Bodrato’s Toom–Cook notes.
// See: http://www.bodrato.it/toom-cook/
// All implementation code and comments in this file are by
// James Pate Williams, Jr., BA, BS, Master of Software Engineering,
// Doctor of Philosophy in Computer Science

class Multiplication
{

public:

	static long long BalancedToomCook3Way(
		long long b, long long U, long long V);
	static long long UnbalancedToomCook3Way(
		long long b, long long U, long long V);
	static long long AsymeticalSquaring(
		long long b, long long U);
};


/*
   Toom–Cook Multiplication (Bodrato-inspired implementation)
   Copyright (C) 2026 James Pate Williams, Jr.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#include "pch.h"
#include "Multiplication.h"

static void BaseBDigits(
	long long B, long long u,
	std::vector<long long>& digits)
{
	while (u > 0)
	{
		long long d = u % B;
		u = u / B;
		digits.push_back(d);
	}
}

long long Multiplication::BalancedToomCook3Way(
	long long b, long long U, long long V)
{
	double k = 3.0;
	double logU = log(U) / log(b);
	double logV = log(V) / log(b);
	double floorU = floor(floor(logU) / k);
	double floorV = floor(floor(logV) / k);
	double i = fmax(floorU, floorV) + 1.0;
	long long B = static_cast<long long>(pow(b, i));
	long long W = 0;
	std::vector<long long> UDigits;
	std::vector<long long> VDigits;

	for (long long j = 0; j < 5; j++)
	{
		long long x = 0, y = 0;

		if (j == 0)
			x = 1;
		else if (j == 1)
			x = -2;
		else if (j == 2)
			x = 1;
		else if (j == 3)
			x = -1;
		else
			x = 1;

		if (j == 0)
			y = 1;
		else if (j == 1)
			y = 1;
		else if (j == 2)
			y = 1;
		else if (j == 3)
			y = 1;
		else
			y = 0;

		UDigits.clear();
		VDigits.clear();

		BaseBDigits(B, U, UDigits);
		BaseBDigits(B, V, VDigits);

		long long U0 = 0, U1 = 0, U2 = 0;
		long long V0 = 0, V1 = 0, V2 = 0;

		if (UDigits.size() == 3)
		{
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (VDigits.size() == 3)
		{
			V2 = VDigits[2];
			V1 = VDigits[1];
			V0 = VDigits[0];
		}

		if (UDigits.size() == 2)
		{
			U2 = 0;
			U1 = UDigits[1];
			U2 = UDigits[0];
		}

		if (VDigits.size() == 2)
		{
			V2 = 0;
			V1 = VDigits[1];
			V0 = VDigits[0];
		}

		if (UDigits.size() == 1)
		{
			U2 = 0;
			U1 = 0;
			U0 = UDigits[0];
		}

		if (VDigits.size() == 1)
		{
			V2 = 0;
			V1 = 0;
			V0 = VDigits[0];
		}

		if (UDigits.size() == 0 ||
			VDigits.size() == 0)
			break;

		long long x2 = x * x;
		long long y2 = y * y;
		long long xy = x * y;
		
		U = U2 * x2 + U1 * xy + U0 * y2;
		V = V2 * x2 + V1 * xy + V0 * y2;
		
		if (U == 0 || V == 0)
			break;

		long long W3 = U0 + U2;
		long long W2 = V0 + V2;
		long long W0 = W3 - U1;
		long long W4 = W2 - V1;
		
		W3 = W3 + U1;
		W2 = W2 + V1;

		long long W1 = W3 * W2;

		W2 = W0 * W4;
		W0 = ((W0 + U2) << 1) - U0;
		W4 = ((W4 + V2) << 1) - V0;
		W3 = W0 * W4;
		W0 = U0 * V0;
		W4 = U2 * V2;

		W3 = (W3 - W1) / k;
		W1 = (W1 - W2) >> 1;
		W2 = W2 - W0;
		W3 = ((W2 - W3) >> 1) + (W4 << 1);
		W2 = W2 + W1;

		W3 = W4 * x + W3 * y;
		W1 = W2 * x + W1 * y;
		W1 = W1 - W3;

		W = W3 * x * x2 + W1 * x * y2 + W0 * y2 * y2;
	}

	return W;
}

long long Multiplication::UnbalancedToomCook3Way(
	long long b, long long U, long long V)
{
	double k = 3.0;
	double logU = log(U) / log(b);
	double logV = log(V) / log(b);
	double floorU = floor(floor(logU) / k);
	double floorV = floor(floor(logV) / k);
	double i = fmax(floorU, floorV) + 1.0;
	long long B = static_cast<long long>(pow(b, i));
	long long W = 0;
	std::vector<long long> UDigits;
	std::vector<long long> VDigits;

	for (long long j = 0; j < 5; j++)
	{
		long long x = 0, y = 0;

		if (j == 0)
			x = 1;
		else if (j == 1)
			x = -2;
		else if (j == 2)
			x = 1;
		else if (j == 3)
			x = -1;
		else
			x = 1;

		if (j == 0)
			y = 1;
		else if (j == 1)
			y = 1;
		else if (j == 2)
			y = 1;
		else if (j == 3)
			y = 1;
		else
			y = 0;

		long long U0 = 0, U1 = 0, U2 = 0;
		long long V0 = 0, V1 = 0, V2 = 0;
		long long U3 = 0;
		long long x2 = x * x, y2 = y * y;
		long long x3 = x * x2, y3 = y * y2;

		UDigits.clear();
		VDigits.clear();

		BaseBDigits(B, U, UDigits);
		BaseBDigits(B, V, VDigits);

		if (UDigits.size() == 4)
		{
			U3 = UDigits[3];
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 3)
		{
			U3 = 0;
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (VDigits.size() == 3)
		{
			V2 = VDigits[2];
			V1 = VDigits[1];
			V0 = VDigits[0];
		}

		if (UDigits.size() == 2)
		{
			U3 = 0;
			U2 = 0;
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (VDigits.size() == 2)
		{
			V2 = 0;
			V1 = VDigits[1];
			V0 = VDigits[0];
		}

		if (UDigits.size() == 1)
		{
			U3 = 0;
			U2 = 0;
			U1 = 0;
			U0 = UDigits[0];
		}

		if (VDigits.size() == 1)
		{
			V2 = 0;
			V1 = 0;
			V0 = VDigits[0];
		}

		if (UDigits.size() == 0 ||
			VDigits.size() == 0)
			break;

		U = U3 * x3 + U2 * x2 * y + U1 * x * y2 + U0 * y3;
		V = V1 * x + V0 * y;

		if (U == 0 || V == 0)
			break;

		long long W3 = U0 + U2;
		long long W2 = V0 + V1;
		long long W1 = U1 + U3;
		long long W4 = V0 - V1;
		long long W0 = W3 - W1;
		
		W3 = W3 + W1;
		W1 = W3 * W2;
		W2 = W0 * W4;
		W0 = U0 - (((U1 - ((U2 - U3) << 1)) << 1) << 1);
		W4 = W4 - V1;
		W3 = W0 * W4;
		W0 = U0 * V0;
		W4 = U3 * V1;

		W3 = (W3 - W1) / k;
		W1 = (W1 - W2) >> 1;
		W2 = W2 - W0;
		W3 = ((W2 - W3) >> 1) + (W4 << 1);
		W2 = W2 + W1;

		W3 = W4 * x + W3 * y;
		W1 = W2 * x + W1 * y;
		W1 = W1 - W3;

		W = W3 * x3 + W1 * x * y2 + W0 * y2 * y2;
	}

	return W;
}

long long Multiplication::AsymeticalSquaring(
	long long b, long long U)
{
	double k = 5.0;
	double logU = log(U) / log(b);
	double floorU = floor(floor(logU) / k);
	double i = floorU + 1.0;
	long long B = static_cast<long long>(pow(b, i));
	long long W = 0, Wp = 0;
	std::vector<long long> UDigits;

	for (long long j = 0; j < 9; j++)
	{
		long long x = 0;

		if (j == 0)
			x = 0;
		else if (j == 1)
			x = 0;
		else if (j == 2)
			x = -1;
		else if (j == 3)
			x = 1;
		else if (j == 4 || j == 5 || j == 6)
			continue;
		else if (j == 7 || j == 8)
			x = 1;

		UDigits.clear();
		BaseBDigits(B, U, UDigits);
		
		long long x2 = x * x, x3 = x * x2, x4 = x2 * x2;
		long long U0 = 0, U1 = 0, U2 = 0, U3 = 0, U4 = 0;

		if (UDigits.size() == 5)
		{
			U4 = UDigits[4];
			U3 = UDigits[3];
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 4)
		{
			U4 = 0;
			U3 = UDigits[3];
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 3)
		{
			U4 = 0;
			U3 = 0;
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 2)
		{
			U4 = 0;
			U3 = 0;
			U2 = 0;
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 1)
		{
			U4 = 0;
			U3 = 0;
			U2 = 0;
			U1 = 0;
			U0 = UDigits[0];
		}

		U = U4 * x4 + U3 * x3 + U2 * x2 + U1 * x + U0;

		long long W0 = U0 - U3;
		long long W1 = U3 - U1;
		long long W6 = U1 - U2;
		long long W4 = U1 + U2;
		long long W5 = W6 - U4;
		long long W3 = W5 + (W0 << 1);
		W0 = W0 - W5;
		W6 = W0 + (W6 << 1);
		long long W7 = W6 + W1;
		W5 = W7 + W1;
		long long W8 = W5 + (W4 << 1);
		W4 = W4 - U4;

		long long W2 = (W4) * (W3);
		W4 = (W6) * (W5);
		W3 = (W7) * (W1);
		W1 = U0 * U1 * 2;
		W7 = U3 * U4 * 2;
		W5 = (W8) * (W8); W6 = (W0) * (W0);
		W0 = U0 * U0;
		W8 = U4 * U4;

		W2 = (W4) * (W3);
		W4 = (W6) * (W5);
		W3 = (W7) * (W1);
		W1 = U0 * U1 * 2;
		W7 = U3 * U4 * 2;
		W5 = (W8) * (W8);
		W6 = (W0) * (W0);
		W0 = U0 * U0;
		W8 = U4 * U4;

		W6 = (W6 + W5) >> 1;
		W5 = W5 - W6;
		W4 = (W4 + W6) >> 1;
		W3 = W3 + (W5 >> 1);
		W6 = W6 - W4;
		W5 = W5 - W3 - W1;
		W4 = W4 - W8 - W0;
		W3 = W3 - W7;
		W2 = W2 - W8 - W1 - W7 + W4 + W5;
		W6 = W6 - W2;

		W = W8 * x4 * x4 + W7 * x3 * x4 + W6 * x3 * x3 +
			W5 * x * x4 + W4 * x4 + W3 * x3 + W2 * x2 + W1 * x + W0;
		
		if (W != 0)
			Wp = W;

		if (Wp != 0)
			break;
	}

	return W;
}


/*
   Toom–Cook Multiplication (Bodrato-inspired implementation)
   Copyright (C) 2026 James Pate Williams, Jr.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#include "pch.h"
#include <iostream>
#include "Multiplication.h"

int main()
{
    bool balanced = false, square = true;

    if (balanced && !square)
    {
        long long U = 7;
        long long V = 8;
        long long W = 0;
        // W = U * V = 56
        W = Multiplication::BalancedToomCook3Way(10, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 9;
        V = 9;
        W = Multiplication::BalancedToomCook3Way(10, U, V);
        // W = U * V = 81
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 18;
        V = 17;
        W = Multiplication::BalancedToomCook3Way(100, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 306
        U = 123;
        V = 456;
        W = Multiplication::BalancedToomCook3Way(1000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 56,088
        U = 1234;
        V = 5678;
        W = Multiplication::BalancedToomCook3Way(10000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 7,006,652
        U = 12345;
        V = 67890;
        W = Multiplication::BalancedToomCook3Way(100000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 838,102,050
        U = 123456;
        V = 789012;
        W = Multiplication::BalancedToomCook3Way(1000000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 838,102,050
    }

    else if (!balanced && !square)
    {
        long long U = 20;
        long long V = 5;
        long long W;
        // W = U * V = 100
        W = Multiplication::UnbalancedToomCook3Way(100, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 112;
        V = 20;
        W = Multiplication::UnbalancedToomCook3Way(1000, U, V);
        // W = U * V = 2240
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 1234;
        V = 567;
        W = Multiplication::UnbalancedToomCook3Way(10000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 699,678
        U = 12345;
        V = 6789;
        W = Multiplication::UnbalancedToomCook3Way(100000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 83,810,205
        U = 123456;
        V = 78901;
        W = Multiplication::UnbalancedToomCook3Way(1000000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 97,408,265,472
        U = 1234567;
        V = 890123;
        W = Multiplication::UnbalancedToomCook3Way(10000000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 9,740,885,481,741
        U = 12345678;
        V = 9012345;
        W = Multiplication::UnbalancedToomCook3Way(100000000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 111,263,509,394,910
    }

    else if (square)
    {
        long long U = 8;
        long long W = Multiplication::AsymeticalSquaring(10, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 12;
        W = Multiplication::AsymeticalSquaring(100, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 321;
        W = Multiplication::AsymeticalSquaring(1000, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 1234;
        W = Multiplication::AsymeticalSquaring(10000, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 54321;
        W = Multiplication::AsymeticalSquaring(100000, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 123456;
        W = Multiplication::AsymeticalSquaring(1000000, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
    }
    return 0;
}

Blog Entry © Wednesday, February 4, 2026, by James Pate Williams, Jr. and the Microsoft Copilot Two Multiplication Algorithms from Wikipedia

[BEGIN COPILOT TEXT]

## Introduction

This small C++ console program implements two classical multiplication algorithms: 

**long multiplication** (the grade‑school method) and **Karatsuba multiplication**, the first sub‑quadratic multiplication algorithm discovered (Anatoly Karatsuba, 1960).

My goal here is not to optimize or modernize the algorithms, but to preserve their structure in a clear, readable form. The implementations follow the versions presented in the Wikipedia “Multiplication” article as faithfully as possible.

The code is intentionally simple and procedural. It uses 1‑based indexing for digit arrays because that mirrors the mathematical notation more closely than 0‑based indexing.

## Long Multiplication

This function implements the classical long multiplication algorithm in an arbitrary base. 

Digits are stored in reverse order (least significant digit first), which simplifies carry propagation.

The implementation below is intentionally direct and mirrors the textbook algorithm step‑by‑step

## Karatsuba Multiplication (Base 10)

This is a minimal, non‑recursive implementation of the Karatsuba method. 

The inputs are split into high and low parts:

x = x_1 * B ^ m + x_0

y = y_1 * B ^ m + y_0

Karatsuba reduces the number of multiplications from 4 to 3 by computing:

z = x * y = z_2 *B ^ 2m + z_1 * B ^ m + z_0

This implementation follows the Wikipedia pseudocode closely.

[END COPILOT TEXT]

#pragma once

class MyAlgorithms
{

public:

	static void LongMultiplication(
		int base, int p, int q, int& pSize,
		std::vector<int> a,
		std::vector<int> b,
		std::vector<int>& product
	);

	static void KaratsubaBase10(
		int x0, int x1, int y0, int y1,
		int B, int m, long long& z);
};

#include "pch.h"
#include "MyAlgorithms.h"

void MyAlgorithms::LongMultiplication(
	int base, int p, int q, int& pSize,
	std::vector<int> a, std::vector<int> b,
	std::vector<int>& product)
{
	pSize = p + q;
	product.resize(pSize + 1LL);

	for (int b_i = 1; b_i <= q; b_i++)
	{
		int carry = 0;

		for (int a_i = 1; a_i <= p; a_i++)
		{
			product[static_cast<long long>(a_i) + b_i - 1LL] +=
				carry + a[a_i] * b[b_i];
			carry = product[static_cast<long long>(a_i) + b_i - 1] / base;
			product[static_cast<long long>(a_i) + b_i - 1] =
				product[static_cast<long long>(a_i) + b_i - 1] % base;
		}

		product[static_cast<long long>(b_i) + p] = carry;
	}
}

void MyAlgorithms::KaratsubaBase10(
	int B, int m, int x0, int x1,
	int y0, int y1, long long& z)
{
	int pb = static_cast<int>(pow(B, m));
	int x = x1 * pb + x0;
	int y = y1 * pb + y0;
	int z2 = x1 * y1;
	int z1 = x1 * y0 + x0 * y1;
	int z0 = x0 * y0;

	z = z2 * static_cast<int>(pow(B, 2 * m)) + z1 * pb + z0;
}


#include "MyAlgorithms.h"

static void DoLongMultiplication()
{
	int base = 0, p = 0, q = 0, pSize = 0;
	char line[128] = "";
	char inputaStr[128] = "";
	char inputbStr[128] = "";
	char* aReverseStr = nullptr;
	char* bReverseStr = nullptr;
	std::cout << "Enter base = ";
	std::cin.getline(line, 128);
	base = atoi(line);
	std::cout << "a = ";
	std::cin.getline(inputaStr, 128);
	std::cout << "b = ";
	std::cin.getline(inputbStr, 128);
	aReverseStr = _strrev(inputaStr);
	bReverseStr = _strrev(inputbStr);
	p = static_cast<int>(strlen(aReverseStr));
	q = static_cast<int>(strlen(bReverseStr));
	pSize = p + q;
	std::vector<int> a(p + 1);
	std::vector<int> b(q + 1);
	std::vector<int> ab(p + q + 1);
	std::vector<int> product;

	for (int i = 1; i <= p; i++)
		a[i] = aReverseStr[i - 1] - '0';

	for (int i = 1; i <= q; i++)
		b[i] = bReverseStr[i - 1] - '0';

	MyAlgorithms::LongMultiplication(
		base, p, q, pSize, a, b, ab);

	size_t i = ab.size() - 1, j = 1;

	while (i >= 0)
	{
		if (ab[i] == 0)
			i--;
		else
			break;
	}

	product.push_back(0);

	for (j = i; j >= 1; j--)
		product.push_back(ab[j]);

	std::cout << "product = ";

	for (int i = 1; i < product.size(); i++)
		std::cout << product[i];

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

static void DoKaratsuba()
{
	char line[128] = "";
	std::cout << "Enter base = ";
	std::cin.getline(line, 128);
	int B = atoi(line);
	std::cout << "Enter m = ";
	std::cin.getline(line, 128);
	int m = atoi(line);
	std::cout << "x1 = ";
	std::cin.getline(line, 128);
	int x1 = atoi(line);
	std::cout << "x0 = ";
	std::cin.getline(line, 128);
	int x0 = atoi(line);
	std::cout << "y1 = ";
	std::cin.getline(line, 128);
	int y1 = atoi(line);
	std::cout << "y0 = ";
	std::cin.getline(line, 128);
	int y0 = atoi(line);
	long long z = 0;
	MyAlgorithms::KaratsubaBase10(
		B, m, x0, x1, y0, y1, z);
	std::cout << "z = " << z << std::endl;
}

int main()
{
	while (true)
	{
		char line[128] = "";
		std::cout << "== Menu ==" << std::endl;
		std::cout << "1 Long Multiplication" << std::endl;
		std::cout << "2 Karatsuba Multiplication" << std::endl;
		std::cout << "3 Exit" << std::endl;
		std::cout << "Option (1 or 2 or 3) = ";
		std::cin.getline(line, 128);
		char option = line[0];

		if (option == '1')
		{
			DoLongMultiplication();
		}

		else if (option == '2')
		{
			DoKaratsuba();
		}

		else
			break;
	}

	return 0;
}

== Menu ==
1 Long Multiplication
2 Karatsuba Multiplication
3 Exit
Option (1 or 2 or 3) = 1
Enter base = 10
a = 506
b = 208
product = 105248
== Menu ==
1 Long Multiplication
2 Karatsuba Multiplication
3 Exit
Option (1 or 2 or 3) = 2
Enter base = 10
Enter m = 2
x1 = 5
x0 = 6
y1 = 2
y0 = 8
z = 105248
== Menu ==
1 Long Multiplication
2 Karatsuba Multiplication
3 Exit
Option (1 or 2 or 3) = 3

D:\Multiplication\x64\Debug\Multiplication.exe (process 30912) exited with code 0 (0x0).
Press any key to close this window . . .

Blog Entry (c) Tuesday February 3, 2026, by James Pate Williams, Jr. Derivation of the Classical Kinetic Energy Formula

Derivation of the Classical Kinetic Energy Equation from EinsteinDownload
// KineticEnergy.cpp
// Author: James Pate Williams, Jr. and NIST
// Copyright Monday February 2, 2026

#include <iomanip>
#include <iostream>
#include <vector>

const double c = 2.99792458E8;
const double mass0 = 1.67492750056E-27;

static double KineticEnergy(
    double v, size_t number,
    double& kilotonsTNT,
    std::vector<double>& terms)
{
    double c2 = c * c, ke = 0;;
    
    if (number >= 0)
        terms[0] = mass0 * c2;
    if (number >= 1)
        terms[1] = 0.5 * mass0 * v * v;
    if (number >= 2)
        terms[2] = 3.0 * mass0 * pow(v, 4.0) / (8.0 * c2);
    if (number >= 3)
        terms[3] = 225.0 * mass0 * pow(v, 6.0) / (720.0 * c2 * c2);

    for (size_t i = 0; i <= number; i++)
        ke += terms[i];

    kilotonsTNT = 2.3900573613767E-13 * ke;
    return ke;
}

int main()
{
    char line[128] = "";

    while (true)
    {
        double ke = 0.0, kilotonsTNT = 0.0, v = 0.0;
        size_t number = 0;

        std::cout << "Enter v = ";
        std::cin.getline(line, 128);
        v = atof(line);

        if (v == 0)
            break;

        std::cout << "# Terms = ";
        std::cin.getline(line, 128);
        number = static_cast<size_t>(atoi(line));
        std::vector<double> terms(number + 1);
        ke = KineticEnergy(v, number, kilotonsTNT, terms);
        std::cout << "KE = " << ke << " joules" << std::endl;
        std::cout << "KE = " << kilotonsTNT << " kilotons TNT";
        std::cout << std::endl;
        std::cout << std::scientific;

        for (size_t i = 0; i <= number; i++)
        {
            std::cout << "KE[" << i << "] = " << terms[i];
            std::cout << " joules" << std::endl;
            kilotonsTNT = 2.3900573613767E-13 * terms[i];
            std::cout << "KE[" << i << "] = " << kilotonsTNT;
            std::cout << " kilotons TNT" << std::endl;
        }
    }
}

Blog Entry (c) February 2, 2026, by James Pate Williams, Jr. and Especially the Microsoft Copilot Three Iterative Dichotomiser 3 (ID3) Examples

Blog Entry Monday February 2 2026Download

Blog Entry © Sunday, January 25, 2026, by James Pate Williams, Jr., Schwarzschild Solution of Einstein’s General Relativity Gravitational Field Equation

Blog Entry Sunday January 25 2026Download