Day 10 - Karatsuba multiplication


Definition(s)

The Karatsuba algorithm is a fast multiplication algorithm.

It reduces the multiplication of two n-digit numbers to at most ${\displaystyle n^{\log _{2}3}\approx n^{1.585}}$ single-digit multiplications in general.


In [1]:
import numpy as np  # used for generating random numbers

Convertion utility functions


In [2]:
def int_to_big(x):
    if x == 0:
        return [0]

    z = []

    while x > 0:
        t = x % 10
        z.append(t)
        x //= 10

    trim(z)
    return z


def big_to_int(x):
    z, p = 0, 1

    for d in x:
        z += p * d
        p *= 10

    return z

Multiplication utility functions


In [3]:
from itertools import zip_longest


def trim(z):
    while len(z) > 1 and z[-1] == 0:
        z.pop(-1)


def add(x, y):
    z, carry = [], 0

    for r, s in zip_longest(x, y, fillvalue=0):
        carry += r + s
        z.append(carry % 10)
        carry //= 10
        
    if carry:
        z.append(carry)

    return z

def subtract(x, y):
    z, carry = [], 0

    for r, s in zip_longest(x, y, fillvalue=0):
        carry += r - s
        z.append(carry % 10)
        carry //= 10
        
    trim(z)
    return z

Karatsuba's algorithm


In [4]:
def karatsuba(x, y):
    # ensure same length
    while len(x) < len(y):
        x.append(0)

    while len(x) > len(y):
        y.append(0)

    # length
    n = len(x)
    half = n // 2

    if n == 1:
        return add([x[0] * y[0]], [])
    
    # cut-off for improved efficiency
    if n <= 50:
        a = big_to_int(x)
        b = big_to_int(y)
        z = a * b
        return int_to_big(z)

    x0, x1 = x[:half], x[half:]
    y0, y1 = y[:half], y[half:]

    # x = x0x1
    # y = y0y1

    # z0 = x0 * y0
    # z1 = x1 * y1
    # z2 = (x0 + x1) * (y0 + y1)
    # z2 = z2 - (z0 + z1)

    z0 = karatsuba(x0, y0)
    z1 = karatsuba(x1, y1)
    z2 = karatsuba(add(x0, x1), add(y0, y1))
    z2 = subtract(z2, add(z0, z1))

    z = add(z0, [0] * (half << 1) + z1)
    z = add(z, [0] * half + z2)

    return z

Multiplication and testing


In [5]:
def multiply(x, y):
    xb = int_to_big(x)
    yb = int_to_big(y)
    
    zb = karatsuba(xb, yb)

    return big_to_int(zb)


def test(x, y):
    z = multiply(x, y)
    assert x * y == z
    print("{} x {} = {}".format(x, y, z))

Generate big integers


In [6]:
def gen_long(n):
    x = ''.join(map(str, np.random.randint(0, 10, n)))
    return int(x)

Run(s)


In [7]:
test(1432423423420, 12321312332131233)
test(8931283129323420, 1233123602345430533)


1432423423420 x 12321312332131233 = 17649336391818484838565676860
8931283129323420 x 1233123602345430533 = 11013376025998265385094936899982860

In [8]:
tests = 30
for _ in range(tests):
    n = np.random.randint(1, 15)
    x, y = gen_long(n), gen_long(n)
    test(int(x), int(y))


755 x 56 = 42280
276 x 5194 = 1433544
8771335752925 x 3258851870994 = 28584483929536201958157450
44640241182 x 985490183003 = 43992519451747237029546
3172239478 x 4188169149 = 13285875514999464222
2357279 x 8530917 = 20109751494843
3133832 x 9434478 = 29566069059696
79790 x 92703 = 7396772370
9076 x 1363 = 12370588
965055 x 721654 = 696435800970
86241813 x 75117433 = 6478263609826029
738 x 402 = 296676
763596619269 x 2216534026665 = 1692537889256097498807885
14488431985761 x 25110754122358 = 363815453212951534707744438
34 x 70 = 2380
1525183 x 2119398 = 3232469799834
521791938414 x 441693054302 = 230471874988240741757028
314061344195 x 180725348526 = 56758845888185421906570
35080644461653 x 88581266358920 = 3107487911100256102874494760
23133 x 89268 = 2065036644
4991 x 5856 = 29227296
3510 x 4347 = 15257970
7021024 x 2249169 = 15791469529056
499860787 x 854644335 = 427203189898191645
6429628 x 7129417 = 45839499166876
322048556457 x 901771307403 = 290414147803477747551171
93592 x 47401 = 4436354392
58888 x 57666 = 3395835408
712227445386 x 838125772789 = 596936178065676542401554
1686242329 x 63929711532 = 107800985666017838028

In [9]:
%%time

a, b = gen_long(1000), gen_long(1000)

z = multiply(a, b)
assert z == a * b


CPU times: user 32 ms, sys: 0 ns, total: 32 ms
Wall time: 33.6 ms

In [10]:
%%time

a, b = gen_long(20000), gen_long(20000)

z = multiply(a, b)
assert z == a * b


CPU times: user 4.77 s, sys: 4 ms, total: 4.77 s
Wall time: 4.77 s

Karatsuba multiplication using Baruchel's implementation

Karatsuba's algorithm is already implemented in Python. Check this package.


In [11]:
from karatsuba import *

In [12]:
def power_of_two(x):
    p = 1
    while p < x:
        p <<= 1

    return p


def reverse(num):
    return int(str(num)[::-1])

In [13]:
def kat_multiply(x, y):
    if x == 0 or y == 0:
        return 0

    xs = list(map(int, str(x)))
    ys = list(map(int, str(y)))

    n = power_of_two(max(len(xs), len(ys)))
    plan = make_plan(range(n), range(n))

    xs = [0] * (n - len(xs)) + xs
    ys = [0] * (n - len(ys)) + ys

    zs = plan(xs, ys)
    zs.pop(-1)
    
    zs = list(reversed(zs))

    while zs[-1] == 0:
        zs.pop(-1)

    ans = 0
    for p, d in enumerate(zs):
        ans += d * 10 ** p

    return ans

In [14]:
tests = 30
for _ in range(tests):
    n = np.random.randint(1, 15)
    x, y = gen_long(n), gen_long(n)
    
    z = kat_multiply(x, y)
    
    assert z == x * y
    print("{} x {} = {}".format(x, y, z))


741760 x 262930 = 195030956800
430408061 x 471695629 = 203021601060065369
2 x 2 = 4
2591169055 x 289097656 = 749100900100235080
80806 x 1670 = 134946020
2706508 x 3428885 = 9280304683580
868464280941 x 454843598229 = 395015418476565586053489
52 x 12 = 624
898196568 x 774385300 = 695550218769650400
91557 x 30470 = 2789741790
1024262 x 90070663 = 92255957425706
925782549 x 582219748 = 539008882381577652
1699487528480 x 309157175085 = 525408763397065283920800
37633788 x 67889859 = 2554952560955892
7559870 x 4880346 = 36894781315020
3070 x 8061 = 24747270
648467 x 281305 = 182417009435
974553 x 2186880 = 2131230464640
855 x 993 = 849015
9565 x 892 = 8531980
89149 x 2457 = 219039093
30 x 61 = 1830
895 x 168 = 150360
6 x 0 = 0
61850 x 96746 = 5983740100
33 x 89 = 2937
859900676683 x 7158636465867 = 6155716341126632932279161
8422334556886 x 2032395236742 = 17117512635662649636305412
2457187 x 1447067 = 3555714220529
374840 x 380954 = 142796797360

In [15]:
%%time

a, b = gen_long(100), gen_long(100)

z = kat_multiply(a, b)
assert z == a * b


CPU times: user 888 ms, sys: 20 ms, total: 908 ms
Wall time: 908 ms