Recently I thought of an algorithm for multiplication and decided to stop dreaming and start writing on paper my ideas, and even implement this to code (in this case – Python 3.9.1).
I do not know if it resembles Karatsuba’s algorithm, but I glanced at it and it seems to work very differently.
The idea behind this multiplication (calculating $x cdot y$) algorithm is to represent them as a power of two, plus some remainder, then use the distributive rule of multiplication to get:
$$x = 2^a + K \ y = 2^b + T$$
$$ x cdot y = (2^a + K) cdot (2^b + T) = 2^{a+b} + T cdot2^a + K cdot 2^b + K cdot T$$
I chose the power to be $2$ as it would help us with bit-manipulation later on.
Calculating $2^{a+b}$ is easy using bitwise operations as so: $$ 2^{a+b} = 1 << (a+b)$$
But how would we find $a$ and $b$?
We want $2^a$ or $2^b$ to be the largest power of $2$ below our $x$ (or $y$ correspondingly), to take as much ‘volume’ from the original number, and thus making the calculations easier with bit manipulations. So, I just used the $lg$ function, which from what I’ve read it can run in $O(1)$ running-time complexity (Or at worst, $lg lg n$). We have:
$$ a = lfloor lg(x) rfloor, ~~~ b = lfloor lg(y) rfloor$$
We then need to find $K$ which is just the remainder when we subtract $2^a$ from $x$: $$K= x – 2^a = x – (1 << a)$$
However, maybe subtraction isn’t the best idea, maybe it takes too much time, and though about another bit manipulation. All I had to do is to flip the most significant bit (left most bit) which represents the greatest power of $2$ this number consists of, and so I had to pad exactly $a$ $1$‘s and use the $&$ bitwise operation to clear the MSB. We now have a code to find $K$ and $T$ respectively:
$$ K = x~~ &~~ text{int(‘1’ * a, 2)} \ T = y~~ &~~ text{int(‘1’ * b, 2)}$$
Finally, we can add all the factors together, calling the function recursively to compute $K cdot T$ to get:
$$ (1 << (a + b)) + (T << a) + (K << b) + overbrace{text{mult(K,T)}}^{text{recursive call}}$$
def mult(x, y):
if x == 1:
return y
elif y == 1:
return x
elif x == 0 or y == 0:
return 0
base_x = int(log2(x))
base_y = int(log2(y))
K = x & int('1' * base_x, 2)
T = y & int('1' * base_y, 2)
return (1 << (base_x + base_y)) + (T << base_x) + (K << base_y) + mult(K, T)
But oh! from what I’ve tested, this algorithm does not seem to get near the time it takes to multiply two numbers by just using the plain-old $text{*}$ operation, Sob!
times = ()
for _ in range(10000):
x = random.randint(10 ** 900, 10 ** 1000)
y = random.randint(10 ** 900, 10 ** 1000)
start = time.time()
mult(x, y)
end = time.time()
times.append(end - start)
print(sum(times)/len(times))
This tests $1,000$ multiplications of $900$ to $1000$ digits long random integers, then printing the average time. On my machine the average is: 0.01391555905342102 seconds. Python’s regular multiplication won’t even show a number, just 0.0 because it is so fast.
From what I know, Python’s algorithm uses Karatsuba’s algorithm, and it is roughly $O(n^{approx 1.58})$ – I did not analyze mine strictly, but in one sense it runs at approximately: $$O(max (text{#Number_of_on_bits_x, #Number_of_on_bits_y}))$$
Because every recursive call, we turn off the $text{MSB}$ – thus the number of recursive calls we make is the maximum number of bits that are on ($=1$) in $x$ and $y$, which is strictly smaller than the numbers themselves.. thus we can surely say it is $O(max (x,y)) sim O(n)$ as all the other operations in the functions are $O(1)$. So it boils down to the question of ‘why?’ – why is it slower? What have I done wrong in my algorithm so it is slower even that from first glance it seems faster?
Thank you!
