Karatsuba algorithm too much recursion

Question

I'm trying to implement the Karatsuba multiplication algorithm in C ++, but right now I'm just trying to get it working in python.

Here is my code:

def mult(x, y, b, m): if max(x, y) < b: return x * y bm = pow(b, m) x0 = x / bm x1 = x % bm y0 = y / bm y1 = y % bm z2 = mult(x1, y1, b, m) z0 = mult(x0, y0, b, m) z1 = mult(x1 + x0, y1 + y0, b, m) - z2 - z0 return mult(z2, bm ** 2, b, m) + mult(z1, bm, b, m) + z0 

What I am not getting should be: how to create z2 , z1 and z0 ? Does mult function mult recursively correctly? If so, I am confused somewhere because the recursion does not stop.

Can someone please indicate where the error is?

Answer 1

NB: the answer below directly addresses the OP question about excessive recursion, but it does not try to provide the correct Karatsuba. Other answers are much more informative in this regard.

Try this version:

 def mult(x, y, b, m): bm = pow(b, m) if min(x, y) <= bm: return x * y # NOTE the following 4 lines x0 = x % bm x1 = x / bm y0 = y % bm y1 = y / bm z0 = mult(x0, y0, b, m) z2 = mult(x1, y1, b, m) z1 = mult(x1 + x0, y1 + y0, b, m) - z2 - z0 retval = mult(mult(z2, bm, b, m) + z1, bm, b, m) + z0 assert retval == x * y, "%d * %d == %d != %d" % (x, y, x * y, retval) return retval 

The biggest problem with your version is that your calculations of x0 and x1, as well as y0 and y1, are inverted. In addition, the algorithm output is not performed if x1 and y1 are 0, since in this case the factorization step becomes invalid. Therefore, you should avoid this possibility by ensuring that both x and y are greater than b ** m.

EDIT: fixed typo in the code; clarifications added

EDIT2:

To be clearer, commenting directly on the original version:

 def mult(x, y, b, m): # The termination condition will never be true when the recursive # call is either # mult(z2, bm ** 2, b, m) # or mult(z1, bm, b, m) # # Since every recursive call leads to one of the above, you have an # infinite recursion condition. if max(x, y) < b: return x * y bm = pow(b, m) # Even without the recursion problem, the next four lines are wrong x0 = x / bm # RHS should be x % bm x1 = x % bm # RHS should be x / bm y0 = y / bm # RHS should be y % bm y1 = y % bm # RHS should be y / bm z2 = mult(x1, y1, b, m) z0 = mult(x0, y0, b, m) z1 = mult(x1 + x0, y1 + y0, b, m) - z2 - z0 return mult(z2, bm ** 2, b, m) + mult(z1, bm, b, m) + z0 

Answer 2

Usually large numbers are stored as arrays of integers. Each integer is one digit. This approach allows you to multiply any number by the base power with a simple left-shift of the array.

Here is my list-based implementation (may contain errors):

 def normalize(l,b): over = 0 for i,x in enumerate(l): over,l[i] = divmod(x+over,b) if over: l.append(over) return l def sum_lists(x,y,b): l = min(len(x),len(y)) res = map(operator.add,x[:l],y[:l]) if len(x) > l: res.extend(x[l:]) else: res.extend(y[l:]) return normalize(res,b) def sub_lists(x,y,b): res = map(operator.sub,x[:len(y)],y) res.extend(x[len(y):]) return normalize(res,b) def lshift(x,n): if len(x) > 1 or len(x) == 1 and x[0] != 0: return [0 for i in range(n)] + x else: return x def mult_lists(x,y,b): if min(len(x),len(y)) == 0: return [0] m = max(len(x),len(y)) if (m == 1): return normalize([x[0]*y[0]],b) else: m >>= 1 x0,x1 = x[:m],x[m:] y0,y1 = y[:m],y[m:] z0 = mult_lists(x0,y0,b) z1 = mult_lists(x1,y1,b) z2 = mult_lists(sum_lists(x0,x1,b),sum_lists(y0,y1,b),b) t1 = lshift(sub_lists(z2,sum_lists(z1,z0,b),b),m) t2 = lshift(z1,m*2) return sum_lists(sum_lists(z0,t1,b),t2,b) 

sum_lists and sub_lists returns an abnormal result - a single digit may be larger than the base value. The normalize function solved this problem.

All functions expect to get a list of numbers in reverse order. For example, 12 in base 10 should be written as [2,1]. Let's take the square 9987654321.

 » a = [1,2,3,4,5,6,7,8,9] » res = mult_lists(a,a,10) » res.reverse() » res [9, 7, 5, 4, 6, 1, 0, 5, 7, 7, 8, 9, 9, 7, 1, 0, 4, 1] 

Answer 3

The goal of Karastuba multiplication is to improve the breeding and overcoming breeding algorithm by making 3 recursive calls instead of four. Therefore, the only lines in your script that should contain a recursive multiplication call are those that assign z0 , z1 and z2 . Everything else will be even worse. You cannot use pow to calculate b^m if you have not yet defined multiplication (and even more so elevation).

To do this, the algorithm critically uses the fact that it uses a positional notation system. If you have a representation of x numbers in the base b , then x*b^m simply obtained by shifting the bits of this representation m times to the left. This switching operation is substantially “free” with any positional notation. It also means that if you want to implement this, you must reproduce this positional notation and the “free” shift. Either you decided to calculate b=2 in the database and use the python bit operators (or the bit operators of the given decimal, hexadecimal, ... base, if they have a test platform), or you decided to implement for educational purposes what works for arbitrary b , and you reproduce this positional arithmetic with something like strings, arrays, or lists .

You already have a solution with lists . I like working with strings in python, since int(s,base) will give you an integer corresponding to the string s , considered as a numeric representation in the base base : this makes testing easier. I posted a heavily commented out string implementation as gist here , including the string for the number and the -to-string number for a good estimate.

You can test it by providing filled strings with a base and their (equal) length as mult arguments:

 In [169]: mult("987654321","987654321",10,9) Out[169]: '966551847789971041' 

If you do not want to determine the length of uppercase or highlighted lines, the fill function can do this for you:

 In [170]: padding("987654321","2") Out[170]: ('987654321', '000000002', 9) 

And of course it works with b>10 :

 In [171]: mult('987654321', '000000002', 16,9) Out[171]: '130eca8642' 

(Refer to wolfram alpha )

Answer 4

I believe that the idea of ​​the method is that the terms z _i are calculated using a recursive algorithm, but the results are not combined together in this way. Since the net result you want is

 z0 B^2m + z1 B^m + z2 

Assuming you choose the appropriate B value (say 2), you can calculate B ^ m without any multiplication. For example, when using B = 2, you can calculate B ^ m using bit shifts rather than multiplications. This means that the last step can be performed without any breeding.

One more thing - I noticed that you chose a fixed value of m for the whole algorithm. As a rule, you implement this algorithm, always having a value m, such that B ^ m is equal to half the number of digits in x and y when they are written in base B. If you use the powers of two, this will be done by choosing m = ceil (( log x) / 2).

Hope this helps!