Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/random_product_start_systems.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 2: with Standard_Random_Numbers; use Standard_Random_Numbers;
! 3: with Standard_Natural_Vectors;
! 4: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
! 5: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 6: with Random_Product_System;
! 7: with Set_Structure;
! 8:
! 9: package body Random_Product_Start_Systems is
! 10:
! 11: type Boolean_Array is array(integer range <>) of boolean;
! 12:
! 13: function leq ( d1,d2 : Degrees ) return boolean is
! 14:
! 15: -- DESCRIPTION :
! 16: -- returns true if for all k: d1(k) <= d2(k)
! 17:
! 18: begin
! 19: for k in d1'range loop
! 20: if d1(k) > d2(k)
! 21: then return false;
! 22: end if;
! 23: end loop;
! 24: return true;
! 25: end leq;
! 26:
! 27: function Dominated ( d : Degrees; p : Poly ) return boolean is
! 28:
! 29: -- DESCRIPTION :
! 30: -- Returns true if the term indicated by the degrees d
! 31: -- is already in the set structure;
! 32: -- the polynomial p contains all terms that are already
! 33: -- belonging to the set structure.
! 34:
! 35: res : boolean := false;
! 36:
! 37: procedure Scan ( t : in Term; continue : out boolean ) is
! 38: begin
! 39: if leq(d,t.dg)
! 40: then res := true;
! 41: continue := false;
! 42: else continue := true;
! 43: end if;
! 44: end Scan;
! 45: procedure Scan_Terms is new Visiting_Iterator (Scan);
! 46:
! 47: begin
! 48: Scan_Terms(p);
! 49: return res;
! 50: end Dominated;
! 51:
! 52: procedure Divide ( dg : in out Standard_Natural_Vectors.Vector;
! 53: i,d : in natural; occupied : in out Boolean_Array ) is
! 54:
! 55: -- DESCRIPTION :
! 56: -- Divides the monomial by every unknown that occurs already
! 57: -- in one set of the ith component of the set structure.
! 58: -- Every set that contains already one unknown of the monomial is
! 59: -- marked by setting the corresponding entry in occupied on true.
! 60:
! 61: j : natural;
! 62:
! 63: begin
! 64: for k in dg'range loop
! 65: j := 1;
! 66: while (dg(k) /= 0) and (j <= d) loop
! 67: if not occupied(j) and Set_Structure.Is_In(i,j,k)
! 68: then occupied(j) := true;
! 69: dg(k) := dg(k)-1;
! 70: end if;
! 71: j := j+1;
! 72: end loop;
! 73: end loop;
! 74: end Divide;
! 75:
! 76: procedure Augment ( ba : in Boolean_Array; index : in out integer ) is
! 77:
! 78: -- DESCRIPTION :
! 79: -- Augments the index till ba(index) = false.
! 80:
! 81: begin
! 82: while ba(index) loop
! 83: index := index + 1;
! 84: exit when index >= ba'last;
! 85: end loop;
! 86: end Augment;
! 87:
! 88: procedure Build_Set_Structure ( i,di : in natural; p : in Poly ) is
! 89:
! 90: -- DESCRIPTION :
! 91: -- This procedure builds the set structure for the polynomial p
! 92:
! 93: -- ON ENTRY :
! 94: -- i the number of the equation
! 95: -- di the degree of the polynomial
! 96: -- p the polynomial for the i-th equation
! 97:
! 98: temp : Poly := Null_Poly;
! 99:
! 100: procedure Process ( t : in Term; continue : out boolean ) is
! 101:
! 102: ind : natural := 1; -- number of the set
! 103: processed : Boolean_Array(1..di) := (1..di => false);
! 104: tdg : Standard_Natural_Vectors.Vector(t.dg'range) := t.dg.all;
! 105:
! 106: begin
! 107: if not Dominated(t.dg,temp)
! 108: then Divide(tdg,i,di,processed);
! 109: Augment(processed,ind);
! 110: for k in tdg'range loop
! 111: for j in 1..tdg(k) loop
! 112: if not Set_Structure.Is_In(i,ind,k)
! 113: then Set_Structure.Add(i,ind,k);
! 114: processed(ind) := true;
! 115: end if;
! 116: Augment(processed,ind);
! 117: end loop;
! 118: end loop;
! 119: Add(temp,t);
! 120: end if;
! 121: continue := true;
! 122: end Process;
! 123: procedure Process_Terms is new Visiting_Iterator (Process);
! 124:
! 125: begin
! 126: Process_Terms(p);
! 127: Clear(temp);
! 128: end Build_Set_Structure;
! 129:
! 130: procedure Build_Set_Structure ( p : in Poly_Sys ) is
! 131:
! 132: n : natural := p'length;
! 133: ns : Standard_Natural_Vectors.Vector(1..n);
! 134:
! 135: begin
! 136: for i in p'range loop
! 137: ns(i) := Degree(p(i));
! 138: end loop;
! 139: Set_Structure.Init(ns);
! 140: for i in p'range loop
! 141: Build_Set_Structure(i,Degree(p(i)),p(i));
! 142: end loop;
! 143: end Build_Set_Structure;
! 144:
! 145: procedure Build_Random_Product_System ( n : in natural ) is
! 146:
! 147: -- DESCRIPTION :
! 148: -- Based upon the set structure of p, a random product system
! 149: -- will be constructed.
! 150:
! 151: h : Vector(0..n);
! 152: first : boolean;
! 153:
! 154: begin
! 155: for i in 1..n loop
! 156: for j in 1..Set_Structure.Number_Of_Sets(i) loop
! 157: h := (0..n => Create(0.0));
! 158: first := true;
! 159: for k in 1..n loop
! 160: if Set_Structure.Is_In(i,j,k)
! 161: then if first
! 162: then h(k) := Create(1.0); first := false;
! 163: else h(k) := Random1;
! 164: end if;
! 165: end if;
! 166: end loop;
! 167: h(0) := Random1;
! 168: Random_Product_System.Add_Hyperplane(i,h);
! 169: end loop;
! 170: end loop;
! 171: end Build_Random_Product_System;
! 172:
! 173: procedure Construct ( p : in Poly_Sys; q : in out Poly_Sys;
! 174: sols : in out Solution_List ) is
! 175:
! 176: nl : natural := 0;
! 177: n : natural := p'length;
! 178:
! 179: begin
! 180: Build_Set_Structure(p);
! 181: Random_Product_System.Init(n);
! 182: Build_Random_Product_System(n);
! 183: q := Random_Product_System.Polynomial_System;
! 184: Random_Product_System.Solve(sols,nl);
! 185: Random_Product_System.Clear;
! 186: end Construct;
! 187:
! 188: end Random_Product_Start_Systems;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>