Python Data Structures and Algorithms
上QQ阅读APP看书,第一时间看更新

Can we do better? A recursive approach

It turns out that in the case of long multiplication the answer is yes, there are in fact several algorithms for multiplying large numbers that require less operations. One of the most well-known alternatives to long multiplication is the Karatsuba algorithm, first published in 1962. This takes a fundamentally different approach: rather than iteratively multiplying single digit numbers, it recursively carries out multiplication operations on progressively smaller inputs. Recursive programs call themselves on smaller subsets of the input. The first step in building a recursive algorithm is to decompose a large number into several smaller numbers. The most natural way to do this is to simply split the number in to two halves, the first half of most significant digits, and a second half of least significant digits. For example, our four-digit number, 2345, becomes a pair of two-digit numbers, 23 and 45. We can write a more general decomposition of any 2 n digit numbers, x and y using the following, where m is any positive integer less than n:

So now we can rewrite our multiplication problem x, y as follows:

When we expand and gather like terms we get the following:

More conveniently, we can write it like this:

Where:

It should be pointed out that this suggests a recursive approach to multiplying two numbers since this procedure does itself involve multiplication. Specifically, the products ac, ad, bc, and bd all involve numbers smaller than the input number and so it is conceivable that we could apply the same operation as a partial solution to the overall problem. This algorithm, so far, consists of four recursive multiplication steps and it is not immediately clear if it will be faster than the classic long multiplication approach.

What we have discussed so far in regards to the recursive approach to multiplication, has been well known to mathematicians since the late 19th century. The Karatsuba algorithm improves on this is by making the following observation. We really only need to know three quantities: z2= ac ; z1=ad +bc, and z0= bd to solve equation 3.1. We need to know the values of a, b, c, d only in so far as they contribute to the overall sum and products involved in calculating the quantities z2, z1, and z0. This suggests the possibility that perhaps we can reduce the number of recursive steps. It turns out that this is indeed the situation.

Since the products ac and bd are already in their simplest form, it seems unlikely that we can eliminate these calculations. We can however make the following observation:

When we subtract the quantities ac and bd, which we have calculated in the previous recursive step, we get the quantity we need, namely (ad + bc):

This shows that we can indeed compute the sum of ad + bc without separately computing each of the individual quantities. In summary, we can improve on equation 3.1 by reducing from four recursive steps to three. These three steps are as follows:

  1. Recursively calculate ac.
  2. Recursively calculate bd.
  3. Recursively calculate (a +b)(c + d) and subtract ac and bd.

The following example shows a Python implementation of the Karatsuba algorithm:

    from math import log10  
def karatsuba(x,y):

# The base case for recursion
if x < 10 or y < 10:
return x*y

#sets n, the number of digits in the highest input number
n = max(int(log10(x)+1), int(log10(y)+1))

# rounds up n/2
n_2 = int(math.ceil(n / 2.0))
#adds 1 if n is uneven
n = n if n % 2 == 0 else n + 1

#splits the input numbers
a, b = divmod(x, 10**n_2)
c, d = divmod(y, 10**n_2)

#applies the three recursive steps
ac = karatsuba(a,c)
bd = karatsuba(b,d)
ad_bc = karatsuba((a+b),(c+d)) - ac - bd

#performs the multiplication
return (((10**n)*ac) + bd + ((10**n_2)*(ad_bc)))

To satisfy ourselves that this does indeed work, we can run the following test function:

    import random 
def test():
for i in range(1000):
x = random.randint(1,10**5)
y = random.randint(1,10**5)
expected = x * y
result = karatsuba(x, y)
if result != expected:
return("failed")
return('ok')