Series: Multiple-Precision Arithmetic

Basic Multiplication

Published on , updated on

After addressing multiple-precision addition and subtraction, we now turn to multiplication of two multiple-precision numbers. Once again, we use the number representation and notation introduced earlier.

Several algorithms exist for doing multiple-precision multiplication. This post will present the basic, pencil-and-paper-like method. Basically, it consists of two parts: Multiplying a number by a single digit and adding together the sub-results, aligned appropriately.

Multiple digits times a single digit

Let us first consider multiplying an nn-digit multiple-precision integer by a single digit. More precisely, we wish to compute z←αv+y+k0z \leftarrow \alpha v + y + k_0 where

v=(vn−1…v1v0)b,1≤α≤b−1,y=(yn−1…y1y0)b,0≤k0≤b−1.v = (v_{n-1} \ldots v_1 v_0)_b, \quad 1 \leq \alpha \leq b-1, \quad y = (y_{n-1} \ldots y_1 y_0)_b, \quad 0 \leq k_0 \leq b-1.

Using this information it is straightforward to show that bn−1≤z≤bn+1−1b^{n-1} \leq z \leq b^{n+1}-1, which implies that the result will fit into nn or n+1n+1 digits, so we set z=(zn…z1z0)bz = (z_n \ldots z_1 z_0)_b, where znz_n may be zero.

We now have the algorithm:

zi←(αvi+yi+ki)  mod  b,ki+1←⌊αvi+yi+kib⌋,\begin{aligned} z_i &\leftarrow (\alpha v_i + y_i + k_i) \;\text{mod}\; b, \\ k_{i+1} &\leftarrow \left\lfloor \frac{\alpha v_i + y_i + k_i}{b} \right\rfloor, \end{aligned}

for i=0,1,…,n−1i = 0, 1, \ldots, n-1 and finally setting zn=knz_n = k_n.

To realize that the algorithm computes what it is supposed to, observe first that

αvi+yi+ki=(αvi+yi+ki)  mod  b+⌊(αvi+yi+ki)/b⌋b=zi+ki+1b.\alpha v_i + y_i + k_i = (\alpha v_i + y_i + k_i) \;\text{mod}\; b + \left\lfloor (\alpha v_i + y_i + k_i)/b \right\rfloor b = z_i + k_{i+1} b.

Then we have

αv+y=∑i=0n−1(αvi+yi)bi=∑i=0n−1(αvi+yi+ki)bi−∑i=0n−1kibi=∑i=0n−1(zi+ki+1b)bi−∑i=0n−1kibi=∑i=0n−1zibi+∑i=0n−1ki+ibi+1−∑i=0n−1kibi=∑i=0n−1zibi+znbn−k0=z−k0,\begin{aligned} \alpha v + y &= \sum_{i=0}^{n-1} (\alpha v_i + y_i) b^i = \sum_{i=0}^{n-1} (\alpha v_i + y_i + k_i) b^i - \sum_{i=0}^{n-1} k_i b^i \\ &= \sum_{i=0}^{n-1} (z_i + k_{i+1} b) b^i - \sum_{i=0}^{n-1} k_i b^i = \sum_{i=0}^{n-1} z_i b^i + \sum_{i=0}^{n-1} k_{i+i} b^{i+1} - \sum_{i=0}^{n-1} k_i b^i \\ &= \sum_{i=0}^{n-1} z_i b^i + z_n b^n - k_0 = z - k_0, \end{aligned}

which is what we wanted.

From the expression for ziz_i we see that they naturally fit into a digit, that is, 0≤zi≤b−10 \leq z_i \leq b-1 for i=0,1,…,n−1i = 0, 1, \ldots, n-1. But what about the kik_i's and thus also znz_n? Assume that 0≤ki≤b−10 \leq k_i \leq b-1. We then have αvi+yi+ki≤(b−1)+(b−1)(b−1)+(b−1)=b2−1,\alpha v_i + y_i + k_i \leq (b-1) + (b-1)(b-1) + (b-1) = b^2-1, which implies that ki+1=⌊(αvi+yi+ki)/b⌋≤b−1k_{i+1} = \lfloor (\alpha v_i + y_i + k_i)/b \rfloor \leq b-1. Since we have also assumed 0≤k0≤b−10 \leq k_0 \leq b-1 we have, by induction, that 0≤zi≤b−10 \leq z_i \leq b-1 for i=0,1,…,ni = 0, 1, \ldots, n.

Multiple digits times multiple digits

We now turn to multiplying two multiple-precision numbers. More specifically, we wish to multiply

u=(um−1…u1u0)b,andv=(vn−1…v1v0)b,u = (u_{m-1} \ldots u_1 u_0)_b, \quad \text{and} \quad v = (v_{n-1} \ldots v_1 v_0)_b,

which implies bm+n−2≤uv<bm+nb^{m+n-2} \leq u v < b^{m+n}. So we set w=(wm+n−1…w1w0)bw = (w_{m+n-1} \ldots w_1 w_0)_b and aim to compute

w←uv=∑i=0m−1uivbi.w \leftarrow u v = \sum_{i=0}^{m-1} u_i v b^i.

We observe two things:

  1. Multiplication by bib^i corresponds to shifting left the number of digits given by ii.
  2. The product uivu_i v is of type single-digit times multiple-digit and can be computed using the algorithm from the previous section.

Putting these pieces together we get:

(wn…w1w0)b←u0v,(wn+1…w2w1)b←(wn…w2w1)+u1v,⋮(wn+m−1…wmwm−1)b←(wn+m−2…wmwm−1)b+um−1v.\begin{aligned} (w_n \ldots w_1 w_0)_b &\leftarrow u_0 v, \\ (w_{n+1} \ldots w_2 w_1)_b &\leftarrow (w_n \ldots w_2 w_1) + u_1 v, \\ &\vdots \\ (w_{n+m-1} \ldots w_m w_{m-1})_b &\leftarrow (w_{n+m-2} \ldots w_m w_{m-1})_b + u_{m-1} v. \end{aligned}

The algorithm can actually be generalized slightly if we compute

w←(wn−1…w1w0)b+u  vw \leftarrow (w_{n-1} \ldots w_1 w_0)_b + u \; v

instead. All we need to do is replace the first step in the algorithm by

(wn…w1w0)b←(wn−1…w1w0)b+u0v.(w_n \ldots w_1 w_0)_b \leftarrow (w_{n-1} \ldots w_1 w_0)_b + u_0 v.

We now move on to multiple-precision division by first adressing the relatively simple case of short division.