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;
}
Numerical Explorations 