Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/degree_structure.adb, Revision 1.1
1.1 ! maekawa 1: with unchecked_deallocation;
! 2: with text_io; use text_io;
! 3: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 4: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 5: with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
! 6: with Standard_Random_Numbers; use Standard_Random_Numbers;
! 7: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
! 8: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
! 9: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 10: with Degrees_in_Sets_of_Unknowns; use Degrees_in_Sets_of_Unknowns;
! 11: with Generate_Unions;
! 12:
! 13: package body Degree_Structure is
! 14:
! 15: -- DECLARATIONS :
! 16:
! 17: type zd(m : natural) is record
! 18: z : Partition(1..m);
! 19: d : Standard_Natural_Vectors.Vector(1..m);
! 20: end record;
! 21: type Link_To_zd is access zd;
! 22: procedure free is new unchecked_deallocation(zd,Link_To_zd);
! 23:
! 24: type dgst is array(positive range <>) of Link_To_zd;
! 25: type Link_To_dgst is access dgst;
! 26: procedure free is new unchecked_deallocation(dgst,Link_To_dgst);
! 27:
! 28: -- INTERNAL DATA :
! 29:
! 30: ds : Link_To_dgst;
! 31:
! 32: -- CREATORS :
! 33:
! 34: procedure Find_Partition ( p : in Poly;
! 35: z : in out Partition; m : in out natural;
! 36: dg : in out Standard_Natural_Vectors.Vector ) is
! 37:
! 38: -- DESCRIPTION :
! 39: -- This routine finds an good partition for the polynomial p
! 40:
! 41: -- ON ENTRY :
! 42: -- p a polynomial.
! 43:
! 44: -- ON RETURN :
! 45: -- z a partition of the set of unknowns of p;
! 46: -- m the number of sets in the partition z;
! 47: -- dg the degrees of the polynomial p in the sets of z,
! 48: -- dg(i) = Degree(p,z(i)).
! 49:
! 50: n : natural := Number_Of_Unknowns(p);
! 51: di : integer;
! 52: added : boolean;
! 53:
! 54: begin
! 55: for i in 1..n loop
! 56: di := Degree(p,i);
! 57: if di > 0
! 58: then added := false;
! 59: for j in 1..m loop
! 60: if di = dg(j)
! 61: then Add(z(j),i);
! 62: if Degree(p,z(j)) = dg(j)
! 63: then added := true;
! 64: else Remove(z(j),i);
! 65: end if;
! 66: end if;
! 67: exit when added;
! 68: end loop;
! 69: if not added
! 70: then m := m + 1;
! 71: Add(z(m),i);
! 72: dg(m) := di;
! 73: end if;
! 74: end if;
! 75: end loop;
! 76: end Find_Partition;
! 77:
! 78: procedure Create ( p : in Poly_Sys ) is
! 79:
! 80: n : natural := p'length;
! 81: z : Partition(1..n);
! 82: d : Standard_Natural_Vectors.Vector(1..n) := (1..n => 0);
! 83: m : natural;
! 84:
! 85: begin
! 86: if ds /= null
! 87: then Clear;
! 88: end if;
! 89: ds := new dgst(1..n);
! 90: for i in p'range loop
! 91: m := 0;
! 92: Create(z,n);
! 93: Find_Partition(p(i),z,m,d);
! 94: ds(i) := new zd(m);
! 95: for j in 1..m loop
! 96: ds(i).z(j) := Create(z(j));
! 97: ds(i).d(j) := d(j);
! 98: end loop;
! 99: Clear(z);
! 100: end loop;
! 101: end Create;
! 102:
! 103: procedure Put ( p : in Poly_Sys;
! 104: i,m : in natural; z : in Partition ) is
! 105:
! 106: n : natural := p'length;
! 107:
! 108: begin
! 109: if ds = null
! 110: then ds := new dgst(1..n);
! 111: end if;
! 112: ds(i) := new zd(m);
! 113: for j in 1..m loop
! 114: ds(i).z(j) := Create(z(j));
! 115: ds(i).d(j) := Degree(p(i),z(j));
! 116: end loop;
! 117: end Put;
! 118:
! 119: -- SELECTORS :
! 120:
! 121: function Empty return boolean is
! 122: begin
! 123: return (ds = null);
! 124: end Empty;
! 125:
! 126: function Get ( i : natural ) return natural is
! 127: begin
! 128: return ds(i).m;
! 129: end Get;
! 130:
! 131: procedure Get ( i : in natural; z : in out Partition;
! 132: d : out Standard_Natural_Vectors.Vector ) is
! 133: begin
! 134: for j in 1..ds(i).m loop
! 135: z(j) := Create(ds(i).z(j));
! 136: d(j) := ds(i).d(j);
! 137: end loop;
! 138: end Get;
! 139:
! 140: -- COMPUTING THE GENERALIZED BEZOUT NUMBER :
! 141:
! 142: function Matrix_Criterion ( z : Partition ) return boolean is
! 143:
! 144: -- DESCRIPTION :
! 145: -- This is the Matrix criterion for testing if
! 146: -- a product is admissible or not.
! 147:
! 148: n : natural := z'last - z'first +1;
! 149: mat : Matrix(1..n,1..n);
! 150: ipvt : Standard_Natural_Vectors.Vector(1..n);
! 151: eps : constant double_float := 10.0**(-10);
! 152: r : Complex_Number;
! 153: rcond : double_float;
! 154:
! 155: begin
! 156: for i in 1..n loop
! 157: r := Create(double_float(i+1));
! 158: for j in 1..n loop
! 159: if Is_In(z(i),j)
! 160: then mat(i,j) := r;
! 161: r := r*r;
! 162: else mat(i,j) := Create(0.0);
! 163: end if;
! 164: end loop;
! 165: end loop;
! 166: lufco(mat,n,ipvt,rcond);
! 167: return (abs(rcond) > eps);
! 168: exception
! 169: when others => return false;
! 170: end Matrix_Criterion;
! 171:
! 172: function Admissible ( z : Partition; n : natural ) return boolean is
! 173:
! 174: temp : Partition(1..n);
! 175: admis : boolean := true;
! 176:
! 177: begin
! 178: temp(1) := Create(z(1));
! 179: for i in 2..(n-1) loop
! 180: temp(i) := Create(z(i));
! 181: admis := Admissible(temp,i,z(i+1));
! 182: exit when not admis;
! 183: end loop;
! 184: Clear(temp);
! 185: return admis;
! 186: end Admissible;
! 187:
! 188: function Admissible ( z : Partition; n : natural; s : Set )
! 189: return boolean is
! 190: begin
! 191: for k in 1..n loop
! 192: if not Admissible(z,k,n,s)
! 193: then return false;
! 194: end if;
! 195: end loop;
! 196: return true;
! 197: end Admissible;
! 198:
! 199: function Admissible ( z : Partition; k,n : natural; s : Set )
! 200: return boolean is
! 201:
! 202: type arr is array (integer range <>) of boolean;
! 203: admis : boolean := true;
! 204:
! 205: procedure check (a : in arr; continue : out boolean) is
! 206:
! 207: u : Set := Create(s);
! 208:
! 209: begin
! 210: for i in a'range loop
! 211: if a(i)
! 212: then Union(u,z(i));
! 213: end if;
! 214: end loop;
! 215: admis := ( Extent_Of(u) >= k+1 );
! 216: continue := admis;
! 217: Clear(u);
! 218: end check;
! 219:
! 220: procedure gen is new Generate_Unions(arr,check);
! 221:
! 222: begin
! 223: gen(k,1,n);
! 224: return admis;
! 225: end Admissible;
! 226:
! 227: procedure Compute ( i,n,sum : in natural; res : in out natural;
! 228: z : in out Partition ) is
! 229: begin
! 230: if i > n
! 231: then res := res + sum;
! 232: else -- Pick out a set and check if it is allowed :
! 233: for j in 1..ds(i).m loop
! 234: if ds(i).d(j) /= 0 and then Admissible(z,i-1,ds(i).z(j))
! 235: then z(i) := Create(ds(i).z(j));
! 236: Compute(i+1,n,sum*ds(i).d(j),res,z);
! 237: Clear(z(i));
! 238: end if;
! 239: end loop;
! 240: end if;
! 241: end Compute;
! 242:
! 243: function Generalized_Bezout_Number return natural is
! 244:
! 245: res : natural := 0;
! 246: n : natural := ds'length;
! 247: z : Partition(1..n);
! 248:
! 249: begin
! 250: Compute(1,n,1,res,z);
! 251: return res;
! 252: end Generalized_Bezout_Number;
! 253:
! 254: function Generalized_Bezout_Number ( p : in Poly_Sys ) return natural is
! 255: begin
! 256: Create(p);
! 257: return Generalized_Bezout_Number;
! 258: end Generalized_Bezout_Number;
! 259:
! 260: -- DESTRUCTOR :
! 261:
! 262: procedure Clear is
! 263: begin
! 264: if ds /= null
! 265: then for i in ds'range loop
! 266: free(ds(i));
! 267: end loop;
! 268: free(ds);
! 269: end if;
! 270: end Clear;
! 271:
! 272: end Degree_Structure;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to degree_structure.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Product