[BACK]Return to greatest_common_divisors.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices

Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/greatest_common_divisors.adb, Revision 1.1

1.1     ! maekawa     1: package body Greatest_Common_Divisors is
        !             2:
        !             3:   use Euclidean_Domain;
        !             4:
        !             5:   function pos_gcd ( a,b : number ) return number is
        !             6:
        !             7:     res,r : number;
        !             8:
        !             9:   begin
        !            10:     if Equal(b,zero)
        !            11:      then Copy(a,res);
        !            12:      else r := Rmd(a,b);
        !            13:           if Equal(r,zero)
        !            14:            then Copy(b,res);
        !            15:            else res := pos_gcd(b,r);
        !            16:           end if;
        !            17:           Clear(r);
        !            18:     end if;
        !            19:     return res;
        !            20:   end pos_gcd;
        !            21:
        !            22:   function gcd ( a,b : number ) return number is
        !            23:
        !            24:     res,mina,minb : number;
        !            25:
        !            26:   begin
        !            27:     if a < zero
        !            28:      then if b < zero
        !            29:           then mina := -a; minb := -b;
        !            30:                 res := pos_gcd(mina,minb);
        !            31:                 Clear(mina); Clear(minb);
        !            32:           else mina := -a;
        !            33:                 res := pos_gcd(mina,b);
        !            34:                 Clear(mina);
        !            35:           end if;
        !            36:      elsif b < zero
        !            37:         then minb := -b;
        !            38:               res := pos_gcd(a,minb);
        !            39:               Clear(minb);
        !            40:         else res := pos_gcd(a,b);
        !            41:     end if;
        !            42:     return res;
        !            43:   end gcd;
        !            44:
        !            45:   function lcm ( a,b : number ) return number is
        !            46:
        !            47:     res : number;
        !            48:     g : number := gcd(a,b);
        !            49:
        !            50:   begin
        !            51:     if Equal(g,zero)
        !            52:      then res := zero;
        !            53:      else res := a/g;
        !            54:           Mul(res,b);
        !            55:     end if;
        !            56:     Clear(g);
        !            57:     return res;
        !            58:   end lcm;
        !            59:
        !            60:   procedure gcd ( a,b : in number; k,l,d : out number ) is
        !            61:
        !            62:     mina,minb,kk,ll : number;
        !            63:
        !            64:     procedure pos_gcd ( a,b : in number; k,l,d : out number ) is
        !            65:
        !            66:       r,aa,bb,k1,l1,k0,l0,h,tmp : number;
        !            67:      -- nk1,nl1 : number;
        !            68:
        !            69:     -- REQUIRED :
        !            70:     --   a >= b and a > 0 and b > 0
        !            71:
        !            72:     begin
        !            73:       r := Rmd(a,b);
        !            74:       if Equal(r,zero)
        !            75:        then k := zero; l := one; Copy(b,d);
        !            76:        else k0 := zero; l0 := one;
        !            77:             Copy(r,bb);
        !            78:             Copy(b,aa);
        !            79:             k1 := one; l1 := - a / b;
        !            80:             loop
        !            81:               r := Rmd(aa,bb);
        !            82:               exit when Equal(r,zero);
        !            83:               h := aa / bb;
        !            84:             -- nk1 := k0 - h*k1; Clear(k0); k0 := k1; k1 := nk1;
        !            85:              -- nl1 := l0 - h*l1; Clear(l0); l0 := l1; l1 := nl1;
        !            86:               tmp := k1; k1 := k0 - h*k1; k0 := tmp;
        !            87:               tmp := l1; l1 := l0 - h*l1; l0 := tmp;
        !            88:               Clear(h);
        !            89:               Copy(bb,aa);
        !            90:              Copy(r,bb);
        !            91:               Clear(r);
        !            92:             end loop;
        !            93:            k := k1; l := l1; d := bb;
        !            94:       end if;
        !            95:     end pos_gcd;
        !            96:
        !            97:   begin
        !            98:     if Equal(a,zero)
        !            99:      then if b < zero
        !           100:            then d := -b; k := zero; l := -one;
        !           101:            else Copy(b,d); k := zero; l := one;
        !           102:           end if;
        !           103:          return;
        !           104:      elsif Equal(b,zero)
        !           105:          then if a < zero
        !           106:                then d := -a; k := -one; l := zero;
        !           107:                else Copy(a,d); k := one; l := zero;
        !           108:               end if;
        !           109:              return;
        !           110:     end if;
        !           111:     if a < b
        !           112:      then if b < zero
        !           113:            then minb := -b; mina := -a;
        !           114:                 pos_gcd(minb,mina,ll,kk,d);
        !           115:                 Clear(minb); Clear(mina);
        !           116:                 Min(kk); k := kk;
        !           117:                 Min(ll); l := ll;
        !           118:           elsif a < zero
        !           119:               then mina := -a;
        !           120:                     pos_gcd(b,mina,l,kk,d);
        !           121:                     Clear(mina);
        !           122:                    Min(kk); k := kk;
        !           123:               else pos_gcd(b,a,l,k,d);
        !           124:           end if;
        !           125:      else -- a >= b
        !           126:           if a < zero
        !           127:            then mina := -a; minb := -b;
        !           128:                 pos_gcd(mina,minb,kk,ll,d);
        !           129:                 Clear(mina); Clear(minb);
        !           130:                 Min(kk); k := kk;
        !           131:                 Min(ll); l := ll;
        !           132:            elsif b < zero
        !           133:                then minb := -b;
        !           134:                     pos_gcd(a,minb,k,ll,d);
        !           135:                     Clear(minb);
        !           136:                    Min(ll); l := ll;
        !           137:               else pos_gcd(a,b,k,l,d);
        !           138:           end if;
        !           139:     end if;
        !           140:   end gcd;
        !           141:
        !           142: end Greatest_Common_Divisors;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>