-
+ 8C306BD23D12315A8D62F760968B2CF32BDA6F49279B3F995DF7A100587C7B273E79C327C295521E42BEDA4D1E0A9F0E7F4C5293D57E92843B82C69A353F4687ffa/libffa/fz_barr.adb(0 . 0)(1 . 300)
216 ------------------------------------------------------------------------------
217 ------------------------------------------------------------------------------
218 -- This file is part of 'Finite Field Arithmetic', aka 'FFA'. --
219 -- --
220 -- (C) 2018 Stanislav Datskovskiy ( www.loper-os.org ) --
221 -- http://wot.deedbot.org/17215D118B7239507FAFED98B98228A001ABFFC7.html --
222 -- --
223 -- You do not have, nor can you ever acquire the right to use, copy or --
224 -- distribute this software ; Should you use this software for any purpose, --
225 -- or copy and distribute it to anyone or in any manner, you are breaking --
226 -- the laws of whatever soi-disant jurisdiction, and you promise to --
227 -- continue doing so for the indefinite future. In any case, please --
228 -- always : read and understand any software ; verify any PGP signatures --
229 -- that you use - for any purpose. --
230 -- --
231 -- See also http://trilema.com/2015/a-new-software-licensing-paradigm . --
232 ------------------------------------------------------------------------------
233
234 -----------------------------------------------------------------------------
235 -- BEFORE YOU EVEN *THINK* ABOUT MODIFYING THIS PROGRAM: --
236 -----------------------------------------------------------------------------
237 -- `dMMd` +NMMMMMMMNo --
238 -- .dM++Md..oNMMMMMMmo` --
239 -- /mM+ +MmmMMMMMMNo. --
240 -- /NM+ +MMMMMMNo` / --
241 -- /Nd- `sNMMMMMy.-oohNNm` --
242 -- `yNd- -yMMMMMMMNNNMMMMs. --
243 -- hMd::. -mMMMMMdNMMMMMMm+` --
244 -- :hNs.`.:yNyyyo--/sNMMMMy` --
245 -- -o.. `.. `.sNh` --
246 -- ..`RRR EEE AA DDD !+:` --
247 -- `: R R E A A D D ! .o- --
248 -- .s RRR EEE AAAA D D ! .:` --
249 -- .. `` R R E A A D D ys. --
250 -- -h /: R R EEE A A DDD !/ :mm- --
251 -- -mm `/- THE PROOFS!!! -y sMm- --
252 -- -mNy `++` YES THAT MEANS YOU! .s. .`-Nm- --
253 -- `oNN-`:/y``:////:` `-////``-o.+. -NNo` --
254 -- `oNN: `+:::hNMMMMNyo. `smNMMMMmy`++ :NNo` --
255 -- `oNy- .: sMMMMMMMMM: -MMMMMMMMMs/s. -yNh- --
256 -- -dNy. `s. sMMMMMMMmo.----mMMMMMMMNo `.` .yMd- --
257 -- .dmo. `o` /mNNNmyo.`sNMMy.+ymNNNNh `-` .omd. --
258 -- .mN/ -o` .---. `oNMNMNs .----. -/. /Nm. --
259 -- +mN/ .hhs:.. ` .hMN-MMy ` `.-+-` /Nm+ --
260 -- +NN: :hMMMs/m`d -y- dy -`:/y :NN+ --
261 -- +Nd: /: `hMMMmm/ y:Ns::.`````.:oh-- :dNs. --
262 -- .sNh. .h+:hMMMy./- -yMMMyyod+dssMM:. `: .hMh. --
263 -- .hMy. +MNMMMMh` ` `yNMhmsNsmhNh: /` +Mh. --
264 -- -hN+ -dMMMMMNso+- :s .ymmNMmyh+- + +Nh- --
265 -- `MN+ /MMMMMMh:- :- :: : .+ +NM` --
266 -- `Md///////+mMMMMys////////sh/- -yy/////////////////////////dM` --
267 -- -ssssssssymssssssssssssssssso- .+ossssssssssssssssssssssssssss- --
268 -- --
269 --Ch. 14A: Barrett’s Modular Reduction: http://www.loper-os.org/?p=2842 --
270 --Ch. 14A-Bis: Barrett’s Physical Bounds: http://www.loper-os.org/?p=2875 --
271 -- --
272 -----------------------------------------------------------------------------
273
274 with W_Pred; use W_Pred;
275 with W_Shifts; use W_Shifts;
276 with FZ_Basic; use FZ_Basic;
277 with FZ_Shift; use FZ_Shift;
278 with FZ_Arith; use FZ_Arith;
279 with FZ_Mul; use FZ_Mul;
280 with FZ_LoMul; use FZ_LoMul;
281 with FZ_Measr; use FZ_Measr;
282 with FZ_QShft; use FZ_QShft;
283
284
285 package body FZ_Barr is
286
287 -- Prepare the precomputed Barrettoid corresponding to a given Modulus
288 procedure FZ_Make_Barrettoid(Modulus : in FZ;
289 Result : out Barretoid) is
290
291 -- Length of Modulus and Remainder
292 Lm : constant Indices := Modulus'Length;
293
294 -- Remainder register, starts as zero
295 Remainder : FZ(1 .. Lm) := (others => 0);
296
297 -- Length of Quotient, with an extra Word for top bit (if Degenerate)
298 Lq : constant Indices := (2 * Lm) + 1;
299
300 -- Valid indices into Quotient, using the above
301 subtype Quotient_Index is Word_Index range 1 .. Lq;
302
303 -- The Quotient we need, i.e. 2^(2 * ModulusBitness) / Modulus
304 Quotient : FZ(Quotient_Index);
305
306 -- Permissible 'cuts' for the Slice operation
307 subtype Divisor_Cuts is Word_Index range 2 .. Lm;
308
309 -- Current bit of Pseudo-Dividend; high bit is 1, all others 0
310 Pb : WBool := 1;
311
312 -- Performs Restoring Division on a given segment
313 procedure Slice(Index : Quotient_Index;
314 Cut : Divisor_Cuts;
315 Bits : Positive) is
316 begin
317
318 declare
319
320 -- Borrow, from comparator
321 C : WBool;
322
323 -- Left-Shift Overflow
324 LsO : WBool;
325
326 -- Current cut of Remainder register
327 Rs : FZ renames Remainder(1 .. Cut);
328
329 -- Current cut of Divisor
330 Ds : FZ renames Modulus(1 .. Cut);
331
332 -- Current Word of Quotient being made, starting from the highest
333 W : Word := 0;
334
335 -- Current bit of Quotient (inverted)
336 nQb : WBool;
337
338 begin
339
340 -- For each bit in the current Pseudo-Dividend Word:
341 for b in 1 .. Bits loop
342
343 -- Advance Rs, shifting in the current Pseudo-Dividend bit:
344 FZ_ShiftLeft_O_I(N => Rs,
345 ShiftedN => Rs,
346 Count => 1,
347 OF_In => Pb, -- Current Pseudo-Dividend bit
348 Overflow => LsO);
349
350 -- Subtract Divisor-Cut from R-Cut; Underflow goes into C:
351 FZ_Sub(X => Rs, Y => Ds, Difference => Rs, Underflow => C);
352
353 -- Negation of current Quotient bit
354 nQb := C and W_Not(LsO);
355
356 -- If C=1, the subtraction underflowed, and we must undo it:
357 FZ_Add_Gated(X => Rs, Y => Ds, Sum => Rs,
358 Gate => nQb);
359
360 -- Save the bit of Quotient that we have found:
361 W := Shift_Left(W, 1) or (1 - nQb); -- i.e. inverse of nQb
362
363 end loop;
364
365 -- We made a complete Word of the Quotient; save it:
366 Quotient(Quotient'Last + 1 - Index) := W; -- Indexed from end
367
368 end;
369
370 end Slice;
371 pragma Inline_Always(Slice);
372
373 -- Measure of the Modulus
374 ModulusMeasure : constant FZBit_Index := FZ_Measure(Modulus);
375
376 begin
377
378 -- First, process the high Word of the Pseudo-Dividend:
379 Slice(1, 2, 1); -- ... it has just one bit: a 1 at the bottom position
380
381 -- Once we ate the top 1 bit of Pseudo-Dividend, rest of its bits are 0
382 Pb := 0;
383
384 -- Process the Modulus-sized segment below the top Word:
385 for i in 2 .. Lm - 1 loop
386
387 Slice(i, i + 1, Bitness); -- stay ahead by a word to handle carry
388
389 end loop;
390
391 -- Process remaining Words:
392 for i in Lm .. Lq loop
393
394 Slice(i, Lm, Bitness);
395
396 end loop;
397
398 -- Record the Quotient (i.e. the Barrettoid proper, Bm)
399 Result.B := Quotient(Result.B'Range);
400
401 -- Record whether we have the Degenerate Case (1 iff Modulus = 1)
402 Result.Degenerate := W_NZeroP(Quotient(Lq));
403
404 -- Record a copy of the Modulus, extended with zero Word:
405 Result.ZXM(Modulus'Range) := Modulus;
406 Result.ZXM(Result.ZXM'Last) := 0;
407
408 -- Record the parameter Jm:
409 Result.J := ModulusMeasure - 1;
410
411 -- Wm - Jm
412 Result.ZSlide :=
413 FZBit_Index(Bitness * Modulus'Length) - ModulusMeasure + 1;
414
415 end FZ_Make_Barrettoid;
416
417
418 -- Reduce X using the given precomputed Barrettoid.
419 procedure FZ_Barrett_Reduce(X : in FZ;
420 Bar : in Barretoid;
421 XReduced : in out FZ) is
422
423 -- Wordness of X, the quantity being reduced
424 Xl : constant Indices := X'Length;
425
426 -- Wordness of XReduced (result), one-half of Xl, and same as of Modulus
427 Ml : constant Indices := XReduced'Length; -- i.e. # of Words in Wm.
428
429 -- The Modulus we will reduce X by
430 Modulus : FZ renames Bar.ZXM(1 .. Ml);
431
432 -- Shifted X
433 Xs : FZ(X'Range);
434
435 -- Z := Xs * Bm (has twice the length of X)
436 Z : FZ(1 .. 2 * Xl);
437
438 -- Upper 3Wm-bit segment of Z that gets shifted to form Zs
439 ZHi : FZ renames Z(Ml + 1 .. Z'Last );
440
441 -- Middle 2Wm-bit segment of Z, that gets multiplied by M to form Q
442 Zs : FZ renames Z(Z'First + Ml .. Z'Last - Ml );
443
444 -- Sub-terms of the Zs * M multiplication:
445 ZsLo : FZ renames Zs(Zs'First .. Zs'Last - Ml );
446 ZsHi : FZ renames Zs(Zs'First + Ml .. Zs'Last );
447 ZsHiM : FZ(1 .. Ml);
448
449 -- Q := Modulus * Zs, i.e. floor(X / M)*M + E*M
450 Q : FZ(1 .. Xl);
451 QHi : FZ renames Q(Q'First + Ml .. Q'Last );
452
453 -- R is made one Word longer than Modulus (see proof re: why)
454 Rl : constant Indices := Ml + 1;
455
456 -- Reduction estimate, overshot by 0, 1, or 2 multiple of Modulus
457 R : FZ(1 .. Rl);
458
459 -- Scratch for the outputs of the gated subtractions
460 S : FZ(1 .. Rl);
461
462 -- Borrow from the gated subtractions
463 C : WBool;
464
465 -- Barring cosmic ray, no underflow can take place in (4) and (5)
466 NoCarry : WZeroOrDie := 0;
467
468 begin
469
470 -- Result is initially zero (and will stay zero if Modulus = 1)
471 FZ_Clear(XReduced);
472
473 -- (1) Ns := X >> Jm
474 FZ_Quiet_ShiftRight(N => X, ShiftedN => Xs, Count => Bar.J);
475
476 -- (2) Z := X * Bm
477 FZ_Multiply_Unbuffered(X => Bar.B, Y => Xs, XY => Z);
478
479 -- (3) Zs := Z >> 2Wm - Jm (already thrown lower Wm, so only Wm - Jm now)
480 FZ_Quiet_ShiftRight(N => ZHi, ShiftedN => ZHi, Count => Bar.ZSlide);
481
482 -- (4) Q := Zs * M (this is split into three operations, see below)
483
484 -- ... first, Q := ZsLo * M,
485 FZ_Multiply_Unbuffered(ZsLo, Modulus, Q);
486
487 -- ... then, compute ZsHiM := ZsHi * M,
488 FZ_Low_Multiply_Unbuffered(ZsHi, Modulus, ZsHiM);
489
490 -- ... finally, add ZsHiM to upper half of Q.
491 FZ_Add_D(X => QHi, Y => ZsHiM, Overflow => NoCarry);
492
493 -- (5) R := X - Q (we only need Rl-sized segments of X and Q here)
494 FZ_Sub(X => X(1 .. Rl), Y => Q(1 .. Rl),
495 Difference => R, Underflow => NoCarry);
496
497 -- (6) S1 := R - M, C1 := Borrow (1st gated subtraction of Modulus)
498 FZ_Sub(X => R, Y => Bar.ZXM, Difference => S, Underflow => C);
499
500 -- (7) R := {C1=0 -> S1, C1=1 -> R}
501 FZ_Mux(X => S, Y => R, Result => R, Sel => C);
502
503 -- (8) S2 := R - M, C2 := Borrow (2nd gated subtraction of Modulus)
504 FZ_Sub(X => R, Y => Bar.ZXM, Difference => S, Underflow => C);
505
506 -- (9) R := {C2=0 -> S2, C2=1 -> R}
507 FZ_Mux(X => S, Y => R, Result => R, Sel => C);
508
509 -- (10) RFinal := {DM=0 -> R, DM=1 -> 0} (handle the degenerate case)
510 FZ_Mux(X => R(1 .. Ml), Y => XReduced, Result => XReduced,
511 Sel => Bar.Degenerate); -- If Modulus = 1, then XReduced is 0.
512
513 end FZ_Barrett_Reduce;
514
515 end FZ_Barr;