(2 . 7)(2 . 7)
------------------------------------------------------------------------------
-- This file is part of 'Finite Field Arithmetic', aka 'FFA'.               --
--                                                                          --
-- (C) 2017 Stanislav Datskovskiy ( www.loper-os.org )                      --
-- (C) 2018 Stanislav Datskovskiy ( www.loper-os.org )                      --
-- http://wot.deedbot.org/17215D118B7239507FAFED98B98228A001ABFFC7.html     --
--                                                                          --
-- You do not have, nor can you ever acquire the right to use, copy or      --
(17 . 17)(17 . 18)
------------------------------------------------------------------------------
------------------------------------------------------------------------------

with Words;    use Words;
with W_Shifts; use W_Shifts;
with FZ_Basic; use FZ_Basic;
with FZ_Pred;  use FZ_Pred;
with FZ_Shift; use FZ_Shift;
with FZ_Mul;   use FZ_Mul;
with FZ_Sqr;   use FZ_Sqr;
with FZ_Divis; use FZ_Divis;
with FZ_Barr;  use FZ_Barr;


package body FZ_ModEx is
   
   -- Modular Multiply: Product := X*Y mod Modulus
   -- (Conventional) Modular Multiply: Product := X*Y mod Modulus
   procedure FZ_Mod_Mul(X        : in  FZ;
                        Y        : in  FZ;
                        Modulus  : in  FZ;
(54 . 7)(55 . 7)
   end FZ_Mod_Mul;
   
   
   -- Modular Squaring: Product := X*X mod Modulus
   -- (Conventional) Modular Squaring: Product := X*X mod Modulus
   procedure FZ_Mod_Sqr(X        : in  FZ;
                        Modulus  : in  FZ;
                        Product  : out FZ) is
(80 . 44)(81 . 68)
   end FZ_Mod_Sqr;
   
   
   -- Modular Exponent: Result := Base^Exponent mod Modulus
   -- (Barrettronic) Modular Exponent: Result := Base^Exponent mod Modulus
   procedure FZ_Mod_Exp(Base     : in  FZ;
                        Exponent : in  FZ;
                        Modulus  : in  FZ;
                        Result   : out FZ) is
      
      -- Working register for the squaring; initially is copy of Base
      B : FZ(Base'Range)     := Base;
      -- Double-width scratch buffer for the modular operations
      D   : FZ(1 .. Base'Length * 2);
      
      -- Copy of Exponent, for cycling through its bits
      E : FZ(Exponent'Range) := Exponent;
      -- Working register for the squaring; initially is copy of Base
      B   : FZ(Base'Range) := Base;
      
      -- Register for the Mux operation
      T : FZ(Result'Range);
      T   : FZ(Result'Range);
      
      -- Buffer register for the Result
      R : FZ(Result'Range);
      R   : FZ(Result'Range);
      
      -- Space for Barrettoid
      Bar : Barretoid(ZXMLength       => Modulus'Length + 1,
                      BarretoidLength => 2 * B'Length);
      
   begin
      
      -- First, pre-compute the Barretoid for the given Modulus:
      FZ_Make_Barrettoid(Modulus => Modulus, Result => Bar);
      
      -- Result := 1
      WBool_To_FZ(1, R);
      
      -- For each bit of R width:
      for i in 1 .. FZ_Bitness(R) loop
         
         -- T := Result * B mod Modulus
         FZ_Mod_Mul(X => R, Y => B, Modulus => Modulus, Product => T);
         
         -- Sel is the current low bit of E;
         --    When Sel=0 -> Result := Result;
         --    When Sel=1 -> Result := T
         FZ_Mux(X => R, Y => T, Result => R, Sel => FZ_OddP(E));
      -- For each Word of the Exponent:
      for i in Exponent'Range loop
         
         -- Advance to the next bit of E
         FZ_ShiftRight(E, E, 1);
         declare
            
            -- The current Word of the Exponent
            Wi : Word := Exponent(i);
            
         begin
            
            -- For each bit of Wi:
            for j in 1 .. Bitness loop
               
               -- T := Result * B mod Modulus
               FZ_Multiply_Unbuffered(X => R, Y => B, XY => D);
               FZ_Barrett_Reduce(X => D, Bar => Bar, XReduced => T);
               
               -- Sel is the current bit of Exponent;
               --    When Sel=0 -> Result := Result;
               --    When Sel=1 -> Result := T
               FZ_Mux(X => R, Y => T, Result => R, Sel => Wi and 1);
               
               -- Advance to the next bit of Wi (i.e. next bit of Exponent)
               Wi := Shift_Right(Wi, 1);
               
               -- B := B^2 mod Modulus
               FZ_Square_Unbuffered(X => B, XX => D);
               FZ_Barrett_Reduce(X => D, Bar => Bar, XReduced => B);
               
            end loop;
         
         -- B := B^2 mod Modulus
         FZ_Mod_Sqr(X => B, Modulus => Modulus, Product => B);
         end;
         
      end loop;