You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
The multiplication is possibly the most important of the 4 basic arithmetic operators in a big integer library. This is because it is the one most likely to work on and generate truly big numbers.
Development-wise, it is an interesting operation because of its time complexity.
The case to consider is the multiplication of 2 numbers $a$ and $b$ of roughly the same size, i.e. each made of a certain number $l$ of limbs. If e.g. $a$ was made of only 1 limb the complexity for the multiplication would simply be $\text{O}(l)$ (each limb of $b$ must be multiplied by $a$). Being limited to such a case would not make things interesting.
In practice, operands rarely are of the exact same size but algorithms can easily cope with any difference and the theoretical considerations can be solved by filling the smaller number with leading zeroes.
With $p = a \times b$, we note $l = \max(\text{size}(a), \text{size}(b))$.
Long multiplication:
This is how we learn multiplications in school, by taking each pair of digits from both operands, and multiply them together, with a carry.
The computational complexity is $\text{Θ}(l^2)$.
The other algorithms (addition, subtraction, ...) having a $\text{O}(l)$ complexity, they are practically free in comparison.
The demonstration that better algorithms exist has triggered a decades-long effort by mathematicians to improve performance further and further. This resulted in several algorithms, each better than the previous and suitable for increasingly big numbers.
Karatsuba:
Each time $l$ doubles, the Karatsuba algorithm trades 1 multiplication for several additions/subtractions, thus making it 3 multiplications instead of 4. The idea is that when $l$ is big enough, the cost of the additional operations is less than what the time saved by skipping the fourth multiplication. As a result, it is not suitable for multiplying (relatively) small numbers.
The complexity is $\text{Θ}(l^{\log_2(3)}) = \text{O}(l^{1.585})$.
Toom-Cook algorithms:
The Toom-Cook algorithms form a family of divide-and-conquer algorithms similar to Karatsuba, except for the fact they are not limited to splitting operands in halves. The smaller the parts, the more multiplications it can trade against simpler operations. For this reason, each Toom-Cook algorithm needs longer and longer numbers to make the trade worth it.
The complexity for Toom-3 is $\text{Θ}(l^{\log_3(5)}) = \text{O}(l^{1.465})$.
Schönhage–Strassen:
This method is not implemented yet.
This is yet another algorithm more efficient that the ones above, but this time, applicable for large numbers only, with millions of digits.
Multiplication function
The multiplication function compares the size of the operands passed to it with hardcoded threshold, which are the best guess as what values give the shortest calculation time:
If $\text{size}(a) > \text{threshold}_\text{Schönhage–Strassen}$
→ return the product calculated using the Schönhage–Strassen algorithm.
If $\text{size}(a) > \text{threshold}_\text{Toom-3}$
→ return the product calculated using the Toom-3 algorithm.
If $\text{size}(a) > \text{threshold}_\text{Karatsuba}$
→ return the product calculated using the Karatsuba algorithm.
Otherwise
→ return the product calculated using the long multiplication algorithm.
This is the same approach as every other library that implement multiplication through several algorithms, the only differences being what algorithms were implemented, small details inside them and the thresholds to decide which one to run.
Long multiplication
The long multiplication algorithm works on the principle that:
The implementation is straight-forward. Indexing limbs of operands starting at index 0:
Initialize the variable $p$ to store the product, so that $2 l = \text{size}(p) = \text{size}(a) + \text{size}(b)$, with all its limbs initialized as $0$.
Start a for-loop: $i$ from $0$ to $\text{size}(a) - 1$
Initialize a variable to store the carry modulo $2^{32}$: $c \leftarrow 0$.
Start a nested for-loop: $j$ from $0$ to $\text{size}(a) - 1$.
Do $p_{i+j} \leftarrow \big(p_{i+j} + a_i \times b_j + c\big) \bmod 2^{32}$ and $c \leftarrow \lfloor {\big(p_{i+j} + a_i \times b_j + c\big)} / {2^{32}} \rfloor$.
In practice, these 2 operations are done by doing $\big(\text{result}_{i+j} + a_i \times b_j + c\big)$ into a 64-bit integer1, then splitting the integer into $\text{result}$ and $c$.
At the end of the outer loop, a carry may remain. Do $p_{2l-1} \leftarrow c$.
Trim $p$ if the carry was 0.
Karatsuba algorithm
Karatsuba, from the name of its inventor, is a divide-and-conquer algorithm. Its principle run the multiplication as:
For convenience, the 4 sums, i.e. substrings of the integers' limbs, appearing in the above formula will be noted with Greek letters: $\alpha_0, \alpha_1, \beta_0, \beta_1$ (we keep letters from the Latin alphabet to represent single limbs of a number).
Although it looks like this just made the formula longer and more complicated, hence longer to calculate, the key to the algorithm is about being able to reuse $\alpha_0 \times \beta_0$ and $2^{64\lceil l/2 \rceil} \times \alpha_1 \times \beta_1$ where necessary.
By doing so, we allow the Karatsuba algorithm to execute in 3 multiplications what the long multiplication would do in 4.
This 1 operation being optimized away at every step is a big deal: if we were to assume that at a given size $s$, the cost of the long multiplication is 1 (of an arbitrary unit) vs. 5 for the Karatsuba algorithm (including the additional operations, allocation of memory, ...), the table below shows how it evolves for each doubling of $s$.
Size
Long multiplication
Cost of Karatsuba
Ratio
1s
1
5
1:5
2s
4
15
1:3.75
4s
16
45
1:2.81
8s
64
135
1:2.11
16s
256
405
1:1.58
32s
1,024
1,215
1:1.19
64s
4,096
3,645
1.12:1
128s
16,384
10,935
1.50:1
256s
65,536
32,805
2.00:1
512s
262,144
98,415
2.66:1
1024s
1,048,576
295,245
3.55:1
The Karatsuba algorithm consists in:
Initialize the variable $p$ to store the product, so that $2l = \text{size}(p) = \text{size}(a) + \text{size}(b)$, with all its limbs initialized as $0$.
Split $a$ a,d $b$ at the limb indexed $k = \lceil l /2 \rceil$.
This allows to obtain the subparts of $a$, $b$ and $p$ as:
$\alpha_0$ and $\alpha_1$,
$\beta_0$ and $\beta_1$,
$\psi_0$, $\psi_1$ and $\psi_2$.
Note that $\psi_0$ and $\psi_2$ do not overlap in $p$. However, $\psi_1$ does overlap with the other two in memory.
Do $\psi_0 \leftarrow \alpha_0 \times \beta_0$, by a recursive call to the multiplication algorithm.
Do $\psi_2 \leftarrow \alpha_1 \times \beta_1$ the same way.
Initialiée a variable $p_\text{mid}$, of size $l$.
Finally, do $\psi_1 \leftarrow \psi_1 + p_\text{mid}$.
Toom-Cook
Toom-Cook is a divide-and-conquer algorithm (rather, an infinite family of algorithms) named after its two discoverers: Andrei Toom who constructed a method for multiplication with and Stephen Cook who turned it in an algorithm. The method produces an infinite family of algorithms.
Just like Karatsuba, Toom-Cook swaps some of the required multiplications for a finite number of simpler operations, all of which are $\text{O}(l)$ in complexity.
The more the operands are broken down (Toom-2 is equivalent to Karatsuba, Toom-3 breaks them down into 3, Toom-4 breaks them down into 4, etc), the better the complexity. However, this also comes at the cost of having more and more of the simpler operations. Hence, doing so requires increasingly bigger numbers to make the trade worth it.
The complexity of Toom-Cook, when breaking down operands of length $l$ into $k$ parts of length $l/k$, is $\text{Θ}(l^{\log_k(2k-1)})$.
Toom-k
Complexity
Approximate
2
$\text{Θ}(l^{\log_2(3)})$
$\text{O}(l^{1.585})$
4
$\text{Θ}(l^{\log_3(5)})$
$\text{O}(l^{1.465})$
5
$\text{Θ}(l^{\log_4(7)})$
$\text{O}(l^{1.404})$
6
$\text{Θ}(l^{\log_5(9)})$
$\text{O}(l^{1.366})$
...
...
...
10
$\text{Θ}(l^{\log_{10}(19)})$
$\text{O}(l^{1.279})$
20
$\text{Θ}(l^{\log_{20}(39)})$
$\text{O}(l^{1.223})$
50
$\text{Θ}(l^{\log_{50}(99)})$
$\text{O}(l^{1.175})$
100
$\text{Θ}(l^{\log_{100}(199)})$
$\text{O}(l^{1.150})$
Because of the diminishing returns, Toom-Cook algorithms are not implemented for large values of $k$, generally not above Toom-4 or Toom-5.
For the sake of comparison, when the size of the operands is multiplied by 4, it generates:
for the long division, 16 times as many multiplications.
for Karatsuba, 9 times as many multiplications.
for Toom-4, 7 times as many multiplications (not counting the multiplications by small constants due to their $\text{O}(l)$ complexity as discussed earlier).
Principle of Toom-Cook
To understand Toom-Cook algorithms, we have to express big integers as the result of polynomials. Since big integers are encoded as limbs, it looks natural to express
$\displaystyle a = \sum_{i \geq 0} \limits a_i \times 2^{32i}$
as the result of a polynomial: $a = A(2^{32})$
where $\displaystyle A(x) = \sum_{i \geq 0} \limits a_i \times x^i$.
This is not a polynomial expression we can use here, though. Toom-k splits operands in $k$ parts so we need to express them as a polynomial with $k$ coefficients, i.e. of degree $k-1$: $\displaystyle A(x) = \large \sum_{i=0}^{k-1} \bigg( \normalsize \sum_{j=0}^{\lceil l/k \rceil-1} a_{j + i \times \lceil l/k \rceil} \times 2^{32j} \bigg) \times x^i$, that we note $\displaystyle A(x) = \sum_{i=0}^{k-1} \alpha_i \times x^i$
The product $p = a \times b$ is now replaced by the polynomial $P$ of degree $2(k-1)$.
Obviously, we are not doing it to just calculate the direct product of $A(x)$ and $B(x)$, each coefficient at a time.
Polynomial interpolation
Thankfully, polynomial interpolation allows to calculate a polynomial of degree $2(k-1)$ using only $2k-1$ points.
We will do exactly that with easy-to-calculate points.
The chosen points are 0 and numbers in the form of $\pm 2^i, i \geq 0$ since they are the easiest to deal with for a computer.
For Toom-3, we need 5 points, that we actually calculate:
We have now replaced the product of 2 operands of length $l$ by 5 products, each with operands of length $l/3$ (hence the complexity).
Linear equations
Remembering that $P(x) = \psi_4 \times x^4 + \psi_3 \times x^3 + \psi_2 \times x^2 + \psi_1 \times x + \psi_0$, we have a system of 5 linear equations:
We can solve it and deduce the values $\psi_i$ but was it the best we could do?
Well... no and while it won't improve the complexity, we can slightly simplify the above equations.
The clue is to note that when $x \rightarrow \infty$, $P(x) \displaystyle \equiv \psi_4 \times x^4$ and $A(x) \times B(x) \equiv (\alpha_2 \times x^2) \times (\beta_2 \times x^2) = \alpha_2 \times \beta_2 \times x^4$.
The asymptotic behavior of $P$ results in a new "point", usually noted $P(\infty)$, defined as $P(\infty) = \psi_4 = \alpha_2 \times \beta_2$.
This is much simpler than the above equations and gives $\psi_4$ directly. As such, it avoids some calculations.
Now, we can replace one of the points chosen above by $P(\infty)$. For instance, if $P(\text{-}2)$ is discarded.
Note that if you were to finish solving the equation system (i.e. get rid of the $\psi_i$ on the right of the equal signs), you would end up with equalities that could be expressed in the form of a product of matrices.
This matrix product is how the solution to Toom-3 is presented on Wikipedia but it is not completely optimal in terms of number of operations to execute.
As big integers are coded with limbs in base $2^{32}$, the actual algorithm works in a base of the form ${2^{32}}^i$. Let alone numbers big enough for the algorithm to be efficient (they would need thousands of digits), even i = 1 would make the numbers too big for an illustration; we are talking about multiplying 2 numbers whose scale is in the octillions (US/short form) or thousand quadrillions (European/long form), with a resulting product in the scale of Septendecillion (US) or Nonillion (European), i.e. over 50 base-10 digits.
For clarity's sake, the illustration is done with numbers under 109, in a base 10k (convenient for human beings), then in a base 2k (convenient for computers).
Applying the formulas we obtained from the linear equations, we deduce that $P(x) = 121{,}401 x^4 + 530{,}514 x^3 + 1{,}116{,}450 x^2 + 662{,}382 x + 253{,}269$ and therefore, thanks to a sum of the shifted coefficients $\psi_i$, equivalent to multiplying them by ${10^3}^i$:
Applying the formulas we obtained from the linear equations, we deduce that $P(x) = 110{,}097 x^4 + 818{,}329 x^3 + 977{,}476 x^2 + 389{,}029 x + 49{,}029$.
The sum of the coefficients, after shifting each by a multiple of 10 bits gives the program: $P(2^{10}) = 121{,}932{,}631{,}112{,}635{,}269$, the very same result as the calculation made in base 103.
Footnotes
As determined in the basics, a 64-bit integer is enough to store the result of any 32-bit combination of the form: a x b + c + d. ↩