------------------------------------------------------------------------------
------------------------------------------------------------------------------
-- This file is part of 'Finite Field Arithmetic', aka 'FFA'.               --
--                                                                          --
-- (C) 2019 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      --
-- distribute this software ; Should you use this software for any purpose, --
-- or copy and distribute it to anyone or in any manner, you are breaking   --
-- the laws of whatever soi-disant jurisdiction, and you promise to         --
-- continue doing so for the indefinite future. In any case, please         --
-- always : read and understand any software ; verify any PGP signatures    --
-- that you use - for any purpose.                                          --
--                                                                          --
-- See also http://trilema.com/2015/a-new-software-licensing-paradigm .     --
------------------------------------------------------------------------------

with Words;    use Words;
with FZ_Basic; use FZ_Basic;
with FZ_Shift; use FZ_Shift;
with FZ_QShft; use FZ_QShft;
with FZ_Arith; use FZ_Arith;
with FZ_Pred;  use FZ_Pred;


package body FZ_GCD is
   
   -- Find Greatest Common Divisor (GCD) of X and Y.
   -- Note that by convention, GCD(0, 0) = 0.
   procedure FZ_Greatest_Common_Divisor(X      : in  FZ;
                                        Y      : in  FZ;
                                        Result : out FZ) is
      
      -- Widths of X, Y, and Result are equal
      subtype Width is Word_Index range X'Range;
      
      -- Working buffers for GCD computation, initially equal to the inputs
      A      : FZ(Width) := X; -- GCD will appear in A in the end
      B      : FZ(Width) := Y;
      
      -- Evenness (negation of lowest bit) of A and B respectively
      Ae, Be : WBool;
      
      -- Common power-of-2 factor
      Twos   : Word := 0;
      
      -- |A - B|
      D      : FZ(Width);
      
      -- This flag is set iff A < B
      A_lt_B : WBool;
      
   begin
      
      -- For convergence, requires number of shots equal to 2 * FZ_Bitness:
      for i in 1 .. 2 * FZ_Bitness(X) loop
         
         -- If A is even, A := A >> 1; otherwise A := A
         Ae := 1 - FZ_OddP(A);
         FZ_ShiftRight(A, A, WBit_Index(Ae));
         
         -- If B is even, B := B >> 1; otherwise B := B
         Be := 1 - FZ_OddP(B);
         FZ_ShiftRight(B, B, WBit_Index(Be));
         
         -- If both A and B were even, increment the common power-of-two
         Twos := Twos + (Ae and Be);
         
         -- D := |A - B|
         FZ_Sub_Abs(X => A, Y => B, Difference => D, Underflow => A_lt_B);
         
         -- B' := min(A, B)
         FZ_Mux(X => B, Y => A, Result => B, Sel => A_lt_B);
         
         -- A' := |A - B|
         A := D;
         
      end loop;
      
      -- Reintroduce the common power-of-2 factor stored in 'Twos'
      FZ_Quiet_ShiftLeft(N => A, ShiftedN => A, Count => Indices(Twos));
      
      -- Output final result
      Result := A;
      
   end FZ_Greatest_Common_Divisor;
   
end FZ_GCD;