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