File : fz_barr.adb
1 ------------------------------------------------------------------------------
2 ------------------------------------------------------------------------------
3 -- This file is part of 'Finite Field Arithmetic', aka 'FFA'. --
4 -- --
5 -- (C) 2019 Stanislav Datskovskiy ( www.loper-os.org ) --
6 -- http://wot.deedbot.org/17215D118B7239507FAFED98B98228A001ABFFC7.html --
7 -- --
8 -- You do not have, nor can you ever acquire the right to use, copy or --
9 -- distribute this software ; Should you use this software for any purpose, --
10 -- or copy and distribute it to anyone or in any manner, you are breaking --
11 -- the laws of whatever soi-disant jurisdiction, and you promise to --
12 -- continue doing so for the indefinite future. In any case, please --
13 -- always : read and understand any software ; verify any PGP signatures --
14 -- that you use - for any purpose. --
15 -- --
16 -- See also http://trilema.com/2015/a-new-software-licensing-paradigm . --
17 ------------------------------------------------------------------------------
18
19 -----------------------------------------------------------------------------
20 -- BEFORE YOU EVEN *THINK* ABOUT MODIFYING THIS PROGRAM: --
21 -----------------------------------------------------------------------------
22 -- `dMMd` +NMMMMMMMNo --
23 -- .dM++Md..oNMMMMMMmo` --
24 -- /mM+ +MmmMMMMMMNo. --
25 -- /NM+ +MMMMMMNo` / --
26 -- /Nd- `sNMMMMMy.-oohNNm` --
27 -- `yNd- -yMMMMMMMNNNMMMMs. --
28 -- hMd::. -mMMMMMdNMMMMMMm+` --
29 -- :hNs.`.:yNyyyo--/sNMMMMy` --
30 -- -o.. `.. `.sNh` --
31 -- ..`RRR EEE AA DDD !+:` --
32 -- `: R R E A A D D ! .o- --
33 -- .s RRR EEE AAAA D D ! .:` --
34 -- .. `` R R E A A D D ys. --
35 -- -h /: R R EEE A A DDD !/ :mm- --
36 -- -mm `/- THE PROOFS!!! -y sMm- --
37 -- -mNy `++` YES THAT MEANS YOU! .s. .`-Nm- --
38 -- `oNN-`:/y``:////:` `-////``-o.+. -NNo` --
39 -- `oNN: `+:::hNMMMMNyo. `smNMMMMmy`++ :NNo` --
40 -- `oNy- .: sMMMMMMMMM: -MMMMMMMMMs/s. -yNh- --
41 -- -dNy. `s. sMMMMMMMmo.----mMMMMMMMNo `.` .yMd- --
42 -- .dmo. `o` /mNNNmyo.`sNMMy.+ymNNNNh `-` .omd. --
43 -- .mN/ -o` .---. `oNMNMNs .----. -/. /Nm. --
44 -- +mN/ .hhs:.. ` .hMN-MMy ` `.-+-` /Nm+ --
45 -- +NN: :hMMMs/m`d -y- dy -`:/y :NN+ --
46 -- +Nd: /: `hMMMmm/ y:Ns::.`````.:oh-- :dNs. --
47 -- .sNh. .h+:hMMMy./- -yMMMyyod+dssMM:. `: .hMh. --
48 -- .hMy. +MNMMMMh` ` `yNMhmsNsmhNh: /` +Mh. --
49 -- -hN+ -dMMMMMNso+- :s .ymmNMmyh+- + +Nh- --
50 -- `MN+ /MMMMMMh:- :- :: : .+ +NM` --
51 -- `Md///////+mMMMMys////////sh/- -yy/////////////////////////dM` --
52 -- -ssssssssymssssssssssssssssso- .+ossssssssssssssssssssssssssss- --
53 -- --
54 --Ch. 14A: Barrett's Modular Reduction: http://www.loper-os.org/?p=2842 --
55 --Ch. 14A-Bis: Barrett's Physical Bounds: http://www.loper-os.org/?p=2875 --
56 -- --
57 -----------------------------------------------------------------------------
58
59 with W_Pred; use W_Pred;
60 with W_Shifts; use W_Shifts;
61 with FZ_Basic; use FZ_Basic;
62 with FZ_Shift; use FZ_Shift;
63 with FZ_Arith; use FZ_Arith;
64 with FZ_Mul; use FZ_Mul;
65 with FZ_LoMul; use FZ_LoMul;
66 with FZ_Measr; use FZ_Measr;
67 with FZ_QShft; use FZ_QShft;
68
69
70 package body FZ_Barr is
71
72 -- Prepare the precomputed Barrettoid corresponding to a given Modulus
73 procedure FZ_Make_Barrettoid(Modulus : in FZ;
74 Result : out Barretoid) is
75
76 -- Length of Modulus and Remainder
77 Lm : constant Indices := Modulus'Length;
78
79 -- Remainder register, starts as zero
80 Remainder : FZ(1 .. Lm) := (others => 0);
81
82 -- Length of Quotient, with an extra Word for top bit (if Degenerate)
83 Lq : constant Indices := (2 * Lm) + 1;
84
85 -- Valid indices into Quotient, using the above
86 subtype Quotient_Index is Word_Index range 1 .. Lq;
87
88 -- The Quotient we need, i.e. 2^(2 * ModulusBitness) / Modulus
89 Quotient : FZ(Quotient_Index);
90
91 -- Permissible 'cuts' for the Slice operation
92 subtype Divisor_Cuts is Word_Index range 2 .. Lm;
93
94 -- Current bit of Pseudo-Dividend; high bit is 1, all others 0
95 Pb : WBool := 1;
96
97 -- Performs Restoring Division on a given segment
98 procedure Slice(Index : Quotient_Index;
99 Cut : Divisor_Cuts;
100 Bits : Positive) is
101 begin
102
103 declare
104
105 -- Borrow, from comparator
106 C : WBool;
107
108 -- Left-Shift Overflow
109 LsO : WBool;
110
111 -- Current cut of Remainder register
112 Rs : FZ renames Remainder(1 .. Cut);
113
114 -- Current cut of Divisor
115 Ds : FZ renames Modulus(1 .. Cut);
116
117 -- Current Word of Quotient being made, starting from the highest
118 W : Word := 0;
119
120 -- Current bit of Quotient (inverted)
121 nQb : WBool;
122
123 begin
124
125 -- For each bit in the current Pseudo-Dividend Word:
126 for b in 1 .. Bits loop
127
128 -- Advance Rs, shifting in the current Pseudo-Dividend bit:
129 FZ_ShiftLeft_O_I(N => Rs,
130 ShiftedN => Rs,
131 Count => 1,
132 OF_In => Pb, -- Current Pseudo-Dividend bit
133 Overflow => LsO);
134
135 -- Subtract Divisor-Cut from R-Cut; Underflow goes into C:
136 FZ_Sub(X => Rs, Y => Ds, Difference => Rs, Underflow => C);
137
138 -- Negation of current Quotient bit
139 nQb := C and W_Not(LsO);
140
141 -- If C=1, the subtraction underflowed, and we must undo it:
142 FZ_Add_Gated(X => Rs, Y => Ds, Sum => Rs,
143 Gate => nQb);
144
145 -- Save the bit of Quotient that we have found:
146 W := Shift_Left(W, 1) or (1 - nQb); -- i.e. inverse of nQb
147
148 end loop;
149
150 -- We made a complete Word of the Quotient; save it:
151 Quotient(Quotient'Last + 1 - Index) := W; -- Indexed from end
152
153 end;
154
155 end Slice;
156 pragma Inline_Always(Slice);
157
158 -- Measure of the Modulus
159 ModulusMeasure : constant FZBit_Index := FZ_Measure(Modulus);
160
161 begin
162
163 -- First, process the high Word of the Pseudo-Dividend:
164 Slice(1, 2, 1); -- ... it has just one bit: a 1 at the bottom position
165
166 -- Once we ate the top 1 bit of Pseudo-Dividend, rest of its bits are 0
167 Pb := 0;
168
169 -- Process the Modulus-sized segment below the top Word:
170 for i in 2 .. Lm - 1 loop
171
172 Slice(i, i + 1, Bitness); -- stay ahead by a word to handle carry
173
174 end loop;
175
176 -- Process remaining Words:
177 for i in Lm .. Lq loop
178
179 Slice(i, Lm, Bitness);
180
181 end loop;
182
183 -- Record the Quotient (i.e. the Barrettoid proper, Bm)
184 Result.B := Quotient(Result.B'Range);
185
186 -- Record whether we have the Degenerate Case (1 iff Modulus = 1)
187 Result.Degenerate := W_NZeroP(Quotient(Lq));
188
189 -- Record a copy of the Modulus, extended with zero Word:
190 Result.ZXM(Modulus'Range) := Modulus;
191 Result.ZXM(Result.ZXM'Last) := 0;
192
193 -- Record the parameter Jm:
194 Result.J := ModulusMeasure - 1;
195
196 -- Wm - Jm
197 Result.ZSlide :=
198 FZBit_Index(Bitness * Modulus'Length) - ModulusMeasure + 1;
199
200 end FZ_Make_Barrettoid;
201
202
203 -- Reduce X using the given precomputed Barrettoid.
204 procedure FZ_Barrett_Reduce(X : in FZ;
205 Bar : in Barretoid;
206 XReduced : in out FZ) is
207
208 -- Wordness of X, the quantity being reduced
209 Xl : constant Indices := X'Length;
210
211 -- Wordness of XReduced (result), one-half of Xl, and same as of Modulus
212 Ml : constant Indices := XReduced'Length; -- i.e. # of Words in Wm.
213
214 -- The Modulus we will reduce X by
215 Modulus : FZ renames Bar.ZXM(1 .. Ml);
216
217 -- Shifted X
218 Xs : FZ(X'Range);
219
220 -- Z := Xs * Bm (has twice the length of X)
221 Z : FZ(1 .. 2 * Xl);
222
223 -- Upper 3Wm-bit segment of Z that gets shifted to form Zs
224 ZHi : FZ renames Z(Ml + 1 .. Z'Last );
225
226 -- Middle 2Wm-bit segment of Z, that gets multiplied by M to form Q
227 Zs : FZ renames Z(Z'First + Ml .. Z'Last - Ml );
228
229 -- Sub-terms of the Zs * M multiplication:
230 ZsLo : FZ renames Zs(Zs'First .. Zs'Last - Ml );
231 ZsHi : FZ renames Zs(Zs'First + Ml .. Zs'Last );
232 ZsHiM : FZ(1 .. Ml);
233
234 -- Q := Modulus * Zs, i.e. floor(X / M)*M + E*M
235 Q : FZ(1 .. Xl);
236 QHi : FZ renames Q(Q'First + Ml .. Q'Last );
237
238 -- R is made one Word longer than Modulus (see proof re: why)
239 Rl : constant Indices := Ml + 1;
240
241 -- Reduction estimate, overshot by 0, 1, or 2 multiple of Modulus
242 R : FZ(1 .. Rl);
243
244 -- Scratch for the outputs of the gated subtractions
245 S : FZ(1 .. Rl);
246
247 -- Borrow from the gated subtractions
248 C : WBool;
249
250 -- Barring cosmic ray, no underflow can take place in (4) and (5)
251 NoCarry : WZeroOrDie := 0;
252
253 begin
254
255 -- Result is initially zero (and will stay zero if Modulus = 1)
256 FZ_Clear(XReduced);
257
258 -- (1) Xs := X >> Jm
259 FZ_Quiet_ShiftRight(N => X, ShiftedN => Xs, Count => Bar.J);
260
261 -- (2) Z := Xs * Bm
262 FZ_Multiply_Unbuffered(X => Bar.B, Y => Xs, XY => Z);
263
264 -- (3) Zs := Z >> 2Wm - Jm (already thrown lower Wm, so only Wm - Jm now)
265 FZ_Quiet_ShiftRight(N => ZHi, ShiftedN => ZHi, Count => Bar.ZSlide);
266
267 -- (4) Q := Zs * M (this is split into three operations, see below)
268
269 -- ... first, Q := ZsLo * M,
270 FZ_Multiply_Unbuffered(ZsLo, Modulus, Q);
271
272 -- ... then, compute ZsHiM := ZsHi * M,
273 FZ_Low_Multiply_Unbuffered(ZsHi, Modulus, ZsHiM);
274
275 -- ... finally, add ZsHiM to upper half of Q.
276 FZ_Add_D(X => QHi, Y => ZsHiM, Overflow => NoCarry);
277
278 -- (5) R := X - Q (we only need Rl-sized segments of X and Q here)
279 FZ_Sub(X => X(1 .. Rl), Y => Q(1 .. Rl),
280 Difference => R, Underflow => NoCarry);
281
282 -- (6) S1 := R - M, C1 := Borrow (1st gated subtraction of Modulus)
283 FZ_Sub(X => R, Y => Bar.ZXM, Difference => S, Underflow => C);
284
285 -- (7) R := {C1=0 -> S1, C1=1 -> R}
286 FZ_Mux(X => S, Y => R, Result => R, Sel => C);
287
288 -- (8) S2 := R - M, C2 := Borrow (2nd gated subtraction of Modulus)
289 FZ_Sub(X => R, Y => Bar.ZXM, Difference => S, Underflow => C);
290
291 -- (9) R := {C2=0 -> S2, C2=1 -> R}
292 FZ_Mux(X => S, Y => R, Result => R, Sel => C);
293
294 -- (10) RFinal := {DM=0 -> R, DM=1 -> 0} (handle the degenerate case)
295 FZ_Mux(X => R(1 .. Ml), Y => XReduced, Result => XReduced,
296 Sel => Bar.Degenerate); -- If Modulus = 1, then XReduced is 0.
297
298 end FZ_Barrett_Reduce;
299
300 end FZ_Barr;