File : fz_gcd.adb

   1 ------------------------------------------------------------------------------
   2 ------------------------------------------------------------------------------
   3 -- This file is part of 'Finite Field Arithmetic', aka 'FFA'.               --
   4 --                                                                          --
   5 -- (C) 2019 Stanislav Datskovskiy ( www.loper-os.org )                      --
   6 -- http://wot.deedbot.org/17215D118B7239507FAFED98B98228A001ABFFC7.html     --
   7 --                                                                          --
   8 -- You do not have, nor can you ever acquire the right to use, copy or      --
   9 -- distribute this software ; Should you use this software for any purpose, --
  10 -- or copy and distribute it to anyone or in any manner, you are breaking   --
  11 -- the laws of whatever soi-disant jurisdiction, and you promise to         --
  12 -- continue doing so for the indefinite future. In any case, please         --
  13 -- always : read and understand any software ; verify any PGP signatures    --
  14 -- that you use - for any purpose.                                          --
  15 --                                                                          --
  16 -- See also http://trilema.com/2015/a-new-software-licensing-paradigm .     --
  17 ------------------------------------------------------------------------------
  18 
  19 with Words;    use Words;
  20 with FZ_Basic; use FZ_Basic;
  21 with FZ_Shift; use FZ_Shift;
  22 with FZ_QShft; use FZ_QShft;
  23 with FZ_Arith; use FZ_Arith;
  24 with FZ_Pred;  use FZ_Pred;
  25 
  26 
  27 package body FZ_GCD is
  28    
  29    -- Find Greatest Common Divisor (GCD) of X and Y.
  30    -- Note that by convention, GCD(0, 0) = 0.
  31    procedure FZ_Greatest_Common_Divisor(X      : in  FZ;
  32                                         Y      : in  FZ;
  33                                         Result : out FZ) is
  34       
  35       -- Widths of X, Y, and Result are equal
  36       subtype Width is Word_Index range X'Range;
  37       
  38       -- Working buffers for GCD computation, initially equal to the inputs
  39       A      : FZ(Width) := X; -- GCD will appear in A in the end
  40       B      : FZ(Width) := Y;
  41       
  42       -- Evenness (negation of lowest bit) of A and B respectively
  43       Ae, Be : WBool;
  44       
  45       -- Common power-of-2 factor
  46       Twos   : Word := 0;
  47       
  48       -- |A - B|
  49       D      : FZ(Width);
  50       
  51       -- This flag is set iff A < B
  52       A_lt_B : WBool;
  53       
  54    begin
  55       
  56       -- For convergence, requires number of shots equal to 2 * FZ_Bitness:
  57       for i in 1 .. 2 * FZ_Bitness(X) loop
  58          
  59          -- If A is even, A := A >> 1; otherwise A := A
  60          Ae := 1 - FZ_OddP(A);
  61          FZ_ShiftRight(A, A, WBit_Index(Ae));
  62          
  63          -- If B is even, B := B >> 1; otherwise B := B
  64          Be := 1 - FZ_OddP(B);
  65          FZ_ShiftRight(B, B, WBit_Index(Be));
  66          
  67          -- If both A and B were even, increment the common power-of-two
  68          Twos := Twos + (Ae and Be);
  69          
  70          -- D := |A - B|
  71          FZ_Sub_Abs(X => A, Y => B, Difference => D, Underflow => A_lt_B);
  72          
  73          -- B' := min(A, B)
  74          FZ_Mux(X => B, Y => A, Result => B, Sel => A_lt_B);
  75          
  76          -- A' := |A - B|
  77          A := D;
  78          
  79       end loop;
  80       
  81       -- Reintroduce the common power-of-2 factor stored in 'Twos'
  82       FZ_Quiet_ShiftLeft(N => A, ShiftedN => A, Count => Indices(Twos));
  83       
  84       -- Output final result
  85       Result := A;
  86       
  87    end FZ_Greatest_Common_Divisor;
  88    
  89 end FZ_GCD;