raw
ffa_ch12_karatsub...    1 ------------------------------------------------------------------------------
ffa_ch12_karatsub... 2 ------------------------------------------------------------------------------
ffa_ch12_karatsub... 3 -- This file is part of 'Finite Field Arithmetic', aka 'FFA'. --
ffa_ch12_karatsub... 4 -- --
ffa_ch12_karatsub... 5 -- (C) 2018 Stanislav Datskovskiy ( www.loper-os.org ) --
ffa_ch12_karatsub... 6 -- http://wot.deedbot.org/17215D118B7239507FAFED98B98228A001ABFFC7.html --
ffa_ch12_karatsub... 7 -- --
ffa_ch12_karatsub... 8 -- You do not have, nor can you ever acquire the right to use, copy or --
ffa_ch12_karatsub... 9 -- distribute this software ; Should you use this software for any purpose, --
ffa_ch12_karatsub... 10 -- or copy and distribute it to anyone or in any manner, you are breaking --
ffa_ch12_karatsub... 11 -- the laws of whatever soi-disant jurisdiction, and you promise to --
ffa_ch12_karatsub... 12 -- continue doing so for the indefinite future. In any case, please --
ffa_ch12_karatsub... 13 -- always : read and understand any software ; verify any PGP signatures --
ffa_ch12_karatsub... 14 -- that you use - for any purpose. --
ffa_ch12_karatsub... 15 -- --
ffa_ch12_karatsub... 16 -- See also http://trilema.com/2015/a-new-software-licensing-paradigm . --
ffa_ch12_karatsub... 17 ------------------------------------------------------------------------------
ffa_ch12_karatsub... 18 ------------------------------------------------------------------------------
ffa_ch12_karatsub... 19
ffa_ch12_karatsub... 20 with Words; use Words;
ffa_ch12_karatsub... 21 with Word_Ops; use Word_Ops;
ffa_ch12_karatsub... 22 with W_Mul; use W_Mul;
ffa_ch12_karatsub... 23 with FZ_Arith; use FZ_Arith;
ffa_ch12_karatsub... 24
ffa_ch12_karatsub... 25
ffa_ch12_karatsub... 26 package body FZ_Sqr is
ffa_ch12_karatsub... 27
ffa_ch12_karatsub... 28 -- Square case of Comba's multiplier. (CAUTION: UNBUFFERED)
ffa_ch12_karatsub... 29 procedure FZ_Sqr_Comba(X : in FZ;
ffa_ch12_karatsub... 30 XX : out FZ) is
ffa_ch12_karatsub... 31
ffa_ch12_karatsub... 32 -- Words in each multiplicand
ffa_ch12_karatsub... 33 L : constant Word_Index := X'Length;
ffa_ch12_karatsub... 34
ffa_ch12_karatsub... 35 -- Length of Product, i.e. double the length of X
ffa_ch12_karatsub... 36 LP : constant Word_Index := 2 * L;
ffa_ch12_karatsub... 37
ffa_ch12_karatsub... 38 -- 3-word Accumulator
ffa_ch12_karatsub... 39 A2, A1, A0 : Word := 0;
ffa_ch12_karatsub... 40
ffa_ch12_karatsub... 41 Lo, Hi : Word; -- Output of WxW multiply/square
ffa_ch12_karatsub... 42
ffa_ch12_karatsub... 43 -- Type for referring to a column of XX
ffa_ch12_karatsub... 44 subtype ColXX is Word_Index range 0 .. LP - 1;
ffa_ch12_karatsub... 45
ffa_ch12_karatsub... 46 procedure Accum is
ffa_ch12_karatsub... 47 C : WBool; -- Carry for the Accumulator addition
ffa_ch12_karatsub... 48 Sum : Word; -- Sum for Accumulator addition
ffa_ch12_karatsub... 49 begin
ffa_ch12_karatsub... 50 -- Now add add double-Word Hi:Lo to accumulator A2:A1:A0:
ffa_ch12_karatsub... 51 -- A0 += Lo; C := Carry
ffa_ch12_karatsub... 52 Sum := A0 + Lo;
ffa_ch12_karatsub... 53 C := W_Carry(A0, Lo, Sum);
ffa_ch12_karatsub... 54 A0 := Sum;
ffa_ch12_karatsub... 55 -- A1 += Hi + C; C := Carry
ffa_ch12_karatsub... 56 Sum := A1 + Hi + C;
ffa_ch12_karatsub... 57 C := W_Carry(A1, Hi, Sum);
ffa_ch12_karatsub... 58 A1 := Sum;
ffa_ch12_karatsub... 59 -- A2 += A2 + C
ffa_ch12_karatsub... 60 A2 := A2 + C;
ffa_ch12_karatsub... 61 end Accum;
ffa_ch12_karatsub... 62 pragma Inline_Always(Accum);
ffa_ch12_karatsub... 63
ffa_ch12_karatsub... 64 procedure SymmDigits(N : in ColXX; From : in ColXX; To : in ColXX) is
ffa_ch12_karatsub... 65 begin
ffa_ch12_karatsub... 66 for j in From .. To loop
ffa_ch12_karatsub... 67 -- Hi:Lo := j-th * (N - j)-th Word, and then,
ffa_ch12_karatsub... 68 Mul_Word(X(X'First + j),
ffa_ch12_karatsub... 69 X(X'First - j + N),
ffa_ch12_karatsub... 70 Lo, Hi);
ffa_ch12_karatsub... 71 Accum;
ffa_ch12_karatsub... 72 Accum; -- Accum += 2 * (Hi:Lo)
ffa_ch12_karatsub... 73 end loop;
ffa_ch12_karatsub... 74 end SymmDigits;
ffa_ch12_karatsub... 75 pragma Inline_Always(SymmDigits);
ffa_ch12_karatsub... 76
ffa_ch12_karatsub... 77 procedure SqrDigit(N : in ColXX) is
ffa_ch12_karatsub... 78 begin
ffa_ch12_karatsub... 79 Sqr_Word(X(X'First + N), Lo, Hi); -- Hi:Lo := Square(N-th digit)
ffa_ch12_karatsub... 80 Accum; -- Accum += Hi:Lo
ffa_ch12_karatsub... 81 end SqrDigit;
ffa_ch12_karatsub... 82 pragma Inline_Always(SqrDigit);
ffa_ch12_karatsub... 83
ffa_ch12_karatsub... 84 procedure HaveDigit(N : in ColXX) is
ffa_ch12_karatsub... 85 begin
ffa_ch12_karatsub... 86 -- Save the Nth (indexed from zero) word of XX:
ffa_ch12_karatsub... 87 XX(XX'First + N) := A0;
ffa_ch12_karatsub... 88 -- Right-Shift the Accumulator by one Word width:
ffa_ch12_karatsub... 89 A0 := A1;
ffa_ch12_karatsub... 90 A1 := A2;
ffa_ch12_karatsub... 91 A2 := 0;
ffa_ch12_karatsub... 92 end HaveDigit;
ffa_ch12_karatsub... 93 pragma Inline_Always(HaveDigit);
ffa_ch12_karatsub... 94
ffa_ch12_karatsub... 95 -- Compute the Nth (indexed from zero) column of the Product
ffa_ch12_karatsub... 96 procedure Col(N : in ColXX; U : in ColXX; V : in ColXX) is
ffa_ch12_karatsub... 97 begin
ffa_ch12_karatsub... 98 -- The branch pattern depends only on FZ wordness
ffa_ch12_karatsub... 99 if N mod 2 = 0 then -- If we're doing an EVEN-numbered column:
ffa_ch12_karatsub... 100 SymmDigits(N, U, V - 1); -- Stop before V: it is the square case
ffa_ch12_karatsub... 101 SqrDigit(V); -- Compute the square case at V
ffa_ch12_karatsub... 102 else -- If we're doing an ODD-numbered column:
ffa_ch12_karatsub... 103 SymmDigits(N, U, V); -- All of them are the symmetric case
ffa_ch12_karatsub... 104 end if; -- After either even or odd column:
ffa_ch12_karatsub... 105 HaveDigit(N); -- We have the N-th digit of the result.
ffa_ch12_karatsub... 106 end Col;
ffa_ch12_karatsub... 107 pragma Inline_Always(Col);
ffa_ch12_karatsub... 108
ffa_ch12_karatsub... 109 begin
ffa_ch12_karatsub... 110 -- First col always exists:
ffa_ch12_karatsub... 111 SqrDigit(ColXX'First);
ffa_ch12_karatsub... 112 HaveDigit(ColXX'First);
ffa_ch12_karatsub... 113
ffa_ch12_karatsub... 114 -- Compute the lower half of the Product:
ffa_ch12_karatsub... 115 for i in 1 .. L - 1 loop
ffa_ch12_karatsub... 116 Col(i, 0, i / 2);
ffa_ch12_karatsub... 117 end loop;
ffa_ch12_karatsub... 118
ffa_ch12_karatsub... 119 -- Compute the upper half (sans last Word) of the Product:
ffa_ch12_karatsub... 120 for i in L .. LP - 2 loop
ffa_ch12_karatsub... 121 Col(i, i - L + 1, i / 2);
ffa_ch12_karatsub... 122 end loop;
ffa_ch12_karatsub... 123
ffa_ch12_karatsub... 124 -- The very last Word of the Product:
ffa_ch12_karatsub... 125 XX(XX'Last) := A0; -- Equiv. of XX(XX'First + 2*L - 1) := A0;
ffa_ch12_karatsub... 126 end FZ_Sqr_Comba;
ffa_ch12_karatsub... 127
ffa_ch12_karatsub... 128
ffa_ch12_karatsub... 129 -- Karatsuba's Squaring. (CAUTION: UNBUFFERED)
ffa_ch12_karatsub... 130 procedure Sqr_Karatsuba(X : in FZ;
ffa_ch12_karatsub... 131 XX : out FZ) is
ffa_ch12_karatsub... 132
ffa_ch12_karatsub... 133 -- L is the wordness of X. Guaranteed to be a power of two.
ffa_ch12_karatsub... 134 L : constant Word_Count := X'Length;
ffa_ch12_karatsub... 135
ffa_ch12_karatsub... 136 -- An 'LSeg' is the same length as X.
ffa_ch12_karatsub... 137 subtype LSeg is FZ(1 .. L);
ffa_ch12_karatsub... 138
ffa_ch12_karatsub... 139 -- K is HALF of the length of X.
ffa_ch12_karatsub... 140 K : constant Word_Index := L / 2;
ffa_ch12_karatsub... 141
ffa_ch12_karatsub... 142 -- A 'KSeg' is the same length as HALF of X.
ffa_ch12_karatsub... 143 subtype KSeg is FZ(1 .. K);
ffa_ch12_karatsub... 144
ffa_ch12_karatsub... 145 -- The three L-sized variables of the product equation, i.e.:
ffa_ch12_karatsub... 146 -- XX = LL + 2^(K*Bitness)(LL + HH - DD) + 2^(2*K*Bitness)HH
ffa_ch12_karatsub... 147 LL, DD, HH : LSeg;
ffa_ch12_karatsub... 148
ffa_ch12_karatsub... 149 -- K-sized term Dx, Dx := abs(XLo - XHi) to then get DD := Dx^2
ffa_ch12_karatsub... 150 Dx : KSeg;
ffa_ch12_karatsub... 151
ffa_ch12_karatsub... 152 -- IGNORED Subtraction borrow, sign of (XL - XH)
ffa_ch12_karatsub... 153 Cx : WBool; -- given as DD := Dx^2, DD is always positive
ffa_ch12_karatsub... 154 pragma Unreferenced(Cx);
ffa_ch12_karatsub... 155
ffa_ch12_karatsub... 156 -- Bottom and Top K-sized halves of X.
ffa_ch12_karatsub... 157 XLo : KSeg renames X( X'First .. X'Last - K );
ffa_ch12_karatsub... 158 XHi : KSeg renames X( X'First + K .. X'Last );
ffa_ch12_karatsub... 159
ffa_ch12_karatsub... 160 -- L-sized middle segment of the product XX (+/- K from the midpoint).
ffa_ch12_karatsub... 161 XXMid : LSeg renames XX( XX'First + K .. XX'Last - K );
ffa_ch12_karatsub... 162
ffa_ch12_karatsub... 163 -- Bottom and Top L-sized halves of the product XX.
ffa_ch12_karatsub... 164 XXLo : LSeg renames XX( XX'First .. XX'Last - L );
ffa_ch12_karatsub... 165 XXHi : LSeg renames XX( XX'First + L .. XX'Last );
ffa_ch12_karatsub... 166
ffa_ch12_karatsub... 167 -- Topmost K-sized quarter segment of the product XX, or 'tail'
ffa_ch12_karatsub... 168 XXHiHi : KSeg renames XXHi( XXHi'First + K .. XXHi'Last );
ffa_ch12_karatsub... 169
ffa_ch12_karatsub... 170 -- Carry from addition of HH and LL terms; borrow from subtraction of DD
ffa_ch12_karatsub... 171 C : WBool;
ffa_ch12_karatsub... 172
ffa_ch12_karatsub... 173 -- Barring a cosmic ray, 0 <= TC <= 2
ffa_ch12_karatsub... 174 subtype TailCarry is Word range 0 .. 2;
ffa_ch12_karatsub... 175
ffa_ch12_karatsub... 176 -- Tail-Carry accumulator, for the final ripple-out into XXHiHi
ffa_ch12_karatsub... 177 TC : TailCarry := 0;
ffa_ch12_karatsub... 178
ffa_ch12_karatsub... 179 -- Barring a cosmic ray, the tail ripple will NOT overflow.
ffa_ch12_karatsub... 180 FinalCarry : WZeroOrDie := 0;
ffa_ch12_karatsub... 181
ffa_ch12_karatsub... 182 begin
ffa_ch12_karatsub... 183
ffa_ch12_karatsub... 184 -- Recurse: LL := XLo^2
ffa_ch12_karatsub... 185 FZ_Square_Unbuffered(XLo, LL);
ffa_ch12_karatsub... 186
ffa_ch12_karatsub... 187 -- Recurse: HH := XHi^2
ffa_ch12_karatsub... 188 FZ_Square_Unbuffered(XHi, HH);
ffa_ch12_karatsub... 189
ffa_ch12_karatsub... 190 -- Dx := |XLo - XHi| , and we don't care about the borrow
ffa_ch12_karatsub... 191 FZ_Sub_Abs(X => XLo, Y => XHi, Difference => Dx, Underflow => Cx);
ffa_ch12_karatsub... 192
ffa_ch12_karatsub... 193 -- Recurse: DD := Dx^2 (always positive, which is why we don't need Cx)
ffa_ch12_karatsub... 194 FZ_Square_Unbuffered(Dx, DD);
ffa_ch12_karatsub... 195
ffa_ch12_karatsub... 196 -- XX := LL + 2^(2 * K * Bitness) * HH
ffa_ch12_karatsub... 197 XXLo := LL;
ffa_ch12_karatsub... 198 XXHi := HH;
ffa_ch12_karatsub... 199
ffa_ch12_karatsub... 200 -- XX += 2^(K * Bitness) * HH, carry goes into Tail Carry accumulator.
ffa_ch12_karatsub... 201 FZ_Add_D(X => XXMid, Y => HH, Overflow => TC);
ffa_ch12_karatsub... 202
ffa_ch12_karatsub... 203 -- XX += 2^(K * Bitness) * LL, ...
ffa_ch12_karatsub... 204 FZ_Add_D(X => XXMid, Y => LL, Overflow => C);
ffa_ch12_karatsub... 205
ffa_ch12_karatsub... 206 -- ... add the carry from LL addition into the Tail Carry accumulator.
ffa_ch12_karatsub... 207 TC := TC + C;
ffa_ch12_karatsub... 208
ffa_ch12_karatsub... 209 -- XX -= 2^(K * Bitness) * DD
ffa_ch12_karatsub... 210 FZ_Sub_D(X => XXMid, -- Because DD is always positive,
ffa_ch12_karatsub... 211 Y => DD, -- this term is always subtractive.
ffa_ch12_karatsub... 212 Underflow => C); -- C is the borrow from DD term subtraction
ffa_ch12_karatsub... 213
ffa_ch12_karatsub... 214 -- Get final Tail Carry for the ripple by subtracting DD term's borrow
ffa_ch12_karatsub... 215 TC := TC - C;
ffa_ch12_karatsub... 216
ffa_ch12_karatsub... 217 -- Ripple the Tail Carry into the tail.
ffa_ch12_karatsub... 218 FZ_Add_D_W(X => XXHiHi, W => TC, Overflow => FinalCarry);
ffa_ch12_karatsub... 219
ffa_ch12_karatsub... 220 end Sqr_Karatsuba;
ffa_ch12_karatsub... 221 -- CAUTION: Inlining prohibited for Sqr_Karatsuba !
ffa_ch12_karatsub... 222
ffa_ch12_karatsub... 223
ffa_ch12_karatsub... 224 -- Squaring. (CAUTION: UNBUFFERED)
ffa_ch12_karatsub... 225 procedure FZ_Square_Unbuffered(X : in FZ;
ffa_ch12_karatsub... 226 XX : out FZ) is
ffa_ch12_karatsub... 227 begin
ffa_ch12_karatsub... 228
ffa_ch12_karatsub... 229 if X'Length <= Sqr_Karatsuba_Thresh then
ffa_ch12_karatsub... 230
ffa_ch12_karatsub... 231 -- Base case:
ffa_ch12_karatsub... 232 FZ_Sqr_Comba(X, XX);
ffa_ch12_karatsub... 233
ffa_ch12_karatsub... 234 else
ffa_ch12_karatsub... 235
ffa_ch12_karatsub... 236 -- Recursive case:
ffa_ch12_karatsub... 237 Sqr_Karatsuba(X, XX);
ffa_ch12_karatsub... 238
ffa_ch12_karatsub... 239 end if;
ffa_ch12_karatsub... 240
ffa_ch12_karatsub... 241 end FZ_Square_Unbuffered;
ffa_ch12_karatsub... 242
ffa_ch12_karatsub... 243
ffa_ch12_karatsub... 244 -- Squaring. Preserves the inputs.
ffa_ch12_karatsub... 245 procedure FZ_Square_Buffered(X : in FZ;
ffa_ch12_karatsub... 246 XX_Lo : out FZ;
ffa_ch12_karatsub... 247 XX_Hi : out FZ) is
ffa_ch12_karatsub... 248
ffa_ch12_karatsub... 249 -- Product buffer.
ffa_ch12_karatsub... 250 P : FZ(1 .. 2 * X'Length);
ffa_ch12_karatsub... 251
ffa_ch12_karatsub... 252 begin
ffa_ch12_karatsub... 253
ffa_ch12_karatsub... 254 FZ_Square_Unbuffered(X, P);
ffa_ch12_karatsub... 255
ffa_ch12_karatsub... 256 XX_Lo := P(P'First .. P'First + X'Length - 1);
ffa_ch12_karatsub... 257 XX_Hi := P(P'First + X'Length .. P'Last);
ffa_ch12_karatsub... 258
ffa_ch12_karatsub... 259 end FZ_Square_Buffered;
ffa_ch12_karatsub... 260
ffa_ch12_karatsub... 261 end FZ_Sqr;