Karatsuba algorithm for multiplying two numbers

Once I had to implement the multiplication of long numbers, through the halves of their representation in C ++. 128-bit numbers were represented as a 64-bit pair. It turned out that to multiply two 128-bit numbers and get all 256 bits of the result by complexity is comparable to 4 pieces of 64-bit halves. How does it work ...



For multiplication, the method of domestic mathematician Anatoly Alekseevich Karatsuba was used . To begin, give your comment on the function. He clarifies a lot by the algorithm.
')

// // Karatsuba multiplication algorithm // // +------+------+ // | x1 | x0 | // +------+------+ // * // +------+------+ // | y1 | y0 | // +------+------+ // = // +-------------+-------------+ // + | x1*y1 | x0*y0 | // +----+-+------+------+------+ // . . . // . . . // +-+------+------+ // + | x0*y1 + x1*y0 | // +-+------+------+ // 

In principle, the commentary makes it clear how the result is calculated.
The essence of the method can be understood if you derive the following formula:



T symbolizes indent, respectively T ^ 2 - double indent.
Here is the code that performs these operations:

 // // Params: // T - is type of x0, x1, y0 and y1 halves // T2 - is type of x, y and half of res // template<typename T, typename T2> inline void Karatsuba_multiply(T * const x, T * const y, T2 * res) { // Define vars (depending from byte order) #define ptrRes ((T*)res) T2 & lowWord = (T2&)(ptrRes[0]); T2 & midWord = (T2&)(ptrRes[1]); T2 & highWord = (T2&)(ptrRes[2]); T & highByte = (T &)(ptrRes[3]); #undef ptrRes const T & x0 = x[0]; const T & x1 = x[1]; const T & y0 = y[0]; const T & y1 = y[1]; // Multiply action lowWord = x0 * y0; highWord = x1 * y1; T2 m1 = x0 * y1; T2 m2 = x1 * y0; midWord += m1; if (midWord < m1) highByte++; midWord += m2; if (midWord < m2) highByte++; } 

Type T is a type of half, type x0, x1, y0, y1.
Type N2 is a type of half of the result, 2 times more than type T.

An example of using the function:

 int main() { typedef unsigned char u8; typedef unsigned short u16; typedef unsigned int u32; u16 a = 1000; u16 b = 2000; u32 r = 0; u8 * a_ptr = (u8*)&a; u8 * b_ptr = (u8*)&b; u16 * r_ptr = (u16*)(void*)&r; Karatsuba_multiply(a_ptr, b_ptr, r_ptr); cout << r; } 

The code with the result of the execution can be found here: codepad.org/1OTeGqhA
The code with the algorithm for different orders of bytes (LE + BE) here: codepad.org/f5Pxtiq1

UPDATE1:
The user mark_ablov noticed that the code lacked #undef.
Corrected code: codepad.org/13U4fuTp
Revised full code (LE + BE): codepad.org/kBazqo8f

UPDATE2:
The user ttim noticed that the above algorithm is not quite the Karatsuba method and indicated the possibility of using 3 multiplications instead of four.

Clarity formula:



Thus, it will be necessary to calculate only 3 products:
1. a * x
2. b * y
3. (Ta + b) * (Tx + y)

Unfortunately, I will not be able to use this formula, since I do not have the ability to multiply the number of bits of 128 bits. Actually, my task was to implement the multiplication of 128-bit numbers. The fact is that the numbers (Ta + b) and (Tx + y) are 128 bits wide ...

UPDATE3:
The susl user continued the discussion of the algorithm and showed that it is not at all necessary to multiply 128-bit numbers.

Another formula:



The function can be rewritten as follows:

 // // Params: // T - is type of x0, x1, y0 and y1 halves // T2 - is type of x, y and half of res // template<typename T, typename T2> inline void Karatsuba_multiply(T * const x, T * const y, T2 * res) { // Define vars (depending from byte order) #define ptrRes ((T*)res) T2 & lowWord = (T2&)(ptrRes[0]); T2 & midWord = (T2&)(ptrRes[1]); T2 & highWord = (T2&)(ptrRes[2]); T & highByte = (T &)(ptrRes[3]); #undef ptrRes const T & x0 = x[0]; const T & x1 = x[1]; const T & y0 = y[0]; const T & y1 = y[1]; // Multiply action lowWord = x0 * y0; highWord = x1 * y1; //T2 m1 = x0 * y1; //T2 m2 = x1 * y0; T2 m = (x0+x1)*(y0+y1) - (lowWord + highWord); //midWord += m1; //if (midWord < m1) highByte++; //midWord += m2; //if (midWord < m2) highByte++; midWord += m; if (midWord < m) highByte++; } 

Corrected example: http://codepad.org/w0INBD77
Corrected example for LE and BE: http://codepad.org/nB9HqWt1