-- Find Greatest Common Divisor (GCD) of X and Y.
   -- It is also defined that GCD(N, 0) = GCD(0, N) = GCD(0, 0) = 0.
   procedure FZ_Greatest_Common_Divisor(X      : in  FZ;
                                        Y      : in  FZ;
                                        Result : out FZ) is
      
      -- Width of X, Y, Result
      subtype Width is Word_Index range X'Range;
      
      -- Result (initially 0, and in degenerate case, stays 0)
      GCD    : FZ(Width) := (others => 0);
      
      -- Working buffers
      A      : FZ(Width) := X;
      B      : FZ(Width) := Y;
      
      -- Eveneties of A and B
      Ea, Eb : WBool;
      
      -- Shared power-of-2 factor
      Twos   : Word := 0;
      
      -- |A - B|
      AsB    : FZ(Width);
      
      -- Whether A < B
      AltB   : WBool;
      
      -- Whether we have the degenerate case
      Degen  : constant WBool := FZ_ZeroP(A) or FZ_ZeroP(B);
      
   begin
      
      for i in 1 .. 2 * FZ_Bitness(X) loop
         
         -- Record whether A, B are even:
         Ea := 1 - FZ_OddP(A);
         Eb := 1 - FZ_OddP(B);
         
         -- If A is even, A := A >> 1
         FZ_ShiftRight(A, A, WBit_Index(Ea));
         
         -- If B is even, B := B >> 1
         FZ_ShiftRight(B, B, WBit_Index(Eb));
         
         -- If both A and B were even, increment the Shared Twos counter
         Twos := Twos + (Ea and Eb);
         
         -- AsB := |A - B|
         FZ_Sub_Abs(X => A, Y => B, Difference => AsB, Underflow => AltB);
         
         -- B' := min(A, B)
         FZ_Mux(X => B,   Y => A,   Result => B, Sel => AltB);
         
         -- A' := |A - B|
         A := AsB;
         
      end loop;
      
      -- If we have the degenerate case (A and/or B is 0), then GCD := 0;
      -- Otherwise, GCD := A.
      FZ_Mux(X => A, Y => GCD, Result => GCD, Sel => Degen);
      
      -- Reintroduce the shared power-of-2
      FZ_Quiet_ShiftLeft(N => GCD, ShiftedN => GCD, Count => Indices(Twos));
      
      -- Output result
      Result := GCD;
      
   end FZ_Greatest_Common_Divisor;