Return to greatest_common_divisors.adb CVS log

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

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;