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>