Quelle  mul.txt   Sprache: Text

Multiplication

This describes the multiplication algorithm used by the MPI library.

This is basically a standard "schoolbook" algorithm.  It is slow --
O(mn) for m = #a, n = #b -- but easy to implement and verify.
Basically, we run two nested loops, as illustrated here (R is the
radix):

k = 0
for j <- 0 to (#b - 1)
  for i <- 0 to (#a - 1)
    w = (a[j] * b[i]) + k + c[i+j]
    c[i+j] = w mod R
    k = w div R
  endfor
  c[i+j] = k;
  k = 0;
endfor

It is necessary that 'w' have room for at least two radix R digits.
The product of any two digits in radix R is at most:

 (R - 1)(R - 1) = R^2 - 2R + 1

Since a two-digit radix-R number can hold R^2 - 1 distinct values,
this insures that the product will fit into the two-digit register.

To insure that two digits is enough for w, we must also show that
there is room for the carry-in from the previous multiplication, and
the current value of the product digit that is being recomputed.
Assuming each of these may be as big as R - 1 (and no larger,
certainly), two digits will be enough if and only if:

 (R^2 - 2R + 1) + 2(R - 1) <= R^2 - 1

Solving this equation shows that, indeed, this is the case:

 R^2 - 2R + 1 + 2R - 2 <= R^2 - 1

 R^2 - 1 <= R^2 - 1

This suggests that a good radix would be one more than the largest
value that can be held in half a machine word -- so, for example, as
in this implementation, where we used a radix of 65536 on a machine
with 4-byte words.  Another advantage of a radix of this sort is that
binary-level operations are easy on numbers in this representation.

Here's an example multiplication worked out longhand in radix-10,
using the above algorithm:

   a =     999
   b =   x 999
  -------------
   p =   98001

w = (a[jx] * b[ix]) + kin + c[ix + jx]
c[ix+jx] = w % RADIX
k = w / RADIX
                                                               product
ix jx a[jx] b[ix] kin w c[i+j] kout 000000
0 0 9 9 0 81+0+0 1 8 000001
0 1 9 9 8 81+8+0 9 8 000091
0 2 9 9 8 81+8+0 9 8 000991
    8   0 008991
1 0 9 9 0 81+0+9 0 9 008901
1 1 9 9 9 81+9+9 9 9 008901
1 2 9 9 9 81+9+8 8 9 008901
    9   0 098901
2 0 9 9 0 81+0+9 0 9 098001
2 1 9 9 9 81+9+8 8 9 098001
2 2 9 9 9 81+9+9 9 9 098001

------------------------------------------------------------------
 This Source Code Form is subject to the terms of the Mozilla Public
 # License, v. 2.0. If a copy of the MPL was not distributed with this
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.