Finite Field Arithmetic

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)
 251       NoCarry : WZeroOrDie := 0;
 252       
 253       -- Borrow from Subtraction in (5) is meaningless, and is discarded
 254       IgnoreC : WBool;
 255       pragma Unreferenced(IgnoreC);
 256       
 257    begin
 258       
 259       -- Result is initially zero (and will stay zero if Modulus = 1)
 260       FZ_Clear(XReduced);
 261       
 262       -- (1) Xs := X >> Jm
 263       FZ_Quiet_ShiftRight(N => X, ShiftedN => Xs, Count => Bar.J);
 264       
 265       -- (2) Z  := Xs * Bm
 266       FZ_Multiply_Unbuffered(X => Bar.B, Y => Xs, XY => Z);
 267       
 268       -- (3) Zs := Z >> 2Wm - Jm (already thrown lower Wm, so only Wm - Jm now)
 269       FZ_Quiet_ShiftRight(N => ZHi, ShiftedN => ZHi, Count => Bar.ZSlide);
 270       
 271       -- (4) Q  := Zs * M (this is split into three operations, see below)
 272       
 273       -- ... first, Q := ZsLo * M,
 274       FZ_Multiply_Unbuffered(ZsLo, Modulus, Q);
 275       
 276       -- ... then, compute ZsHiM := ZsHi * M,
 277       FZ_Low_Multiply_Unbuffered(ZsHi, Modulus, ZsHiM);
 278       
 279       -- ... finally, add ZsHiM to upper half of Q.
 280       FZ_Add_D(X => QHi, Y => ZsHiM, Overflow => NoCarry);
 281       
 282       -- (5) R  := X - Q (we only need Rl-sized segments of X and Q here)
 283       FZ_Sub(X => X(1 .. Rl), Y => Q(1 .. Rl),
 284              Difference => R, Underflow => IgnoreC); -- Borrow is discarded
 285       
 286       -- (6) S1 := R - M, C1 := Borrow (1st gated subtraction of Modulus)
 287       FZ_Sub(X => R, Y => Bar.ZXM, Difference => S, Underflow => C);
 288       
 289       -- (7) R := {C1=0 -> S1, C1=1 -> R}
 290       FZ_Mux(X => S, Y => R, Result => R, Sel => C);
 291       
 292       -- (8) S2 := R - M, C2 := Borrow (2nd gated subtraction of Modulus)
 293       FZ_Sub(X => R, Y => Bar.ZXM, Difference => S, Underflow => C);
 294       
 295       -- (9) R := {C2=0 -> S2, C2=1 -> R}
 296       FZ_Mux(X => S, Y => R, Result => R, Sel => C);
 297       
 298       -- (10) RFinal := {DM=0 -> R, DM=1 -> 0} (handle the degenerate case)
 299       FZ_Mux(X => R(1 .. Ml), Y => XReduced, Result => XReduced,
 300              Sel => Bar.Degenerate); -- If Modulus = 1, then XReduced is 0.
 301       
 302    end FZ_Barrett_Reduce;
 303    
 304 end FZ_Barr;