Marina Mele's site

3rd-grade & Karatsuba multiplication Algorithms

In this post we’re going to study the third-grade algorithm to compute the product of two numbers, and we’re going to compare it with a much more efficient algorithm: The Karatsuba multiplication algorithm.

Did you know that these two algorithms are the ones used in the built in Python multiplication?

We will talk about the Order of both algorithms and give you Python implementations of both of them.

Maybe you’re thinking, why is she writing now about algorithms? Some time ago, I took the course Algorithms: Design and Analysis, Part 1&2, by Tim Roughgarden,  and it was a lot of fun.

Now I’m taking the course again, but I’m spending more time to review the algorithms and play with them. I really encourage you to take the course if you have some time. But first, read this post

This is the outline:

• Third grade multiplication algorithm
• Python code for the third grade product algorithm
• The Karatsuba Algorithm
• Python code for the Karatsuba algorithm

Let’s start!

Third grade multiplication algorithm

First, we’re going to review the third grade algorithm, which all of you already know

Let’s start with these two numbers: 5678 x 1234. In order to compute their product you start with 4*5678, represented as:

  (2)(3)(3)
   5  6  7  8
x  1  2  3 |4|
-------------
2  2  7  1  2

Let’s count the number of operations performed in this step:

– If n is the number of digits of the first number, there are n products and at most n sums (carried numbers).

– So in total, you need 2n operations.

If we continue with the product, we have to repeat n times the same operation (where we assume that n is also the number of digits of the second number):

         5  6  7  8
      x  1  2  3  4
      -------------
      2  2  7  1  2
   1  7  0  3  4
1  1  3  5  6
5  6  7  8

which gives a total of 2n² operations.

Finally, you need to sum all these numbers,

            5  6  7  8
         x  1  2  3  4
         -------------
         2  2  7  1  2
      1  7  0  3  4
   1  1  3  5  6
+  5  6  7  8
----------------------
   7  0  0  6  6  5  2

which takes around n² operations more.

So in total, the number of operations is ~3n², which means that it’s quadratic with the input (proportional to n²).

One point I would like to mention about this algorithm is that it’s complete: no matter what x and y you start with, if you perform correctly the algorithm, the algorithm terminates and finds the correct solution.

Python code for the third grade product algorithm

Here I give you an example of an implementation of the third grade algorithm, where I have included a Counter for the sum and product operations:

# third_grade_algorithm.py

def counted(fn):
    # Counter Decorator
    def wrapper(*args, **kwargs):
        if "" in args or " " in args:
            return "".join(map(lambda s: s.strip(), args))
        wrapper.called += 1
        return fn(*args, **kwargs)
    wrapper.called = 0
    wrapper.__name__ = fn.__name__
    return wrapper


@counted
def prod(x, y):
    # x, y are strings --> returns a string of x*y
    return str(eval("%s * %s" % (x, y)))


@counted
def suma(x, y):
    # x, y are strings --> returns a string of x+y
    return str(eval("%s + %s" % (x, y)))


def one_to_n_product(d, x):
    """d is a single digit, x is n-digit --> returns a string of d*x
    """
    result = ""
    carry = "0"
    for i, digit in enumerate(reversed(x)):
        r = suma(prod(d, digit), carry)
        carry, digit = r[:-1], r[-1]
        result = digit + result
    return carry + result


def sum_middle_products(middle_products):
    # middle_products is a list of strings --> returns a string
    max_length = max([len(md) for md in middle_products])
    for i, md in enumerate(middle_products):
        middle_products[i] = " " * (max_length - len(md)) + md
    carry = "0"
    result = ""
    for i in range(1, max_length + 1):
        row = [carry] + [md[-i] for md in middle_products]
        r = reduce(suma, row)
        carry, digit = r[:-1], r[-1]
        result = digit + result
    return carry + result


def algorithm(x, y):
    # x, y are integers --> returns an integer, x*y
    x, y = str(x), str(y)
    middle_products = []
    for i, digit in enumerate(reversed(y)):
        middle_products.append(one_to_n_product(digit, x) + " " * i)
    return int(sum_middle_products(middle_products))

Using that algorithm, if you run

$ python -i third_grade_algorithm.py

where (third_grade_algorithm.py is the name of the file), you will run the previous code and terminate with the Python console open. This way you can call your algorithm function and try:

>>> algorithm(6885, 1600)
1101600
>>> print("Suma was called %i times" % suma.called)
Suma was called 20 times 
>>> print("Prod was called %i times" % prod.called)
Prod was called 16 times

So it took 20 + 16 = 36 operations to compute the product of these two numbers.

Once we have this code, we can average the number of operations used in the product of n-digit numbers:

def random_prod(n):
    suma.called = 0
    prod.called = 0
    x = randint(pow(10, n - 1), pow(10, n))
    y = randint(pow(10, n - 1), pow(10, n))
    algorithm(str(x), str(y))
    return suma.called + prod.called


def average():
    ntimes = 200
    nmax = 10
    result = []
    for n in range(nmax):
        avg = sum([random_prod(n + 1) for i in range(ntimes)]) / float(ntimes)
        result.append([n + 1, avg])
    return result

In the following figure, we plot the result of the previous experiment when ntimes = 200 samples and nmax = 10 digits:

We can see that these points fit the curve f(n) = 2.7 n².

Note that 2.7 < 3, the proportionality factor we deduced before. This is because we were assuming the worst case scenario, whereas the average product takes into account random numbers.

But can we do better? Is there an algorithm that performs the product of two numbers quicker?

The answer is yes, and in the following section we will study one of them: The Karatsuba Algorithm.

The Karatsuba Algorithm

Let’s first explain this algorithm through an example:

Imagine you want to compute again the product of these two numbers:

x = 5678
y = 1234

In order to do so, we will first decompose them in a singular way:

x = 5678
a = 56; b = 78
--> x = 100 * a + b

and the same for y:

y = 1234
c = 12; d = 34
--> y = 100 * c + d

If we want to know the product x * y:

xy = (100 * a + b)(100 * c + d) = 100^2 * ac + 100 * (ad + bc) + bd

Where we can calculate the following 3 parts separately:

A = ac 
B = ad + bc
C = bd

However, we need two compute two products for the B term. Let’s see how we can compute only one product to reduce calculations

If we expand the product:

D = (a+b) * (c+d) = ac + ad + bc + bd

we see than the right hand side of the equation contains all A,B and C terms.

In particular, if we isolate B = ad + bc, we get:

B = D - ac - bd = D - A - C

Great! Now we only need three smaller products in order to compute x * y:

xy = 100^2 A + 100(D - A - C) + C
A = ac
B = (a+b)(c+d)
C = bd

Let’s put some numbers into the previous expressions:

xy = 100^2 * A + 100 * (D - A - C) + C
A = ac = 56 * 12 = 672
C = bd = 78 * 34 = 2652
D = (a+b)(c+d) = (56 + 78)(12 + 34) = 134 * 46 = 6164

xy = 100^2 * 672 + 100 * (6164 - 672 - 2652) + 2652 
   = 6720000 + 284000 + 2652
   = 7006652

Yes! the result of 1234 * 5678 = 7006652!

In the next section, we’ll see the pseudo code for the Karatsuba Algorithm and we’re going to translate it into Python code.

Python code for the Karatsuba Algorithm

In general though, we don’t have 4-digit numbers but n-digit numbers. In that case, the decomposition to be made is:

x = n-digit number
m = n/2 if n is even
m = (n+1)/2 if n is odd

a = 10^m * x1 + x2
--> x1 = First (n-m) digits of x
--> x2 = Last m digits of x

which is slightly different if n is even or odd.

In pseudocode, the Karatsuba Algorithm is:

procedure karatsuba(x, y)
  /* Base Case */
  if (x < 10) or (y < 10)
    return x*y
  /* calculates the number of digits of the numbers */
  m = max(size_base10(x), size_base10(y))
  m2 = m/2
  /* split the digit sequences about the middle */
  a, b = split_at(x, m2)
  c, d = split_at(y, m2)
  /* 3 products of numbers with half the size */
  A = karatsuba(a, c)
  C = karatsuba(b, d)
  D = karatsuba(a+b, c+d)
  return A*10^(2*m2) + (D-A-C)*10^(m2) + C

And the order of this algorithm Ο(n^log₂3) ≈ Ο(n^1.585) is (more info).

Let’s write now the Python code

def karatsuba(x, y):
    # Base case
    if x < 10 or y < 10:
        return x * y
    # Calculate the number of digits of the numbers
    sx, sy = str(x), str(y)
    m2 = max(len(sx), len(sy)) / 2
    # Split the digit sequences about the middle
    ix, iy = len(sx) - m2, len(sy) - m2
    a, b = int(sx[:ix]), int(sx[ix:])
    c, d = int(sy[:iy]), int(sy[iy:])
    # 3 products of numbers with half the size
    A = karatsuba(a, c)
    C = karatsuba(b, d)
    D = karatsuba(a + b, c + d)
    return A * 10**(2 * m2) + (D - A - C) * 10**m2 + C

assert(karatsuba(1234, 5678) == 7006652)

Which is much simpler than the third-grade algorithm right?

Moreover, if we can compare the number of calculations of this algorithm with respect to the third-grade one:

where the red line is the Karatsuba algorithm, and the blue lines the 3rd grade one (see above).

Finally, I want to write here one of the quotes mentioned in the Algorithms course, which might encourage you to find a better solution

The post 3rd-grade & Karatsuba multiplication Algorithms appeared first on Marina Mele's site.