[BACK]Return to fewnomial_system_solvers.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift

Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Implift/fewnomial_system_solvers.adb, Revision 1.1

1.1     ! maekawa     1: with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
        !             2: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
        !             3: with Standard_Natural_Vectors;
        !             4: with Standard_Integer_Vectors;
        !             5: with Standard_Integer_VecVecs;           use Standard_Integer_VecVecs;
        !             6: with Standard_Complex_Vectors;           use Standard_Complex_Vectors;
        !             7: with Standard_Complex_Matrices;          use Standard_Complex_Matrices;
        !             8: with Standard_Complex_Linear_Solvers;    use Standard_Complex_Linear_Solvers;
        !             9: with Standard_Complex_Laur_Polys;        use Standard_Complex_Laur_Polys;
        !            10: with Binomial_System_Solvers;            use Binomial_System_Solvers;
        !            11:
        !            12: package body Fewnomial_System_Solvers is
        !            13:
        !            14:   function Equal ( v : Standard_Integer_Vectors.Link_to_Vector; d : Degrees )
        !            15:                  return boolean is
        !            16:   begin
        !            17:     return Standard_Integer_Vectors.Equal
        !            18:              (v,Standard_Integer_Vectors.Link_to_Vector(d));
        !            19:   end Equal;
        !            20:
        !            21:   function Is_Fewnomial_System ( p : Laur_Sys ) return boolean is
        !            22:
        !            23:     da : VecVec(p'range);
        !            24:     nb : natural;
        !            25:     n : natural := p'last - p'first + 1;
        !            26:     cnst : Standard_Integer_Vectors.Link_to_Vector;
        !            27:     res : boolean;
        !            28:
        !            29:     procedure Scan_Term ( t : in Term; cont : out boolean ) is
        !            30:
        !            31:       found : boolean;
        !            32:
        !            33:     begin
        !            34:       found := false;
        !            35:       for i in da'range loop
        !            36:        exit when i > nb;
        !            37:        if Equal(da(i),t.dg)
        !            38:         then found := true;
        !            39:              exit;
        !            40:         end if;
        !            41:       end loop;
        !            42:       if not found
        !            43:        then if nb < n
        !            44:             then nb := nb + 1;
        !            45:                  found := true;
        !            46:                  da(nb) := new Standard_Integer_Vectors.Vector(t.dg'range);
        !            47:                  for j in t.dg'range loop
        !            48:                    da(nb)(j) := t.dg(j);
        !            49:                   end loop;
        !            50:              elsif nb < n + 1
        !            51:                 then nb := nb + 1;
        !            52:                      cnst := new Standard_Integer_Vectors.Vector(t.dg'range);
        !            53:                      found := true;
        !            54:                      for j in t.dg'range loop
        !            55:                        cnst(j) := t.dg(j);
        !            56:                       end loop;
        !            57:                 elsif Equal(cnst,t.dg)
        !            58:                     then found := true;
        !            59:                     else res := false;
        !            60:             end if;
        !            61:       end if;
        !            62:       cont := found;
        !            63:     end Scan_Term;
        !            64:     procedure Scan_Terms is new Visiting_Iterator (Scan_Term);
        !            65:
        !            66:   begin
        !            67:     res := true;
        !            68:     nb := 0;
        !            69:     for i in p'range loop
        !            70:       Scan_Terms(p(i));
        !            71:       exit when not res;
        !            72:     end loop;
        !            73:     Clear(da);
        !            74:     Standard_Integer_Vectors.Clear(cnst);
        !            75:     return res;
        !            76:   end Is_Fewnomial_System;
        !            77:
        !            78:   procedure Scale ( cnst : in Standard_Integer_Vectors.Link_to_Vector;
        !            79:                     da : in out VecVec ) is
        !            80:
        !            81:   -- DESCRIPTION :
        !            82:   --   If cnst is not equal to the zero degrees,
        !            83:   --   then all entries in da will be shifted along cnst.
        !            84:
        !            85:     zero : boolean;
        !            86:
        !            87:   begin
        !            88:     zero := true;
        !            89:     for i in cnst'range loop
        !            90:       if cnst(i) /= 0
        !            91:        then zero := false; exit;
        !            92:       end if;
        !            93:     end loop;
        !            94:     if not zero
        !            95:      then for i in da'range loop
        !            96:            Standard_Integer_Vectors.Sub(da(i),cnst);
        !            97:           end loop;
        !            98:     end if;
        !            99:   end Scale;
        !           100:
        !           101:   procedure Initialize ( p : in Laur_Sys; n : in natural;
        !           102:                         a : in out matrix; b : in out vector;
        !           103:                         da : in out VecVec; fail : out boolean ) is
        !           104:
        !           105:   -- DESCRIPTION :
        !           106:   --   initializes the data needed for the solver;
        !           107:   --   fail becomes true if the system p has too many different terms.
        !           108:
        !           109:     cnst : Standard_Integer_Vectors.Link_to_Vector;
        !           110:     da_numb,cnt : natural;
        !           111:     fl : boolean;
        !           112:
        !           113:     use Standard_Integer_Vectors;
        !           114:
        !           115:     procedure Search_Term ( t : in Term; cont : out boolean ) is
        !           116:
        !           117:       found : boolean;
        !           118:
        !           119:     begin
        !           120:       found := false;
        !           121:       for i in da'range loop
        !           122:         exit when i > da_numb;
        !           123:         if Equal(da(i),t.dg)
        !           124:          then found := true;
        !           125:               a(cnt,i) := t.cf;
        !           126:               exit;
        !           127:         end if;
        !           128:       end loop;
        !           129:       if not found
        !           130:        then
        !           131:          if da_numb < n
        !           132:           then da_numb := da_numb + 1;
        !           133:                da(da_numb) := new Standard_Integer_Vectors.Vector(t.dg'range);
        !           134:                for k in da(da_numb)'range loop
        !           135:                  da(da_numb)(k) := t.dg(k);
        !           136:                end loop;
        !           137:                found := true;
        !           138:                a(cnt,da_numb) := t.cf;
        !           139:           elsif da_numb < n+1
        !           140:               then da_numb := da_numb + 1;
        !           141:                   cnst := new Standard_Integer_Vectors.Vector(t.dg'range);
        !           142:                   for k in cnst'range loop
        !           143:                      cnst(k) := t.dg(k);
        !           144:                    end loop;
        !           145:                    found := true;
        !           146:                    b(cnt) := -t.cf;
        !           147:               elsif Equal(cnst,t.dg)
        !           148:                   then found := true;
        !           149:                        b(cnt) := -t.cf;
        !           150:          end if;
        !           151:       end if;
        !           152:       if not found
        !           153:        then cont := false;
        !           154:            fl := true;
        !           155:        else cont := true;
        !           156:            fl := false;
        !           157:       end if;
        !           158:     end Search_Term;
        !           159:     procedure Search is new Visiting_Iterator(Search_Term);
        !           160:
        !           161:   begin
        !           162:     da_numb := 0;
        !           163:     for i in p'range loop
        !           164:       for j in a'range(2) loop
        !           165:        a(i,j) := Create(0.0);
        !           166:       end loop;
        !           167:       b(i) := Create(0.0);
        !           168:       cnt := i;
        !           169:       Search(p(i));
        !           170:       exit when fl;
        !           171:     end loop;
        !           172:     if not fl and then (cnst /= null)
        !           173:      then Scale(cnst,da);
        !           174:          Standard_Integer_Vectors.Clear(cnst);
        !           175:          fail := false;
        !           176:      else fail := true;
        !           177:     end if;
        !           178:   end Initialize;
        !           179:
        !           180:   procedure Solve ( p : in Laur_Sys; sols : in out Solution_List;
        !           181:                    fail : out boolean ) is
        !           182:
        !           183:     n : natural := p'length;
        !           184:     da : VecVec(p'range);
        !           185:     a : Matrix(1..n,1..n);
        !           186:     b : Vector(1..n);
        !           187:     fl : boolean;
        !           188:
        !           189:   begin
        !           190:     Initialize(p,n,a,b,da,fl);
        !           191:     if fl
        !           192:      then fail := true;
        !           193:      else fail := false;
        !           194:          declare
        !           195:            piv : Standard_Natural_Vectors.Vector(1..n);
        !           196:            info : integer;
        !           197:           begin
        !           198:            lufac(a,n,piv,info);
        !           199:             if info = 0
        !           200:             then lusolve(a,n,piv,b);
        !           201:                  fl := false;
        !           202:                  for i in b'range loop
        !           203:                    if AbsVal(b(i)) + 1.0 = 1.0
        !           204:                     then fl := true;
        !           205:                          exit;
        !           206:                     end if;
        !           207:                   end loop;
        !           208:                  if not fl
        !           209:                    then Solve(da,b,n,sols);
        !           210:                   end if;
        !           211:             end if;
        !           212:           end;
        !           213:     end if;
        !           214:     Clear(da);
        !           215:   end Solve;
        !           216:
        !           217: end Fewnomial_System_Solvers;

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