Karatsuba multiplication with integers of size 3

I understand how to apply Karatsuba multiplication in 2 digit integers.

$$\begin{array}
\quad & \quad & x & y \\
\times & \quad & z & w \\
\hline
\quad &?&?&?
\end{array}$$

$$\begin{align}
i & = zx \\
ii & = wy \\
iii & = (x+y)(w+z)
\end{align}$$

and the result is
$$i \cdot 10^2n + [(iii – ii – i) \cdot 10^n] + ii$$

So then my question is, how do I apply this in 3 digit integers?

i.e., this is the problem:

I have 2 three digit numbers and I want to split this into 5 multiplications and they’re of size $n/3$. I’m pretty sure I have to use Karatsuba multiplication, but I’m not entirely sure… I’m pretty new to here so my phrasing might be very poor. I hope I conveyed myself well enough. Thank you for any help.

Solutions Collecting From Web of "Karatsuba multiplication with integers of size 3"

Well I’m also pretty new to multiplication and I’ve stumbled on your question after looking for multiplication algorithm over solving multiplication of two 3-digit numbers with only 5 multiplications.

So…

  1. First thing – with Karatsuba you probably can’t do this. Deriving from 2-digit version of algorithm:
    $$
    (10a + b) (10x + y) = 100ax + by + 10\Big( (a+b)(x+y) – ax – by\Big)
    $$
    and applying it recursively we’ll get:
    $$
    \Big(
    \; 10 (10a + b) + c \;
    \Big) \cdot
    \Big(
    \; 10 (10x + y) + z \;
    \Big) =
    100 (10a+b)(10x+y) \;+\; cz \;+\;
    10\Big(
    (10a+b+c)(10x+y+z) \; – \; (10a+b)(10x+y) \; – \; cz
    \Big) =
    100
    \Bigg(
    100ax + by + 10 \Big( (a+b)(x+y) – ax – by \Big)
    \Bigg) \;+\; cz \;+\;
    10\Bigg\{
    100ax \;+\; (b+c)(y+z) \;+\;
    10\bigg[
    \; (a+b+c)(x+y+z) \;-\; ax \;-\; (b+c)(y+z) \;
    \bigg] \;-\;
    \bigg[ \;
    100ax \;+\; by \;+\;
    10\Big(
    (a+b)(x+y) – ax – by
    \Big)
    \bigg] \;-\; cz
    \Bigg\} \;\;\;\; =
    \;\;\;\;10^4ax +10^2by +10^3(a+b)(x+y)-10^3ax-10^3by+cz+10^3ax+10(b+c)(y+z)+10^2(a+b+c)(x+y+z)-10^2ax-10^2(b+c)(y+z)-10^3ax-10by-10^2(a+b)(x+y)+10^2ax+10^2by-10cz \;\;\;\; =
    \;\;\;\;
    10^4ax +10^3\Big( (a+b)(x+y)-ax-by \Big) + 10^2\Big( (a+b+c)(x+y+z)-(b+c)(y+z)-(a+b)(x+y)+2by\Big)+10\Big((b+c)(y+z)-by-cz\Big)+cz
    $$
    So it reduces number of necessary multiplications to 6:

    • ax
    • by
    • cz
    • (a+b)(x+y)
    • (b+c)(y+z)
    • (a+b+c)(x+y+z)
  2. But it should be possible to do this with Toom-Cook multiplication algorithm.

What is idea behind that? Simply – instead of doing calculations we treat each set of numbers as polynomial:
$$
p(t) = at^2 + bt + c
$$
$$
q(t) = xt^2 + yt + z
$$

And we are looking for polynomial:

$$
w(t) = p(q) \cdot q(t) = w_4 t^4 + w_3 t^3 + w_2 t^2 + w_1 t + w_0
$$

As it can be easily seen each coefficient $w_i$ of our $w(t)$ polynomial is one of products from multiplication. To avoid multiplication of polynomials (which would kill the case) we will try to calculate coefficients of $w(t)$ by using interpolation on our $p(t)$ and $q(t)$ polynomials. To do so we need some interpolation points. For 5 variables ($w_0,w_1,w_2,w_3,w_4$) we need 5 points. Let those be (as in wikipedia article) $t \in \{-2,-1,0,1,\infty\} $.

Points of polynomial evaluation may be whatever you wish. But it is common to chose $0$ and $\infty$ among them. In case of polynomials value at $\infty$ doesn’t really make sense so instead of this we are using
$$
p(\infty) = \lim_{t \to \infty} {{p(t)} \over {t^{deg(p)}}}
$$

And thus it is always coefficient of highest power (like $w_4$ in case of $w(t)$)

Having this explained let’s try it:

$$
\begin{array}[b]{l}
t=0 &&p(0) = c \\
&&q(0) = z \\
t=1 &&p(1) = a +b +c \\
&&q(1) = x +y +z \\
t=-1 &&p(-1) = a -b +c \\
&&q(-1) = x -y +z \\
t=-2 &&p(-2) = 4a-2b +c \\
&&q(-2) = 4x-2y +z \\
t=\infty &&p(\infty) = a \\
&&q(\infty) = x
\end{array}
$$
Which leads us to:
$$
\begin{array}[b]{l}
w(0) &&= w_0 &&= cz \\
w(1) &&= w_4 +w_3+w_2+w_1+w_0 &&= (a+b+c)(x+y+z) \\
w(-1) &&= w_4 -w_3+w_2-w_1+w_0 &&= (a-b+c)(x-y+z) \\
w(-2) &&= 16w_4-8w_3+4w_2-2w_1+w_0 &&= (4a-2b+c)(4x-2y+z)\\
w(\infty) &&= w_4 &&= ax
\end{array}
$$
So we see that there is only needed 5 multiplications of those specific sums and we’re almost there.

To speed things up let’s use matrix representation of equations above:

$$
\begin{pmatrix}
w(0) \\
w(1) \\
w(-1)\\
w(-2)\\
w(\infty)
\end{pmatrix} = \begin{pmatrix}
1 && 0 && 0 && 0 && 0 \\
1 && 1 && 1 && 1 && 1 \\
1 && -1 && 1 && -1 && 1 \\
16 && -8 && 4 && -2 && 1 \\
0 && 0 && 0 && 0 && 1
\end{pmatrix} \cdot \begin{pmatrix}
w_0 \\
w_1 \\
w_2 \\
w_3 \\
w_4
\end{pmatrix}
$$

Using inverse of this matrix (taken directly from Wikipedia):

$$
\begin{pmatrix}
w_0 \\
w_1 \\
w_2 \\
w_3 \\
w_4
\end{pmatrix} =
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 \\
1/2 & 1/3 & -1 & 1/6 & -2 \\
-1 & 1/2 & 1/2 & 0 & -1 \\
-1/2 & 1/6 & 1/2 & -1/6 & 2 \\
0 & 0 & 0 & 0 & 1
\end{pmatrix} \cdot \begin{pmatrix}
w(0) \\ w(1) \\ w(-1) \\ w(-2) \\ w(\infty)
\end{pmatrix}
$$
Which leads to:
$$
\begin{array}[b]{l}
w_0 = cz \\
w_1 = \frac {w(0)} 2 + \frac {w(1)} 3 – w(-1) + \frac {w(-2)} 6 -2 w(\infty) \\
w_2 = -w(0) + \frac {w(1)} 2 + \frac {w(-1)} 2 – w(\infty) \\
w_3 = -\frac{w(0)} 2 + \frac{w(1)} 6 + \frac{w(-1)} 2 – \frac{w(-2)} 6 + 2 w(\infty) \\
w_4 = ax
\end{array}
$$
And finally:
$$
\begin{array}[b]{l}
w_0 = cz \\
w_1 = \frac {cz} 2 + \frac {(a+b+c)(x+y+z)} 3 – (a-b+c)(x-y+z) + \frac {(4a-2b+c)(4x-2y+z)} 6 -2 ax \\
w_2 = -cz + \frac {(a+b+c)(x+y+z)} 2 + \frac {(a-b+c)(x-y+z)} 2 – ax \\
w_3 = -\frac{cz} 2 + \frac{(a+b+c)(x+y+z)} 6 + \frac{(a-b+c)(x-y+z)} 2 – \frac{(4a-2b+c)(4x-2y+z)} 6 + 2 ax \\
w_4 = ax
\end{array}
$$

Edit:

For those looking to use those equations as support for some programming (as myself – I needed this for my own 96bit integer class):

  1. Take overflow into account – especially with addition as it may be easily overlooked but will happen – especially when adding 3 segments with multipliers like 16 etc.

  2. Quite important cost (algorithmically speaking) is division here – although it will be “exact division” (see in GMP math library) so it means that it will always come up without any reminder (so it’s cheaper) and it can be precalculated (either by compiler or you in assembly or whatever)

  3. You may try to experiment with different interpolation points to get “more fitting” division/multiplication numbers – although going any further than 8/-8 seems impractical – you will fast hit overflow limits or hurt precision (depends if you are going for either integer or floating point numbers)

PS. Please point me out any mistakes I’ve made – both in english as well in equations