janmr blog

Basic Multiple-Precision Multiplication

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=(vn1v1v0)b,1αb1,y=(yn1y1y0)b,0k0b1.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 bn1zbn+11b^{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=(znz1z0)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,,n1i = 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)/bb=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=0n1(αvi+yi)bi=i=0n1(αvi+yi+ki)bii=0n1kibi=i=0n1(zi+ki+1b)bii=0n1kibi=i=0n1zibi+i=0n1ki+ibi+1i=0n1kibi=i=0n1zibi+znbnk0=zk0,\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, 0zib10 \leq z_i \leq b-1 for i=0,1,,n1i = 0, 1, \ldots, n-1. But what about the kik_i's and thus also znz_n? Assume that 0kib10 \leq k_i \leq b-1. We then have αvi+yi+ki(b1)+(b1)(b1)+(b1)=b21,\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)/bb1k_{i+1} = \lfloor (\alpha v_i + y_i + k_i)/b \rfloor \leq b-1. Since we have also assumed 0k0b10 \leq k_0 \leq b-1 we have, by induction, that 0zib10 \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=(um1u1u0)b,andv=(vn1v1v0)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+n2uv<bm+nb^{m+n-2} \leq u v < b^{m+n}. So we set w=(wm+n1w1w0)bw = (w_{m+n-1} \ldots w_1 w_0)_b and aim to compute

wuv=i=0m1uivbi.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:

(wnw1w0)bu0v,(wn+1w2w1)b(wnw2w1)+u1v,(wn+m1wmwm1)b(wn+m2wmwm1)b+um1v.\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(wn1w1w0)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

(wnw1w0)b(wn1w1w0)b+u0v.(w_n \ldots w_1 w_0)_b \leftarrow (w_{n-1} \ldots w_1 w_0)_b + u_0 v.