fn karatsuba(acc: &mut [u64], x: &[u64], y: &[u64])Expand description
Karatsuba multiplication:
The idea is that we break x and y up into two smaller numbers that each have about half as many digits, like so (note that multiplying by b is just a shift):
x = x0 + x1 * b y = y0 + y1 * b
With some algebra, we can compute x * y with three smaller products, where the inputs to each of the smaller products have only about half as many digits as x and y:
x * y = (x0 + x1 * b) * (y0 + y1 * b)
x * y = x0 * y0 + x0 * y1 * b + x1 * y0 * b + x1 * y1 * b^2
Let p0 = x0 * y0 and p2 = x1 * y1:
x * y = p0 + (x0 * y1 + x1 * y0) * b + p2 * b^2
The real trick is that middle term:
x0 * y1 + x1 * y0 = x0 * y1 + x1 * y0 - p0 + p0 - p2 + p2 = x0 * y1 + x1 * y0 - x0 * y0 - x1 * y1 + p0 + p2
Now we complete the square:
= -(x0 * y0 - x0 * y1 - x1 * y0 + x1 * y1) + p0 + p2 = -((x1 - x0) * (y1 - y0)) + p0 + p2
Let p1 = (x1 - x0) * (y1 - y0), and substitute back into our original formula:
x * y = p0 + (p0 + p2 - p1) * b + p2 * b^2
Where the three intermediate products are:
p0 = x0 * y0 p1 = (x1 - x0) * (y1 - y0) p2 = x1 * y1
In doing the computation, we take great care to avoid unnecessary temporary variables (since creating a BigUint requires a heap allocation): thus, we rearrange the formula a bit so we can use the same temporary variable for all the intermediate products:
x * y = p2 * b^2 + p2 * b + p0 * b + p0 - p1 * b
The other trick we use is instead of doing explicit shifts, we slice acc at the appropriate offset when doing the add.