Thu, 22 Jan 2015 13:21:57 +0100
Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6
michael@0 | 1 | Multiplication |
michael@0 | 2 | |
michael@0 | 3 | This describes the multiplication algorithm used by the MPI library. |
michael@0 | 4 | |
michael@0 | 5 | This is basically a standard "schoolbook" algorithm. It is slow -- |
michael@0 | 6 | O(mn) for m = #a, n = #b -- but easy to implement and verify. |
michael@0 | 7 | Basically, we run two nested loops, as illustrated here (R is the |
michael@0 | 8 | radix): |
michael@0 | 9 | |
michael@0 | 10 | k = 0 |
michael@0 | 11 | for j <- 0 to (#b - 1) |
michael@0 | 12 | for i <- 0 to (#a - 1) |
michael@0 | 13 | w = (a[j] * b[i]) + k + c[i+j] |
michael@0 | 14 | c[i+j] = w mod R |
michael@0 | 15 | k = w div R |
michael@0 | 16 | endfor |
michael@0 | 17 | c[i+j] = k; |
michael@0 | 18 | k = 0; |
michael@0 | 19 | endfor |
michael@0 | 20 | |
michael@0 | 21 | It is necessary that 'w' have room for at least two radix R digits. |
michael@0 | 22 | The product of any two digits in radix R is at most: |
michael@0 | 23 | |
michael@0 | 24 | (R - 1)(R - 1) = R^2 - 2R + 1 |
michael@0 | 25 | |
michael@0 | 26 | Since a two-digit radix-R number can hold R^2 - 1 distinct values, |
michael@0 | 27 | this insures that the product will fit into the two-digit register. |
michael@0 | 28 | |
michael@0 | 29 | To insure that two digits is enough for w, we must also show that |
michael@0 | 30 | there is room for the carry-in from the previous multiplication, and |
michael@0 | 31 | the current value of the product digit that is being recomputed. |
michael@0 | 32 | Assuming each of these may be as big as R - 1 (and no larger, |
michael@0 | 33 | certainly), two digits will be enough if and only if: |
michael@0 | 34 | |
michael@0 | 35 | (R^2 - 2R + 1) + 2(R - 1) <= R^2 - 1 |
michael@0 | 36 | |
michael@0 | 37 | Solving this equation shows that, indeed, this is the case: |
michael@0 | 38 | |
michael@0 | 39 | R^2 - 2R + 1 + 2R - 2 <= R^2 - 1 |
michael@0 | 40 | |
michael@0 | 41 | R^2 - 1 <= R^2 - 1 |
michael@0 | 42 | |
michael@0 | 43 | This suggests that a good radix would be one more than the largest |
michael@0 | 44 | value that can be held in half a machine word -- so, for example, as |
michael@0 | 45 | in this implementation, where we used a radix of 65536 on a machine |
michael@0 | 46 | with 4-byte words. Another advantage of a radix of this sort is that |
michael@0 | 47 | binary-level operations are easy on numbers in this representation. |
michael@0 | 48 | |
michael@0 | 49 | Here's an example multiplication worked out longhand in radix-10, |
michael@0 | 50 | using the above algorithm: |
michael@0 | 51 | |
michael@0 | 52 | a = 999 |
michael@0 | 53 | b = x 999 |
michael@0 | 54 | ------------- |
michael@0 | 55 | p = 98001 |
michael@0 | 56 | |
michael@0 | 57 | w = (a[jx] * b[ix]) + kin + c[ix + jx] |
michael@0 | 58 | c[ix+jx] = w % RADIX |
michael@0 | 59 | k = w / RADIX |
michael@0 | 60 | product |
michael@0 | 61 | ix jx a[jx] b[ix] kin w c[i+j] kout 000000 |
michael@0 | 62 | 0 0 9 9 0 81+0+0 1 8 000001 |
michael@0 | 63 | 0 1 9 9 8 81+8+0 9 8 000091 |
michael@0 | 64 | 0 2 9 9 8 81+8+0 9 8 000991 |
michael@0 | 65 | 8 0 008991 |
michael@0 | 66 | 1 0 9 9 0 81+0+9 0 9 008901 |
michael@0 | 67 | 1 1 9 9 9 81+9+9 9 9 008901 |
michael@0 | 68 | 1 2 9 9 9 81+9+8 8 9 008901 |
michael@0 | 69 | 9 0 098901 |
michael@0 | 70 | 2 0 9 9 0 81+0+9 0 9 098001 |
michael@0 | 71 | 2 1 9 9 9 81+9+8 8 9 098001 |
michael@0 | 72 | 2 2 9 9 9 81+9+9 9 9 098001 |
michael@0 | 73 | |
michael@0 | 74 | ------------------------------------------------------------------ |
michael@0 | 75 | This Source Code Form is subject to the terms of the Mozilla Public |
michael@0 | 76 | # License, v. 2.0. If a copy of the MPL was not distributed with this |
michael@0 | 77 | # file, You can obtain one at http://mozilla.org/MPL/2.0/. |