michael@0: Multiplication
michael@0: 
michael@0: This describes the multiplication algorithm used by the MPI library.
michael@0: 
michael@0: This is basically a standard "schoolbook" algorithm.  It is slow --
michael@0: O(mn) for m = #a, n = #b -- but easy to implement and verify.
michael@0: Basically, we run two nested loops, as illustrated here (R is the
michael@0: radix):
michael@0: 
michael@0: k = 0
michael@0: for j <- 0 to (#b - 1)
michael@0:   for i <- 0 to (#a - 1)
michael@0:     w = (a[j] * b[i]) + k + c[i+j]
michael@0:     c[i+j] = w mod R
michael@0:     k = w div R
michael@0:   endfor
michael@0:   c[i+j] = k;
michael@0:   k = 0;
michael@0: endfor
michael@0: 
michael@0: It is necessary that 'w' have room for at least two radix R digits.
michael@0: The product of any two digits in radix R is at most:
michael@0: 
michael@0: 	(R - 1)(R - 1) = R^2 - 2R + 1
michael@0: 
michael@0: Since a two-digit radix-R number can hold R^2 - 1 distinct values,
michael@0: this insures that the product will fit into the two-digit register.
michael@0: 
michael@0: To insure that two digits is enough for w, we must also show that
michael@0: there is room for the carry-in from the previous multiplication, and
michael@0: the current value of the product digit that is being recomputed.
michael@0: Assuming each of these may be as big as R - 1 (and no larger,
michael@0: certainly), two digits will be enough if and only if:
michael@0: 
michael@0: 	(R^2 - 2R + 1) + 2(R - 1) <= R^2 - 1
michael@0: 
michael@0: Solving this equation shows that, indeed, this is the case:
michael@0: 
michael@0: 	R^2 - 2R + 1 + 2R - 2 <= R^2 - 1
michael@0: 
michael@0: 	R^2 - 1 <= R^2 - 1
michael@0: 
michael@0: This suggests that a good radix would be one more than the largest
michael@0: value that can be held in half a machine word -- so, for example, as
michael@0: in this implementation, where we used a radix of 65536 on a machine
michael@0: with 4-byte words.  Another advantage of a radix of this sort is that
michael@0: binary-level operations are easy on numbers in this representation.
michael@0: 
michael@0: Here's an example multiplication worked out longhand in radix-10,
michael@0: using the above algorithm:
michael@0: 
michael@0:    a =     999
michael@0:    b =   x 999
michael@0:   -------------
michael@0:    p =   98001
michael@0: 
michael@0: w = (a[jx] * b[ix]) + kin + c[ix + jx]
michael@0: c[ix+jx] = w % RADIX
michael@0: k = w / RADIX
michael@0:                                                                product
michael@0: ix	jx	a[jx]	b[ix]	kin	w	c[i+j]	kout	000000
michael@0: 0	0	9	9	0	81+0+0	1	8	000001
michael@0: 0	1	9	9	8	81+8+0	9	8	000091
michael@0: 0	2	9	9	8	81+8+0	9	8	000991
michael@0: 				8			0	008991
michael@0: 1	0	9	9	0	81+0+9	0	9	008901
michael@0: 1	1	9	9	9	81+9+9	9	9	008901
michael@0: 1	2	9	9	9	81+9+8	8	9	008901
michael@0: 				9			0	098901
michael@0: 2	0	9	9	0	81+0+9	0	9	098001
michael@0: 2	1	9	9	9	81+9+8	8	9	098001
michael@0: 2	2	9	9	9	81+9+9	9	9	098001
michael@0: 
michael@0: ------------------------------------------------------------------
michael@0:  This Source Code Form is subject to the terms of the Mozilla Public
michael@0:  # License, v. 2.0. If a copy of the MPL was not distributed with this
michael@0:  # file, You can obtain one at http://mozilla.org/MPL/2.0/.