raw
ffa_ch16_miller_r...    1 ------------------------------------------------------------------------------
ffa_ch16_miller_r...    2 ------------------------------------------------------------------------------
ffa_ch16_miller_r...    3 -- This file is part of 'Finite Field Arithmetic', aka 'FFA'.               --
ffa_ch16_miller_r...    4 --                                                                          --
ffa_ch16_miller_r...    5 -- (C) 2019 Stanislav Datskovskiy ( www.loper-os.org )                      --
ffa_ch16_miller_r...    6 -- http://wot.deedbot.org/17215D118B7239507FAFED98B98228A001ABFFC7.html     --
ffa_ch16_miller_r...    7 --                                                                          --
ffa_ch16_miller_r...    8 -- You do not have, nor can you ever acquire the right to use, copy or      --
ffa_ch16_miller_r...    9 -- distribute this software ; Should you use this software for any purpose, --
ffa_ch16_miller_r...   10 -- or copy and distribute it to anyone or in any manner, you are breaking   --
ffa_ch16_miller_r...   11 -- the laws of whatever soi-disant jurisdiction, and you promise to         --
ffa_ch16_miller_r...   12 -- continue doing so for the indefinite future. In any case, please         --
ffa_ch16_miller_r...   13 -- always : read and understand any software ; verify any PGP signatures    --
ffa_ch16_miller_r...   14 -- that you use - for any purpose.                                          --
ffa_ch16_miller_r...   15 --                                                                          --
ffa_ch16_miller_r...   16 -- See also http://trilema.com/2015/a-new-software-licensing-paradigm .     --
ffa_ch16_miller_r...   17 ------------------------------------------------------------------------------
ffa_ch16_miller_r...   18 
ffa_ch16_miller_r...   19 with Word_Ops; use Word_Ops;
ffa_ch16_miller_r...   20 with W_Pred;   use W_Pred;
ffa_ch16_miller_r...   21 with W_Shifts; use W_Shifts;
ffa_ch16_miller_r...   22 
ffa_ch16_miller_r...   23 with FZ_Basic; use FZ_Basic;
ffa_ch16_miller_r...   24 with FZ_QShft; use FZ_QShft;
ffa_ch16_miller_r...   25 with FZ_Arith; use FZ_Arith;
ffa_ch16_miller_r...   26 with FZ_Divis; use FZ_Divis;
ffa_ch16_miller_r...   27 with FZ_Pred;  use FZ_Pred;
ffa_ch16_miller_r...   28 with FZ_Cmp;   use FZ_Cmp;
ffa_ch16_miller_r...   29 with FZ_Barr;  use FZ_Barr;
ffa_ch16_miller_r...   30 with FZ_ModEx; use FZ_ModEx;
ffa_ch16_miller_r...   31 
ffa_ch16_miller_r...   32 
ffa_ch16_miller_r...   33 package body FZ_Prime is
ffa_ch16_miller_r...   34    
ffa_ch16_miller_r...   35    -- Find the highest power of 2 which divides N. ( or 0 if N is zero. )
ffa_ch16_miller_r...   36    function FZ_Count_Bottom_Zeros(N : in FZ) return FZBit_Index is
ffa_ch16_miller_r...   37       
ffa_ch16_miller_r...   38       -- The result (default : 0, will remain 0 if N is in fact zero)
ffa_ch16_miller_r...   39       Index     : Word := 0;
ffa_ch16_miller_r...   40       
ffa_ch16_miller_r...   41       -- The currently-examined Word, starting from the highest
ffa_ch16_miller_r...   42       Ni        : Word;
ffa_ch16_miller_r...   43       
ffa_ch16_miller_r...   44       -- The most recently-seen nonzero Word, if indeed any exist
ffa_ch16_miller_r...   45       W         : Word := 0;
ffa_ch16_miller_r...   46       
ffa_ch16_miller_r...   47       -- 1 if currently-examined Word is zero, otherwise 0
ffa_ch16_miller_r...   48       NiZ       : WBool;
ffa_ch16_miller_r...   49       
ffa_ch16_miller_r...   50    begin
ffa_ch16_miller_r...   51       
ffa_ch16_miller_r...   52       -- Find lowest non-zero Word (or simply stop at lowest, if N = 0)
ffa_ch16_miller_r...   53       for i in reverse 0 .. Indices(N'Length - 1) loop
ffa_ch16_miller_r...   54          Ni    := N(N'First + i);              -- Ni := current Word;
ffa_ch16_miller_r...   55          NiZ   := W_ZeroP(Ni);                 -- ... is this Word = 0?
ffa_ch16_miller_r...   56          Index := W_Mux(Word(i), Index, NiZ);  -- If NO, save its Index,
ffa_ch16_miller_r...   57          W     := W_Mux(Ni,      W,     NiZ);  -- ... and save the Word.
ffa_ch16_miller_r...   58       end loop;
ffa_ch16_miller_r...   59       
ffa_ch16_miller_r...   60       -- Set Index to be the bit-position of the lowest non-zero Word:
ffa_ch16_miller_r...   61       Index := Shift_Left(Index, BitnessLog2); -- Index := Index * Bitness
ffa_ch16_miller_r...   62       
ffa_ch16_miller_r...   63       -- Handle degenerate case: if W is equal to 0, Index is not changed;
ffa_ch16_miller_r...   64       -- Otherwise, start by advancing Index by an entire Word Bitness:
ffa_ch16_miller_r...   65       Index := Index + ((0 - W_NZeroP(W)) and Word(Bitness));
ffa_ch16_miller_r...   66       
ffa_ch16_miller_r...   67       -- Now crank back the Index by the number of high bits of W (i.e.
ffa_ch16_miller_r...   68       -- starting from its top) that must be discarded for W to become zero:
ffa_ch16_miller_r...   69       for b in 1 .. Bitness loop
ffa_ch16_miller_r...   70          
ffa_ch16_miller_r...   71          -- If W is non-zero at this time, decrement the Index:
ffa_ch16_miller_r...   72          Index := Index - W_NZeroP(W);
ffa_ch16_miller_r...   73          
ffa_ch16_miller_r...   74          -- Advance W left, i.e. discard the top bit of it:
ffa_ch16_miller_r...   75          W     := Shift_Left(W, 1);
ffa_ch16_miller_r...   76          
ffa_ch16_miller_r...   77       end loop;
ffa_ch16_miller_r...   78       
ffa_ch16_miller_r...   79       -- If N = 0, result will be 0; otherwise: length of bottom zeros in N.
ffa_ch16_miller_r...   80       return FZBit_Index(Index);
ffa_ch16_miller_r...   81       
ffa_ch16_miller_r...   82    end FZ_Count_Bottom_Zeros;
ffa_ch16_miller_r...   83    
ffa_ch16_miller_r...   84    
ffa_ch16_miller_r...   85    -- Constant-Time Miller-Rabin Test on N using the given Witness.
ffa_ch16_miller_r...   86    -- Witness will be used as-is if it conforms to the valid bounds,
ffa_ch16_miller_r...   87    -- i.e. 2 <= Witness <= N - 2; otherwise will be transformed into a
ffa_ch16_miller_r...   88    -- valid Witness via modular arithmetic.
ffa_ch16_miller_r...   89    -- Outputs ONE if N WAS FOUND composite; ZERO if NOT FOUND.
ffa_ch16_miller_r...   90    -- Handles degenerate cases of N that M-R per se cannot eat:
ffa_ch16_miller_r...   91    -- N=0, N=1: ALWAYS 'FOUND COMPOSITE'; 2, 3 - ALWAYS 'NOT FOUND'.
ffa_ch16_miller_r...   92    -- If N is Even and not equal to 2, N is ALWAYS 'FOUND COMPOSITE.'
ffa_ch16_miller_r...   93    -- For ALL other N, the output is equal to that of the M-R test.
ffa_ch16_miller_r...   94    function FZ_MR_Composite_On_Witness(N       : in  FZ;
ffa_ch16_miller_r...   95                                        Witness : in  FZ) return WBool is
ffa_ch16_miller_r...   96       
ffa_ch16_miller_r...   97       -- Widths of N, Witness, and Result are equal :
ffa_ch16_miller_r...   98       subtype Width is Word_Index range N'Range;
ffa_ch16_miller_r...   99       
ffa_ch16_miller_r...  100       -- Whether N is 'small' (i.e. 1 Word) and therefore possibly degenerate :
ffa_ch16_miller_r...  101       N_Is_Small       : constant WBool := FZ_OneWordP(N);
ffa_ch16_miller_r...  102       
ffa_ch16_miller_r...  103       -- Head of N, for detecting (and handling) the degenerate-N case :
ffa_ch16_miller_r...  104       N_Head           : constant Word  := FZ_Get_Head(N);
ffa_ch16_miller_r...  105       
ffa_ch16_miller_r...  106       -- Zero and One are defined as COMPOSITE degenerate cases of N :
ffa_ch16_miller_r...  107       N_Is_Zero_Or_One : constant WBool
ffa_ch16_miller_r...  108         := N_Is_Small and (W_EqP(N_Head, 0) or W_EqP(N_Head, 1));
ffa_ch16_miller_r...  109       
ffa_ch16_miller_r...  110       -- Two and Three are PRIME degenerate cases of N :
ffa_ch16_miller_r...  111       N_Is_Two         : constant WBool := N_Is_Small and W_EqP(N_Head, 2);
ffa_ch16_miller_r...  112       N_Is_Three       : constant WBool := N_Is_Small and W_EqP(N_Head, 3);
ffa_ch16_miller_r...  113       
ffa_ch16_miller_r...  114       -- The COMPOSITE degenerate case of N where N != 2 and N is Even :
ffa_ch16_miller_r...  115       N_Ne_2_Is_Even   : constant WBool :=
ffa_ch16_miller_r...  116         (1 - N_Is_Two) and (1 - W_OddP(N_Head));
ffa_ch16_miller_r...  117       
ffa_ch16_miller_r...  118       -- Degeneracy latch. If N is Zero or One, then set immediately :
ffa_ch16_miller_r...  119       Degen_Composite  : constant WBool := N_Is_Zero_Or_One or N_Ne_2_Is_Even;
ffa_ch16_miller_r...  120       
ffa_ch16_miller_r...  121       -- Possible-primality latch. If N is Two or Three, then set immediately :
ffa_ch16_miller_r...  122       Possibly_Prime   : WBool := N_Is_Two or N_Is_Three;
ffa_ch16_miller_r...  123       
ffa_ch16_miller_r...  124       -- The writable copy of N that we will operate on :
ffa_ch16_miller_r...  125       X                : FZ(Width) := N;
ffa_ch16_miller_r...  126       
ffa_ch16_miller_r...  127       -- Potentially-replaced head of X (if degenerate N) :
ffa_ch16_miller_r...  128       X_Head           : Word := N_Head;
ffa_ch16_miller_r...  129       
ffa_ch16_miller_r...  130       -- Will be Barrettoid(X), for operations modulo X :
ffa_ch16_miller_r...  131       XBar             : Barretoid(ZXMLength       => X'Length + 1,
ffa_ch16_miller_r...  132                                    BarretoidLength => 2 * X'Length);
ffa_ch16_miller_r...  133       
ffa_ch16_miller_r...  134       -- The Bound Witness, constrained to valid range 2 <= BW <= X - 2 :
ffa_ch16_miller_r...  135       BW               : FZ(Width);
ffa_ch16_miller_r...  136       
ffa_ch16_miller_r...  137       -- Head of BW, for detecting crossing of the lower bound :
ffa_ch16_miller_r...  138       BW_Head          : Word;
ffa_ch16_miller_r...  139       
ffa_ch16_miller_r...  140       -- X - 1 (for M-R proper, and to constrain BW) :
ffa_ch16_miller_r...  141       X_Minus_One      : FZ(Width);
ffa_ch16_miller_r...  142       
ffa_ch16_miller_r...  143       -- X - 1 = 2^R * K, with K odd :
ffa_ch16_miller_r...  144       K                : FZ(Width);
ffa_ch16_miller_r...  145       R                : FZBit_Index;
ffa_ch16_miller_r...  146       
ffa_ch16_miller_r...  147       -- Working register for all M-R modular operations :
ffa_ch16_miller_r...  148       T                : FZ(Width);
ffa_ch16_miller_r...  149       
ffa_ch16_miller_r...  150       -- For subtraction where no underflow can happen, barring cosmic ray:
ffa_ch16_miller_r...  151       NoCarry          : WZeroOrDie := 0;
ffa_ch16_miller_r...  152       
ffa_ch16_miller_r...  153    begin
ffa_ch16_miller_r...  154       
ffa_ch16_miller_r...  155       -- First: X, which will remain equal to N unless N is degenerate:
ffa_ch16_miller_r...  156       
ffa_ch16_miller_r...  157       -- If N is one of the two prohibited small primes (2,3) X will become 5:
ffa_ch16_miller_r...  158       X_Head := W_Mux(A => X_Head, B => 5, Sel => Possibly_Prime);
ffa_ch16_miller_r...  159       
ffa_ch16_miller_r...  160       -- If one of the two prohibited small composites (0,1), or even, then 9:
ffa_ch16_miller_r...  161       X_Head := W_Mux(A => X_Head, B => 9, Sel => Degen_Composite);
ffa_ch16_miller_r...  162       
ffa_ch16_miller_r...  163       -- Given as we're stuck carrying out M-R even if N is degenerate case,
ffa_ch16_miller_r...  164       -- then let the M-R do The Right Thing, for added cosmic ray resistance.
ffa_ch16_miller_r...  165       
ffa_ch16_miller_r...  166       -- X gets a new head, if N was a degenerate case; else keeps old head:
ffa_ch16_miller_r...  167       FZ_Set_Head(X, X_Head);
ffa_ch16_miller_r...  168       
ffa_ch16_miller_r...  169       -- Compute X - 1. As now X > 3, underflow is not permitted:
ffa_ch16_miller_r...  170       FZ_Sub_W(X => X, W => 1, Difference => X_Minus_One,
ffa_ch16_miller_r...  171                Underflow => NoCarry);
ffa_ch16_miller_r...  172       
ffa_ch16_miller_r...  173       -- Now, compute BW (Bound Witness), which satisfies 2 <= BW <= X - 2:
ffa_ch16_miller_r...  174       
ffa_ch16_miller_r...  175       -- First, BW := Witness mod (X - 1). After this, 0 <= BW <= X - 2:
ffa_ch16_miller_r...  176       FZ_Mod(Dividend => Witness, Divisor => X_Minus_One, Remainder => BW);
ffa_ch16_miller_r...  177       
ffa_ch16_miller_r...  178       -- Now, adjust BW if it is found to be below Two (the lower bound) :
ffa_ch16_miller_r...  179       
ffa_ch16_miller_r...  180       -- Get head of BW:
ffa_ch16_miller_r...  181       BW_Head := FZ_Get_Head(BW);
ffa_ch16_miller_r...  182       
ffa_ch16_miller_r...  183       -- If BW is equal to Zero or One, it will be forced to equal Two:
ffa_ch16_miller_r...  184       BW_Head := W_Mux(A   => BW_Head,
ffa_ch16_miller_r...  185                        B   => 2,
ffa_ch16_miller_r...  186                        Sel => FZ_OneWordP(BW) and
ffa_ch16_miller_r...  187                          (W_EqP(BW_Head, 0) or W_EqP(BW_Head, 1)));
ffa_ch16_miller_r...  188       
ffa_ch16_miller_r...  189       -- BW gets the new head, if it must; otherwise keeps its old head:
ffa_ch16_miller_r...  190       FZ_Set_Head(BW, BW_Head);
ffa_ch16_miller_r...  191       
ffa_ch16_miller_r...  192       -- We finished adjusting X and BW for degeneracies, and can now M-R:
ffa_ch16_miller_r...  193       
ffa_ch16_miller_r...  194       -- Generate Barrettoid(X) to use in all of the modulo-X operations:
ffa_ch16_miller_r...  195       FZ_Make_Barrettoid(Modulus => X, Result => XBar);
ffa_ch16_miller_r...  196       
ffa_ch17_peh.kv       197       -- Find R >= 1, and odd K, where X - 1 = 2^R * K :
ffa_ch16_miller_r...  198       
ffa_ch16_miller_r...  199       -- ... first, find R, the largest power of two which divides X - 1 :
ffa_ch16_miller_r...  200       R := FZ_Count_Bottom_Zeros(X_Minus_One);
ffa_ch16_miller_r...  201       
ffa_ch16_miller_r...  202       -- ... and now obtain K := X / 2^R, i.e., K := X >> R :
ffa_ch16_miller_r...  203       FZ_Quiet_ShiftRight(N => X_Minus_One, ShiftedN => K, Count => R);
ffa_ch16_miller_r...  204       
ffa_ch16_miller_r...  205       -- T := BW^K mod X
ffa_ch16_miller_r...  206       FZ_Mod_Exp_Barrett(Base => BW, Exponent => K, Bar => XBar, Result => T);
ffa_ch16_miller_r...  207       
ffa_ch16_miller_r...  208       -- if T = 1 or T = X - 1, then X is possibly-prime:
ffa_ch16_miller_r...  209       Possibly_Prime := Possibly_Prime or
ffa_ch16_miller_r...  210         FZ_EqP_W(T, 1) or FZ_EqP(T, X_Minus_One);
ffa_ch16_miller_r...  211       
ffa_ch16_miller_r...  212       -- Needs R - 1 times, but perform for max possible count (gated):
ffa_ch16_miller_r...  213       for i in 1 .. FZ_Bitness(N) - 1 loop
ffa_ch16_miller_r...  214          
ffa_ch16_miller_r...  215          -- T := T^2 mod X
ffa_ch16_miller_r...  216          FZ_Mod_Sqr_Barrett(X => T, Bar => XBar, Product => T);
ffa_ch16_miller_r...  217          
ffa_ch16_miller_r...  218          -- if T = X - 1, and i < R, then X is possibly-prime:
ffa_ch16_miller_r...  219          Possibly_Prime := Possibly_Prime or 
ffa_ch16_miller_r...  220            (FZ_EqP(T, X_Minus_One) and W_LtP(Word(i), Word(R)));
ffa_ch16_miller_r...  221          
ffa_ch16_miller_r...  222       end loop;
ffa_ch16_miller_r...  223       
ffa_ch16_miller_r...  224       -- The output, which will be 1 iff X WAS FOUND to be composite via BW,
ffa_ch16_miller_r...  225       -- ... or if X was found to equal Zero or One (without regard to BW.)
ffa_ch16_miller_r...  226       return (1 - Possibly_Prime) or Degen_Composite;
ffa_ch16_miller_r...  227       -- If X was found to equal Two or Three, output will be 0 (under any BW).
ffa_ch16_miller_r...  228       
ffa_ch16_miller_r...  229    end FZ_MR_Composite_On_Witness;
ffa_ch16_miller_r...  230    
ffa_ch16_miller_r...  231 end FZ_Prime;