karatsuba.a68
1 COMMENT
2
3 @section Synopsis
4
5 KARATSUBA method, a divide-and-conquer multiprecision multiplication algorithm.
6
7 This multiplication algorithm is described for example in Knuth section 4.3.3 part A.
8
9 KARATSUBA computes the product of two long mubers x and y, by splitting-in-half and recursion.
10 Let R be the radix where the halfway split occurs then x = Rx1 + x0 and y = Ry1 + y0, and
11
12 xy = R (R + 1) x1y1 - R (x1-x0)(y1-y0) + (R + 1) x0y0
13
14 This asks three multiplications of N/2 limbs, instead of four such multiplications as in
15 grammar school multiplication. Recursive splitting reduces complexity from O(N^2) to Q(N^~1.6).
16 KARATSUBA's overhead outweighs its gains for smaller values of N. The threshold is platform dependent.
17
18 Note that this scheme can be further optimized when x = y.
19
20 In this example, long numbers (NUMs) are strings of integers. Each integer is a digit in base 'radix'.
21 The least significant digit (with weight 1) is the rightmost one, which is natural.
22
23 COMMENT
24
25
26 OP = = (NUM x, y) BOOL:
27 IF INT len = UPB x; len /= UPB y
28 THEN FALSE
29 ELSE BOOL ok := TRUE;
30 FOR i TO len WHILE ok
31 DO ok := x[i] = y[i]
32 OD;
33 ok
34 FI;
35
36
37 OP /= = (NUM x, y) BOOL:
38 NOT (x = y);
39
40
41 PROC nil = (REF NUM z) VOID:
42 BEGIN INT len = UPB z;
43 FOR i TO len
44 DO z[i] := 0
45 OD
46 END;
47
48
49 PROC show = (STRING txt, NUM x) VOID: printf (($"<"g">"$, txt, $-3d$, x, $l$));
50
51
52 PROC mul = (REF NUM z, NUM x, y) VOID:
53 IF INT len = UPB x;
54 len = 1
55 THEN # Trivial case #
56 INT p := x[1] * y[1];
57 z[1] := p OVER radix; z[2] := p MOD radix
58 ELSE # Grammar school multiplication #
59 nil (z);
60 FOR j FROM UPB x DOWNTO 1
61 DO INT carry := 0;
62 FOR i FROM UPB y DOWNTO 1
63 DO INT p = z[i + j] + y[j] * x[i] + carry;
64 carry := p OVER radix; z[i + j] := p MOD radix
65 OD;
66 z[j] +:= carry
67 OD
68 FI;
69
70
71 PROC add = (REF NUM z, NUM x, INT pos) BOOL:
72 BEGIN INT i := pos, INT xlen = UPB x;
73
74 # Add #
75 FOR j FROM xlen DOWNTO 1 WHILE i > 0
76 DO z[i] +:= x[j]; i -:= 1
77 OD;
78
79 # Normalize the sum #
80 INT carry := 0;
81 FOR j FROM pos DOWNTO 1
82 DO z[j] +:= carry; carry := z[j] OVER radix; z[j] %*:= radix
83 OD;
84
85 carry > 0
86 END;
87
88
89 PROC sub = (REF NUM z, NUM x, INT pos) BOOL:
90 BEGIN INT i := pos, INT zlen = UPB z, xlen = UPB x;
91
92 # Subtract #
93 FOR j FROM xlen DOWNTO 1 WHILE i > 0
94 DO z[i] -:= x[j]; i -:= 1
95 OD;
96
97 # Correct negative difference #
98 INT sign := 0;
99 FOR j FROM 1 TO zlen
100 DO IF z[j] /= 0
101 THEN IF sign = 0
102 THEN sign := (z[j] >= 0 | 1 | -1)
103 FI;
104 (sign = -1 | z[j] := -z[j])
105 FI
106 OD;
107
108 # Normalize the difference #
109 INT carry := 0;
110 FOR j FROM pos DOWNTO 1
111 DO z[j] -:= carry;
112 IF z[j] < 0
113 THEN z[j] +:= radix; carry := 1
114 ELSE carry := 0
115 FI
116 OD;
117
118 sign = -1
119 END;
120
121
122 PROC karat = (REF NUM z, NUM x, y) VOID:
123 IF INT len = UPB x;
124 ODD len
125 THEN mul (z, x, y)
126 ELSE COMMENT
127 Numbers x and y are split into two parts of equal length, like so
128
129 High Low
130 +----------+----------+
131 | x1 | x0 |
132 +----------+----------+
133
134 +----------+----------+
135 | y1 | y0 |
136 +----------+----------+
137
138 COMMENT
139
140 INT half = len OVER 2;
141 [half] INT x1, x0, y1, y0;
142 FOR i TO len
143 DO IF i > half
144 THEN x0[i - half] := x[i]; y0[i - half] := y[i]
145 ELSE x1[i] := x[i]; y1[i] := y[i]
146 FI
147 OD;
148
149 # Compute products #
150 [len] INT x1y1, x0y0, xmym;
151 karat (x0y0, x0, y0); karat (x1y1, x1, y1);
152 BOOL p = sub (x1, x0, half), q = sub (y1, y0, half);
153 karat (xmym, x1, y1);
154
155 COMMENT
156 Factors R(R + 1) etc represent the positions where the three products must be added.
157
158 High Low
159 +--------+--------+ +--------+--------+
160 | x1*y1 | | x0*y0 |
161 +--------+--------+ +--------+--------+
162 +--------+--------+
163 + | x1*y1 |
164 +--------+--------+
165 +--------+--------+
166 + | x0*y0 |
167 +--------+--------+
168 +--------+--------+
169 - | (x1-x0)*(y1-y0) |
170 +--------+--------+
171
172 COMMENT
173
174 FOR j FROM 2 * len DOWNTO 1
175 DO IF j > len
176 THEN z[j] := x0y0[j - len]
177 ELSE z[j] := x1y1[j]
178 FI
179 OD;
180
181 VOID (add (z, x1y1, len + half));
182 VOID (add (z, x0y0, len + half));
183
184 # Term (x1-x0)(y1-y0) was calculated as an absolute value,
185 and we add or subtract according to the sign. #
186
187 VOID ((p = q | sub | add) (z, xmym, len + half));
188
189 # Sanity check, should not fail #
190 IF [2 * len] INT w;
191 mul (w, x, y);
192 z /= w
193 THEN print (("Fail at length ", whole (len, 0), new line));
194 show ("x", x);
195 show ("y", y);
196 show ("w", w);
197 show ("z", z);
198 stop
199 FI
200 FI;
201
202
203 PROC karat test = (NUM x, y) VOID:
204 BEGIN [2 * UPB x] INT z;
205 printf ($67("-")l$);
206 show ("x", x); show ("y", y);
207 mul (z, x, y); show ("*", z);
208 karat (z, x, y); show ("^", z)
209 END;
210
211
212 INT radix = 1000; # Radix^2 must fit an INT. #
213
214 MODE NUM = [16] INT;
215
216 [] NUM tab = ((000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000),
217 (000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 000, 001),
218 (141, 421, 356, 237, 309, 504, 880, 168, 872, 420, 969, 807, 856, 967, 187, 538),
219 (161, 803, 398, 874, 989, 484, 820, 458, 683, 436, 563, 811, 772, 030, 917, 981),
220 (271, 828, 182, 845, 904, 523, 536, 028, 747, 135, 266, 249, 775, 724, 709, 370),
221 (314, 159, 265, 358, 979, 323, 846, 264, 338, 327, 950, 288, 419, 716, 939, 938),
222 (577, 215, 664, 901, 532, 860, 606, 512, 090, 082, 402, 431, 042, 159, 335, 940),
223 (999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999, 999));
224
225 FOR i TO UPB tab
226 DO FOR j TO UPB tab
227 DO karat test (tab[i], tab[j])
228 OD
229 OD
© J.M. van der Veer • jmvdveer@algol68genie.nl
