Toom-k Polynomial Multiplication

I’ve been pretty busy with classes the past two weeks, but I’ve learned a few neat things and I’d like to share one of them.
Multiplying big numbers is a problem when the numbers are really, really big. How big is really, really big? Depends on the hardware your using.

But to multiply really big numbers, we represent the two multiplicands as polynomials where each term is in the form cxⁿ, where x is the base in the number system, c is its coefficent modifier and n is the power to which the base is modified. For example, we can represent 4587 as 4x³ + 5x² + 8x¹ + 7x⁰, where x = 10.
Now, two numbers in this polynomial base representation can be multiplied to give the result. Here’s a small example:

321 = 3*10^2 + 2*10^1 + 1*10^0 ==> 3x^2 + 2x + 1

123 = 1*10^2 + 2*10^1 + 3*10^0 ==> 1x^2 + 2x + 3

321 * 123 = (3x^2 + 2x + 1)*(1x^2 + 2x + 3) = 3x^4 + 6x^3 + 9x^2 + 2x^3 + 4x^2 + 6x + x^2 +2x + 3 
= 3x^4 + 8x^3 + 14x^2 + 8x + 3

Now plugging in x = 10 ==>  3(10)^4 + 8(10)^3 + 14(10)^2 + 8(10) + 3 = 39483

Now verify with a calculator that 123*321 = 39483

Cool, huh?

We can use this knowledge to multiply some really, really big numbers. Since we have a way to calculate them product using polynomials, the question now is this: how do we multiply these polynomials when the order is large (like order-50 or order-n)? This is where the Toom-k algorithm becomes helpful.
First, let’s represent the polynomials as arrays where the position i indicates the power and the value at i represents the coefficient.

3x^2 + 2x + 1

will be represented as

[1, 2, 3]

Toom-k states that we can divide the polynomial array into sub-arrays of length (n/k) and perform (2*k)-1 recursive calls on those sub-arrays, then recombine them in O(n) time, and have a solution. For this article I want to focus on Toom-3, or otherwise called the Toom-Cook algorithm. To simplfiy the algorithm, we’ll require the length of our input arrays to be divisible by 3. Call the two input polynomials P and Q. P and Q get partitioned into n/3 or otherwise stated, they are divided into thirds. Call these sub-arrays A, B, C, D, E, and F where these are each of size n/3. Create 5 more arrays of size n/3 called G, H, J, L, M. Their values are as follows:
G = A+C, H = D+F, J = G-B, K = H-E, L = 2(J+A)-C, M = 2(K+D)-F
Here comes some serious math.
In this method of representing polynomials as arrays, an order-n polynomial takes up an array of length (n+1). Like in our previous example:

3x^2 + 2x + 1 is order-2

so n = 2, but its array representation takes up (n+1) = 3 length

[1, 2, 3]

Since this is known, then we can determine that multiplying two polynomials together of any order-n will generate in our system of representation an array of length (2*n)-1.
Ok, back to the Toom-3. We now create 5 more arrays of length (2*n)-1 named R, S, T, U, V. They are defined as R = AD, S = CF, T = (G+B)(H+E), U = JK, V = LM.
One more set of arrays of length (2*n/3)-1 named W, X, Y, Z and are defined as W = (V-T)/3, X = (T-U)/2, Y = U-S, Z = (Y-W)/2 + 2R.
Now we need to recombine these arrays into the final solution form. The product P*Q is defined as follows:
P*Q = Rx^(4n/3) + Zx^n + (X+Y-R)x^(2n/3) + (X-Z)x^(n/3) + S
This formula can look confusing, but it makes sense with some code to go with it:

BigInt* make_array(BigInt size)
{
	BigInt*res = new BigInt[size];
	for(BigInt i = 0; i < size; i ++)
	{
		res[i] = 0;
	}
	return res;
}
BigInt* Toom3(BigInt* P, BigInt*Q, BigInt n)
{
	BigInt pq_size = (2*n) -1 ;
	BigInt *PQ_Res = make_array(pq_size);
	BigInt sub_eq_size = n/3;
	BigInt *A = new BigInt[sub_eq_size], *B =new BigInt[sub_eq_size], *C=new BigInt[sub_eq_size], *D=new BigInt[sub_eq_size], *E=new BigInt[sub_eq_size], *F=new BigInt[sub_eq_size];
	memcpy(C,&P[0], sizeof(BigInt) * sub_eq_size);
	memcpy(B,&P[sub_eq_size], sizeof(BigInt) * sub_eq_size);
	memcpy(A,&P[sub_eq_size*2], sizeof(BigInt) * sub_eq_size);
	memcpy(F,&Q[0], sizeof(BigInt) * sub_eq_size);
	memcpy(E,&Q[sub_eq_size], sizeof(BigInt) * sub_eq_size);
	memcpy(D,&Q[sub_eq_size*2], sizeof(BigInt) * sub_eq_size);

	BigInt *G = new BigInt[sub_eq_size], *H =new BigInt[sub_eq_size], *J =new BigInt[sub_eq_size], *K =new BigInt[sub_eq_size], *L =new BigInt[sub_eq_size],*M =new BigInt[sub_eq_size];
	BigInt* GpB=new BigInt[sub_eq_size], *HpE=new BigInt[sub_eq_size];
	for(BigInt i = 0; i < sub_eq_size; i++)
	{
		G[i] = A[i] + C[i];
		H[i] = D[i] + F[i];
		J[i] = G[i] - B[i];
		K[i] = H[i] - E[i];
		L[i] = (2* (J[i] + A[i])) - C[i];
		M[i] = (2* (K[i] + D[i])) - F[i];
		GpB[i] = G[i] + B[i];
		HpE[i] = H[i] + E[i];
	}
	BigInt *R, *S, *T, *U, *V; //finally!
	R = mult(A, D, sub_eq_size, sub_eq_size);
	S = mult(C, F, sub_eq_size, sub_eq_size);
	U = mult(J, K, sub_eq_size, sub_eq_size);
	V = mult(L, M, sub_eq_size, sub_eq_size);
	T = mult(GpB, HpE, sub_eq_size, sub_eq_size);
	
	BigInt *W=new BigInt[2*n/3 -1], *X=new BigInt[2*n/3 -1], *Y=new BigInt[2*n/3 -1], *Z=new BigInt[2*n/3 -1]; // i just thought i was done.
	for(BigInt i = 0; i < 2*n/3 -1; i++)
	{
		W[i] = (V[i] - T[i])/3;
		X[i] = (T[i] - U[i])/2;
		Y[i] = U[i] - S[i];
		Z[i] = ((Y[i] - W[i])/2) + (2*R[i]);
	}
	for(BigInt i = 0; i < 2*n/3 -1; i++)
	{
		PQ_Res[i] 		+= S[i];
		PQ_Res[i + (n/3)]	+= (X[i] - Z[i]);
		PQ_Res[i + (2*n/3)] 	+= (X[i] + Y[i] - R[i]);
		PQ_Res[i + (n)] 	+= Z[i];
		PQ_Res[i + (4*n)/ 3] 	+= R[i];
	}
	return PQ_Res;
}

This method can be analyzed in its runtime using the Master Recurrence Theorem. Take Toom-3 for example. Its running time is T(n) = θ(1) + 5T(n/3) + θ(n). Using the Recurrence theorem we can show that the running time of Toom-3 is n^log₃5.For any general Toom-k, it can be shown that since we do (2*k) – 1 recursions and we split the arrays in to sub-arrays of size (n/k), the running time of Toom-k is T(n) = θ(1) + ((2*k)-1)T(n/k) + θ(n) which when used against the Recurrence theorem we show that it is of n^log_k(2*k)-1 time complexity. I found this very interesting and I hope you do too!!