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;