Return to greatest_common_divisors.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices

Import the second public release of PHCpack.

OKed by Jan Verschelde.

package body Greatest_Common_Divisors is

  use Euclidean_Domain;

  function pos_gcd ( a,b : number ) return number is

    res,r : number;

  begin
    if Equal(b,zero)
     then Copy(a,res);
     else r := Rmd(a,b);
          if Equal(r,zero)
           then Copy(b,res);
           else res := pos_gcd(b,r);
          end if;
          Clear(r);
    end if;
    return res;
  end pos_gcd;

  function gcd ( a,b : number ) return number is

    res,mina,minb : number;

  begin
    if a < zero
     then if b < zero
	   then mina := -a; minb := -b;
                res := pos_gcd(mina,minb);
                Clear(mina); Clear(minb);
	   else mina := -a;
                res := pos_gcd(mina,b);
                Clear(mina);
          end if;
     elsif b < zero
	 then minb := -b;
              res := pos_gcd(a,minb);
              Clear(minb);
	 else res := pos_gcd(a,b);
    end if;
    return res;
  end gcd;

  function lcm ( a,b : number ) return number is

    res : number;
    g : number := gcd(a,b);

  begin
    if Equal(g,zero)
     then res := zero;
     else res := a/g;
          Mul(res,b);
    end if;
    Clear(g);
    return res;
  end lcm;

  procedure gcd ( a,b : in number; k,l,d : out number ) is

    mina,minb,kk,ll : number;

    procedure pos_gcd ( a,b : in number; k,l,d : out number ) is

      r,aa,bb,k1,l1,k0,l0,h,tmp : number;
     -- nk1,nl1 : number;

    -- REQUIRED :
    --   a >= b and a > 0 and b > 0

    begin
      r := Rmd(a,b);
      if Equal(r,zero)
       then k := zero; l := one; Copy(b,d);
       else k0 := zero; l0 := one;
            Copy(r,bb);
            Copy(b,aa);
            k1 := one; l1 := - a / b;
            loop
              r := Rmd(aa,bb);
              exit when Equal(r,zero);
              h := aa / bb;
  	     -- nk1 := k0 - h*k1; Clear(k0); k0 := k1; k1 := nk1;
             -- nl1 := l0 - h*l1; Clear(l0); l0 := l1; l1 := nl1;
              tmp := k1; k1 := k0 - h*k1; k0 := tmp;
              tmp := l1; l1 := l0 - h*l1; l0 := tmp;
              Clear(h);
              Copy(bb,aa);
	      Copy(r,bb);
              Clear(r);
            end loop;
	    k := k1; l := l1; d := bb;
      end if;
    end pos_gcd;

  begin
    if Equal(a,zero)
     then if b < zero
           then d := -b; k := zero; l := -one;
           else Copy(b,d); k := zero; l := one;
          end if;
	  return;
     elsif Equal(b,zero)
         then if a < zero
               then d := -a; k := -one; l := zero;
               else Copy(a,d); k := one; l := zero;
              end if;
	      return;
    end if;
    if a < b
     then if b < zero
           then minb := -b; mina := -a;
                pos_gcd(minb,mina,ll,kk,d);
                Clear(minb); Clear(mina);
                Min(kk); k := kk;
                Min(ll); l := ll;
  	   elsif a < zero
	       then mina := -a;
                    pos_gcd(b,mina,l,kk,d);
                    Clear(mina);
	            Min(kk); k := kk;
	       else pos_gcd(b,a,l,k,d);
          end if;
     else -- a >= b
          if a < zero
           then mina := -a; minb := -b;
                pos_gcd(mina,minb,kk,ll,d);
                Clear(mina); Clear(minb);
                Min(kk); k := kk;
                Min(ll); l := ll;
           elsif b < zero
               then minb := -b;
                    pos_gcd(a,minb,k,ll,d);
                    Clear(minb);
          	    Min(ll); l := ll;
  	       else pos_gcd(a,b,k,l,d);
          end if;
    end if;
  end gcd;

end Greatest_Common_Divisors;