- 5CB9ECB938F842B7C34CA7EDF99FB92502986D7F660B223659C9B19C6008A24C3BE6DE8135DB37E29E7FF278A451B68C61F6B57D8E404372DCD3FD6D4320EA0A
+ 98A9315BADAC1530D96E214873BE3EFCED853C4A1D922CA1D707D93E62993ADDFFE906A64525204C9FB83B38C9AC04A170711D8D10D33D3EE2699953DDED18B9
ffa/libffa/fz_mul.adb
(20 . 15)(20 . 15)
204 with Words; use Words;
205 with Word_Ops; use Word_Ops;
206 with W_Mul; use W_Mul;
207 with FZ_Arith; use FZ_Arith;
208
209
210 package body FZ_Mul is
211
212 -- Comba's multiplier.
213
214 -- Comba's multiplier. (CAUTION: UNBUFFERED)
215 procedure FZ_Mul_Comba(X : in FZ;
216 Y : in FZ;
217 XY_Lo : out FZ;
218 XY_Hi : out FZ) is
219 XY : out FZ) is
220
221 -- Words in each multiplicand
222 L : constant Word_Index := X'Length;
(39 . 11)(39 . 8)
224 -- 3-word Accumulator
225 A2, A1, A0 : Word := 0;
226
227 -- Register holding Product; indexed from zero
228 XY : FZ(0 .. LP - 1);
229
230 -- Type for referring to a column of XY
231 subtype ColXY is Word_Index range XY'Range;
232 subtype ColXY is Word_Index range 0 .. LP - 1;
233
234 -- Compute the Nth (indexed from zero) column of the Product
235 procedure Col(N : in ColXY; U : in ColXY; V : in ColXY) is
(86 . 7)(83 . 7)
237 end loop;
238
239 -- We now have the Nth (indexed from zero) word of XY
240 XY(N) := A0;
241 XY(XY'First + N) := A0;
242
243 -- Right-Shift the Accumulator by one Word width:
244 A0 := A1;
(115 . 11)(112 . 163)
246 -- The very last Word of the Product:
247 XY(XY'Last) := A0;
248
249 -- Output the Product's lower and upper FZs:
250 XY_Lo := XY(0 .. L - 1);
251 XY_Hi := XY(L .. XY'Last);
252
253 end FZ_Mul_Comba;
254 pragma Inline_Always(FZ_Mul_Comba);
255
256
257 -- Karatsuba's Multiplier. (CAUTION: UNBUFFERED)
258 procedure Mul_Karatsuba(X : in FZ;
259 Y : in FZ;
260 XY : out FZ) is
261
262 -- L is the wordness of a multiplicand. Guaranteed to be a power of two.
263 L : constant Word_Count := X'Length;
264
265 -- An 'LSeg' is the same length as either multiplicand.
266 subtype LSeg is FZ(1 .. L);
267
268 -- K is HALF of the length of a multiplicand.
269 K : constant Word_Index := L / 2;
270
271 -- A 'KSeg' is the same length as HALF of a multiplicand.
272 subtype KSeg is FZ(1 .. K);
273
274 -- The three L-sized variables of the product equation, i.e.:
275 -- XY = LL + 2^(K*Bitness)(LL + HH + (-1^DD_Sub)*DD) + 2^(2*K*Bitness)HH
276 LL, DD, HH : LSeg;
277
278 -- K-sized terms of Dx * Dy = DD
279 Dx, Dy : KSeg; -- Dx = abs(XLo - XHi) , Dy = abs(YLo - YHi)
280
281 -- Subtraction borrows, signs of (XL - XH) and (YL - YH),
282 Cx, Cy : WBool; -- so that we can calculate (-1^DD_Sub)
283
284 -- Bottom and Top K-sized halves of the multiplicand X.
285 XLo : KSeg renames X( X'First .. X'Last - K );
286 XHi : KSeg renames X( X'First + K .. X'Last );
287
288 -- Bottom and Top K-sized halves of the multiplicand Y.
289 YLo : KSeg renames Y( Y'First .. Y'Last - K );
290 YHi : KSeg renames Y( Y'First + K .. Y'Last );
291
292 -- L-sized middle segment of the product XY (+/- K from the midpoint).
293 XYMid : LSeg renames XY( XY'First + K .. XY'Last - K );
294
295 -- Bottom and Top L-sized halves of the product XY.
296 XYLo : LSeg renames XY( XY'First .. XY'Last - L );
297 XYHi : LSeg renames XY( XY'First + L .. XY'Last );
298
299 -- Topmost K-sized quarter segment of the product XY, or 'tail'
300 XYHiHi : KSeg renames XYHi( XYHi'First + K .. XYHi'Last );
301
302 -- Whether the DD term is being subtracted.
303 DD_Sub : WBool;
304
305 -- Carry from individual term additions.
306 C : WBool;
307
308 -- Tail-Carry accumulator, for the final ripple
309 TC : Word;
310
311 begin
312
313 -- Recurse: LL := XL * YL
314 FZ_Multiply(XLo, YLo, LL);
315
316 -- Recurse: HH := XH * YH
317 FZ_Multiply(XHi, YHi, HH);
318
319 -- Dx := |XL - XH| , Cx := Borrow (i.e. 1 iff XL < XH)
320 FZ_Sub_Abs(X => XLo, Y => XHi, Difference => Dx, Underflow => Cx);
321
322 -- Dy := |YL - YH| , Cy := Borrow (i.e. 1 iff YL < YH)
323 FZ_Sub_Abs(X => YLo, Y => YHi, Difference => Dy, Underflow => Cy);
324
325 -- Recurse: DD := Dx * Dy
326 FZ_Multiply(Dx, Dy, DD);
327
328 -- Whether (XL - XH)(YL - YH) is positive, and so DD must be subtracted:
329 DD_Sub := 1 - (Cx xor Cy);
330
331 -- XY := LL + 2^(2 * K * Bitness) * HH
332 XYLo := LL;
333 XYHi := HH;
334
335 -- XY += 2^(K * Bitness) * HH, but carry goes in Tail Carry accum.
336 FZ_Add_D(X => XYMid, Y => HH, Overflow => TC);
337
338 -- XY += 2^(K * Bitness) * LL, ...
339 FZ_Add_D(X => XYMid, Y => LL, Overflow => C);
340
341 -- ... but the carry goes into the Tail Carry accumulator.
342 TC := TC + C;
343
344 -- XY += 2^(K * Bitness) * (-1^DD_Sub) * DD
345 FZ_Not_Cond_D(N => DD, Cond => DD_Sub); -- invert DD if 2s-complementing
346 FZ_Add_D(OF_In => DD_Sub, -- ... and then must increment
347 X => XYMid,
348 Y => DD,
349 Overflow => C); -- carry will go in Tail Carry accumulator
350
351 -- Compute the final Tail Carry for the ripple
352 TC := TC + C - DD_Sub;
353
354 -- Barring a cosmic ray, 0 <= TC <= 2 .
355 pragma Assert(TC <= 2);
356
357 -- Ripple the Tail Carry into the tail.
358 FZ_Add_D_W(X => XYHiHi, W => TC, Overflow => C);
359
360 -- Barring a cosmic ray, the tail ripple will NOT overflow.
361 pragma Assert(C = 0);
362
363 end Mul_Karatsuba;
364 -- CAUTION: Inlining prohibited for Mul_Karatsuba !
365
366
367 -- Multiplier. (CAUTION: UNBUFFERED)
368 procedure FZ_Multiply(X : in FZ;
369 Y : in FZ;
370 XY : out FZ) is
371
372 -- The length of either multiplicand
373 L : constant Word_Count := X'Length;
374
375 begin
376
377 if L <= Karatsuba_Thresh then
378
379 -- Base case:
380 FZ_Mul_Comba(X, Y, XY);
381
382 else
383
384 -- Recursive case:
385 Mul_Karatsuba(X, Y, XY);
386
387 end if;
388
389 end FZ_Multiply;
390 pragma Inline_Always(FZ_Multiply);
391
392
393 -- Multiplier. Preserves the inputs.
394 procedure FZ_Mult(X : in FZ;
395 Y : in FZ;
396 XY_Lo : out FZ;
397 XY_Hi : out FZ) is
398
399 -- Product buffer.
400 P : FZ(1 .. 2 * X'Length);
401
402 begin
403
404 FZ_Multiply(X, Y, P);
405
406 XY_Lo := P(P'First .. P'First + X'Length - 1);
407 XY_Hi := P(P'First + X'Length .. P'Last);
408
409 end FZ_Mult;
410 pragma Inline_Always(FZ_Mult);
411
412 end FZ_Mul;