-- Perform a Miller-Rabin Primality Test on N using the given Witness.
   -- Result will equal 1 if N WAS FOUND to be composite; if NOT FOUND, 0.
   -- This is a piecewise-defined function. Always in constant spacetime:
   --  0, 1: always 'FOUND composite'; 2, 3 - always 'NOT found composite.'
   --  If N is Even and not equal to 2: always 'FOUND composite.'
   --  For ALL other N: the result of the M-R test will be the output.
   procedure FZ_Miller_Rabin_Test(N       : in  FZ;
                                  Witness : in  FZ;
                                  Result  : out FZ) is
      
      -- Widths of N, Witness, and Result are equal :
      subtype Width is Word_Index range N'Range;
      
      -- Whether N is 'small' (i.e. 1 Word) and therefore possibly degenerate :
      N_Is_Small      : constant WBool := FZ_OneWordP(N);
      N_Head          : constant Word  := FZ_Get_Head(N);
      
      -- Zero and One are defined as COMPOSITE degenerate cases of N :
      N_Is_Zero_Or_One : constant WBool
        := N_Is_Small and (W_EqP(N_Head, 0) or W_EqP(N_Head, 1));
      
      -- Two and Three are PRIME degenerate cases of N :
      N_Is_Two        : constant WBool := N_Is_Small and W_EqP(N_Head, 2);
      N_Is_Three      : constant WBool := N_Is_Small and W_EqP(N_Head, 3);
      
      -- The remaining degenerate case of N, where N != 2 and N is Even :
      N_Ne_2_Is_Even  : constant WBool := (1 - N_Is_Two) and (1 - FZ_OddP(N));
      
      -- Degeneracy latch. If N is Zero or One, then set immediately :
      Degen_Composite : constant WBool := N_Is_Zero_Or_One or N_Ne_2_Is_Even;
      
      -- Possible-primality latch. If N is Two or Three, then set immediately :
      Possibly_Prime  : WBool := N_Is_Two or N_Is_Three;
      
      -- The writable copy of N that we will operate on :
      X               : FZ(Width) := N;
      
      -- Potentially-replaced head of X (if degenerate N) :
      X_Head          : Word := N_Head;
      
      -- Space for Barrettoid, for operations modulo X
      Bar             : Barretoid(ZXMLength       => X'Length + 1,
                                  BarretoidLength => 2 * X'Length);
      
      -- The Bound Witness, constrained to range {2 .. X - 2} :
      BW              : FZ(Width);
      
      -- Head of BW, for detecting crossing of the lower bound :
      BW_Head         : Word;
      
      -- X - 1 (for M-R proper) :
      X_Minus_One     : FZ(Width);
      
      -- X - 2 (for upper bound of the Bound Witness) :
      X_Minus_Two     : FZ(Width);
      
      -- N - 1 = 2^R * K, with K odd :
      K               : FZ(Width);
      R               : FZBit_Index;
      
      -- Working register :
      T               : FZ(Width);
      
      -- For subtractions where no underflow can happen, barring cosmic ray:
      NoCarry         : WZeroOrDie := 0;
      
   begin
      
      -- If N is one of the two prohibited small primes (2,3) X will become 5:
      X_Head := W_Mux(A => X_Head, B => 5, Sel => Possibly_Prime);
      
      -- If N is one of the two prohibited small composites (0,1), then 4:
      X_Head := W_Mux(A => X_Head, B => 4, Sel => N_Is_Zero_Or_One);
      
      -- X gets a new head, if N was a degenerate case; else keeps old head:
      FZ_Set_Head(X, X_Head);
      
      -- Compute X - 1. As now X > 3, underflow is not permitted:
      FZ_Sub_W(X => X, W => 1, Difference => X_Minus_One,
               Underflow => NoCarry);
      
      -- Now, compute the Bound Witness, which satisfies 2 <= BW <= X - 2:
      
      -- Compute X - 2, the upper bound for the Witness:
      FZ_Sub_W(X => X, W => 2, Difference => X_Minus_Two,
               Underflow => NoCarry);
      
      -- First, BW := Witness mod (X - 2) :
      FZ_Mod(Dividend => Witness, Divisor => X_Minus_Two, Remainder => BW);
      
      -- Now, adjust BW if it were below the lower bound (2) :
      BW_Head := FZ_Get_Head(BW);
      
      -- If BW is now Zero or One, it will be forced to equal Two:
      BW_Head := W_Mux(A   => BW_Head,
                       B   => 2,
                       Sel => FZ_OneWordP(BW) and
                         (W_EqP(BW_Head, 0) or W_EqP(BW_Head, 1)));
      
      -- BW gets a new head, if it must; otherwise keeps its old head:
      FZ_Set_Head(BW, BW_Head);
      
      -- Now we are done with adjusting X and BW for degeneracies, can M-R:
      
      -- Generate Barretoid(X) to use for all of the modulo-X operations:
      FZ_Make_Barrettoid(Modulus => X, Result => Bar);
      
      -- Find R >= 1, and odd K, where X − 1 = 2^R * K:
      -- ... first, find R, the largest power of two which divides X - 1:
      R := FZ_Count_Bottom_Zeros(X_Minus_One);
      
      -- ... K := X / 2^R.
      FZ_Quiet_ShiftRight(N => X_Minus_One, ShiftedN => K, Count => R);
      
      -- T := W^K mod X
      FZ_Mod_Exp_Barrett(Base => BW, Exponent => K, Bar => Bar, Result => T);
      
      -- if T = 1 or T = X - 1, then possibly-prime:
      Possibly_Prime := Possibly_Prime or FZ_EqP_W(T, 1);
      Possibly_Prime := Possibly_Prime or FZ_EqP(T, X_Minus_One);
      
      -- Needs R - 1 times, but perform for max possible count (will cycle):
      for I in 1 .. FZ_Bitness(N) - 1 loop
         
         -- T := T^2 mod X
         FZ_Mod_Sqr_Barrett(X => T, Bar => Bar, Product => T);
         
         -- if T = X - 1, then possibly-prime:
         Possibly_Prime := Possibly_Prime or FZ_EqP(T, X_Minus_One);
         
      end loop;
      
      -- Output the Result:
      WBool_To_FZ((1 - Possibly_Prime) or Degen_Composite, Result);
      
   end FZ_Miller_Rabin_Test;