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

Blog Entry © Friday, January 16, 2026, by James Pate Williams, Jr., Another Update of My Iowa Class Battleship Artillery Exterior Ballistics Application

Blog Entry Friday January 16 2026Download

Blog Entry © Wednesday, January 14, 2026, by James Pate Williams, Jr. Curvature of the Earth Table

Blog Entry Wednesday January 14 2026Download
// CurvatureOfTheEarth.cpp : Defines the entry point for the application.
//

#include "pch.h"
#include "framework.h"
#include "CurvatureOfTheEarth.h"
#include "GreatCircleDistance.h"
#include "Vincenty.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[128];                                // general purpose buffer
std::wstring outputText;                        // output wide character text

// 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_PTR CALLBACK    TableDialog(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_CURVATUREOFTHEEARTH, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    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_CURVATUREOFTHEEARTH));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_CURVATUREOFTHEEARTH);
    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)
{
    switch (message)
    {
    case WM_CREATE:
        DialogBox(hInst, MAKEINTRESOURCE(IDD_TABLE_DIALOG), hWnd, TableDialog);
        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);
            // 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;
}

INT_PTR CALLBACK TableDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    static double deltaTime = 0;
    static int height = 20, width = 80;
    static size_t delta = 0;
    static HFONT hFont = NULL;
    static HWND hEditMultiline = NULL;
    static GreatCircleDistance gcd;
    static Vincenty vincenty;

    switch (message)
    {
    case WM_INITDIALOG:
        hFont = CreateFont(
            -MulDiv(7, GetDeviceCaps(GetDC(hDlg), LOGPIXELSY), 72),
            0, 0, 0, FW_BOLD, FALSE, FALSE, FALSE,
            DEFAULT_CHARSET, OUT_DEFAULT_PRECIS,
            CLIP_DEFAULT_PRECIS, DEFAULT_QUALITY,
            FIXED_PITCH | FF_MODERN,
            TEXT("Courier New")
        );

        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_LEFT | ES_MULTILINE | ES_READONLY,
            10, 0, 10 * width, 30 * height,         // Position and size
            hDlg,                                   // Parent window handle
            (HMENU)IDC_EDIT_MULTILINE,              // Unique control ID
            hInst,                                  // Application instance
            NULL                                    // Extra parameter
        );

        SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, TRUE);
        return (INT_PTR)TRUE;

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

        if (LOWORD(wParam) == IDC_BUTTON_COMPUTE)
        {
            outputText += L"Range in Yards Versus Curvature of the Earth in Feet\r\n\r\n";
            outputText += L"Yards\t 0\t 100\t 200\t 300\t 400\t 500\t 600\t 700\t 800\t 900\r\n";
                
            for (double row = 1000.0; row <= 40000.0; row += 1000.0)
            {
                swprintf_s(line, L"%6.1lf\t", row);
                outputText += line;

                for (double col = 0.0; col <= 900.0; col += 100.0)
                {
                    double r = row + col;
                    double distance = 3.0 * gcd.Interpolate(r);

                    swprintf_s(line, L"%4.1lf\t", distance);
                    outputText += line;
                }

                outputText += L"\r\n";
            }

            outputText += L"Yards\t 0\t 100\t 200\t 300\t 400\t 500\t 600\t 700\t 800\t 900\r\n";
            SetWindowText(hEditMultiline, outputText.c_str());
            return (INT_PTR)TRUE;
        }
    }

    return (INT_PTR)FALSE;
}


#pragma once
#include "Vincenty.h"

struct PointRaR
{
    double ra;   // curvature of the Earth correction in yards
    double r;    // flat Earth distance (chord) in yards

    PointRaR(double ra, double r)
    {
        this->ra = ra;
        this->r = r;
    }
};

class GreatCircleDistance
{

private:

    Vincenty vincenty;
    std::vector<PointRaR> pts;

    bool binarySearch(double x, int& lt, int& rt);

public:

    friend int compare(
        const PointRaR& lt,
        const PointRaR& rt)
    {
        if (lt.r < rt.r)
            return -1;
        if (lt.r > rt.r)
            return +1;

        return 0;
    };

    // R is in yards, returns in yards
    // Curvature of the Earth correction

    double Interpolate(double R);

    // construction of the curvature of the Earth table

    GreatCircleDistance();
};


#pragma once

class Vincenty
{

public:

    static double Re; // Radius of Earth in meters
    double deltaSigma(
        double phi1, double lambda1,
        double phi2, double lambda2);
    double distance(
        double phi1, double lambda1,
        double phi2, double lambda2);
    double distanceConstantLambda(
        double phi1, double phi2);
    double distanceConstantXY(double z);
    double x(double phi, double lambda);
    double y(double phi, double lambda);
    double z(double phi);
    double phi(double z);
    double lambda(double x, double y);
};


#include "pch.h"
#include "GreatCircleDistance.h"

bool GreatCircleDistance::binarySearch(
    double x, int& lt, int& rt)
{
    int n = static_cast<int>(pts.size()), L = 0, M, R = n - 1;

Label10:

    if (x == pts[L].r)
    {
        lt = rt = L;
        return true;
    }

    if (x == pts[R].r)
    {
        lt = rt = R;
        return true;
    }

    if (x > pts[L].r && x < pts[R].r && R - L == 1)
    {
        lt = L;
        rt = R;
        return true;
    }

    if (x > pts[L].r && x < pts[R].r)
    {
        M = (L + R) / 2;

        if (x > pts[M].r)
        {
            L = M;
            goto Label10;
        }

        if (x < pts[M].r)
        {
            R = M;
            goto Label10;
        }
    }

    lt = rt = -1;
    return false;
}

double GreatCircleDistance::Interpolate(double R)
{
    int lt, rt;

    if (binarySearch(R, lt, rt))
    {
        double x0 = pts[lt].ra, x1 = pts[rt].ra;
        double y0 = pts[lt].r, y1 = pts[rt].r;
        double deltaX = x1 - x0, deltaY = y1 - y0;
        double ra = deltaX * (R - y0) / deltaY + x0;

        return ra;
    }

    return -1.0;
}

GreatCircleDistance::GreatCircleDistance()
{
    double deltaPhi = 0.000001, phi1 = 0.0, phi2 = deltaPhi, delta;
    double deltaRa, d0 = vincenty.z(phi1), d1, r, ra;
    int cnt = 0;

    pts.push_back(PointRaR(0.0, 0.0));

    while (cnt < 10000)
    {
        d1 = vincenty.z(phi2);

        if (d0 >= d1)
        {
            delta = d0 - d1;
            deltaRa = d0 * d0 - d1 * d1;
        }

        else
        {
            delta = d1 - d0;
            deltaRa = d1 * d1 - d0 * d0;
        }

        r = delta;
        ra = sqrt(deltaRa);
        ra = r >= ra ? r - ra : ra - r;
        pts.push_back(PointRaR(1.0936 * r, 1.0936 * ra));
        phi2 += deltaPhi;
        cnt++;
    }
}


#include "pch.h"
#include "Vincenty.h"

double Vincenty::Re = 6378137.0;	// radius of Earth in meters

double Vincenty::deltaSigma(
    double phi1, double lambda1,
    double phi2, double lambda2)
{
    double deltaPhi = phi1 - phi2, deltaLambda = lambda1 - lambda2;
    double cosPhi1 = cos(phi1), cosPhi2 = cos(phi2);
    double sinPhi1 = sin(phi1), sinPhi2 = sin(phi2);
    double cosDeltaLambda = cos(deltaLambda), sinDeltaLambda = sin(deltaLambda);
    double numer1 = cosPhi2 * sinDeltaLambda;
    double numer2 = cosPhi1 * sinPhi2 - sinPhi1 * cosPhi2 * cosDeltaLambda;
    double numer = sqrt(numer1 * numer1 + numer2 * numer2);
    double denom = sinPhi1 * sinPhi2 + cosPhi1 * cosPhi2 * cosDeltaLambda;

    return atan2(numer, denom);
}

double Vincenty::distance(
    double phi1, double lambda1,
    double phi2, double lambda2)
{
    return Re * deltaSigma(phi1, lambda1, phi2, lambda2);
}

double Vincenty::distanceConstantLambda(
    double phi1, double phi2)
{
    return Re * distance(phi1, 0.0, phi2, 0.0);
}

double Vincenty::distanceConstantXY(double z)
{
    return Re * distance(phi(0.0), 0.0, phi(z), 0.0);
}

double Vincenty::x(double phi, double lambda)
{
    return Re * sin(phi) * cos(lambda);
}

double Vincenty::y(double phi, double lambda)
{
    return Re * sin(phi) * sin(lambda);
}

double Vincenty::z(double phi)
{
    return Re * cos(phi);
}

double Vincenty::phi(double z)
{
    return acos(z / sqrt(Re));
}

double Vincenty::lambda(double x, double y)
{
    return asin(y / sqrt(x * x + y * y));
}

Blog Entry © Thursday, January 8, 2026, by James Pate Williams, Jr., Revised United States Navy Fast Battleship Iowa Class Artillery Ballistics Tables

Blog Entry Thursday January 8 2026Download

Blog Entry © Sunday, January 4, 2026, by James Pate Williams, Jr. Iterative Deepening A* Search to Solve the Fifteen Tile Puzzle (Win32 C/C++ Release x64 Configuration)

Blog Entry Sunday January 4 2026Download

Blog Entry © Thursday, January 1, 2026, by James Pate Williams, Jr., Win32 C/C++ Fast Battleship Class Iowa Ballistics Calculator (BB-61 Iowa, BB-62 New Jersey, BB-63 Missouri, BB-64 Wisconsin)

Blog Entry Thursday January 1 2026Download