Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Product/random_product_system.adb, Revision 1.1
1.1 ! maekawa 1: with unchecked_deallocation;
! 2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 3: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 4: with Standard_Natural_Vectors;
! 5: with Standard_Integer_Vectors;
! 6: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
! 7: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
! 8: with Generic_Lists;
! 9:
! 10: package body Random_Product_System is
! 11:
! 12: -- DATA STRUCTURES :
! 13:
! 14: package List_of_Vectors is new Generic_Lists(Link_to_Vector);
! 15: type Equation_List is new List_of_Vectors.List;
! 16:
! 17: type Equation is record
! 18: first,last : Equation_List;
! 19: end record;
! 20:
! 21: type Equations is array(positive range <>) of Equation;
! 22: type Link_To_Equations is access Equations;
! 23: procedure free is new unchecked_deallocation (Equations,Link_To_Equations);
! 24:
! 25: -- INTERNAL DATA :
! 26:
! 27: rps : Link_To_Equations;
! 28:
! 29: Bad_Condition : constant double_float := 10.0**(-12);
! 30:
! 31: -- OPERATIONS :
! 32:
! 33: procedure Init ( n : in natural ) is
! 34: begin
! 35: rps := new Equations(1..n);
! 36: end Init;
! 37:
! 38: procedure Clear ( eql : in out Equation_List ) is
! 39:
! 40: temp : Equation_List := eql;
! 41: lv : Link_To_Vector;
! 42:
! 43: begin
! 44: while not Is_Null(temp) loop
! 45: lv := Head_Of(temp);
! 46: Clear(lv);
! 47: temp := Tail_of(temp);
! 48: end loop;
! 49: List_Of_Vectors.Clear(List_Of_Vectors.List(eql));
! 50: end Clear;
! 51:
! 52: procedure Clear ( eq : in out Equation ) is
! 53: begin
! 54: Clear(eq.first);
! 55: -- eq.last is just a pointer to the last element of eq.first;
! 56: -- if eq.first disappears, then also eq.last does
! 57: end Clear;
! 58:
! 59: procedure Clear ( eqs : in out Equations ) is
! 60: begin
! 61: for i in eqs'range loop
! 62: Clear(eqs(i));
! 63: end loop;
! 64: end Clear;
! 65:
! 66: procedure Clear is
! 67: begin
! 68: if rps /= null
! 69: then for i in rps'range loop
! 70: Clear(rps(i));
! 71: end loop;
! 72: free(rps);
! 73: end if;
! 74: end Clear;
! 75:
! 76: procedure Add_Hyperplane ( i : in natural; h : in Vector ) is
! 77:
! 78: eqi : Equation renames rps(i);
! 79: lh : Link_To_Vector := new Vector'(h);
! 80:
! 81: begin
! 82: if Is_Null(eqi.first)
! 83: then Construct(lh,eqi.first);
! 84: eqi.last := eqi.first;
! 85: else declare
! 86: temp : Equation_List;
! 87: begin
! 88: Construct(lh,temp);
! 89: Swap_Tail(eqi.last,temp);
! 90: eqi.last := Tail_Of(eqi.last);
! 91: end;
! 92: end if;
! 93: end Add_Hyperplane;
! 94:
! 95: function Dimension return natural is
! 96: begin
! 97: if rps = null
! 98: then return 0;
! 99: else return rps'last;
! 100: end if;
! 101: end Dimension;
! 102:
! 103: function Number_of_Hyperplanes ( i : natural ) return natural is
! 104: begin
! 105: if rps = null
! 106: then return 0;
! 107: else return Length_Of(rps(i).first);
! 108: end if;
! 109: end Number_Of_Hyperplanes;
! 110:
! 111: function Get_Hyperplane ( i,j : in natural ) return Link_to_Vector is
! 112:
! 113: nulvec : Link_to_Vector := null;
! 114:
! 115: begin
! 116: if rps = null
! 117: then return nulvec;
! 118: else declare
! 119: eqi : Equation_List := rps(i).first;
! 120: count : natural := 1;
! 121: begin
! 122: while not Is_Null(eqi) loop
! 123: if count = j
! 124: then return Head_Of(eqi);
! 125: else count := count + 1;
! 126: eqi := Tail_Of(eqi);
! 127: end if;
! 128: end loop;
! 129: end;
! 130: return nulvec;
! 131: end if;
! 132: end Get_Hyperplane;
! 133:
! 134: function Get_Hyperplane ( i,j : in natural ) return Vector is
! 135:
! 136: lres : Link_to_Vector := Get_Hyperplane(i,j);
! 137: nulvec : Vector(0..0) := (0..0 => Create(0.0));
! 138:
! 139: begin
! 140: if lres = null
! 141: then return nulvec;
! 142: else return lres.all;
! 143: end if;
! 144: end Get_Hyperplane;
! 145:
! 146: procedure Change_Hyperplane ( i,j : in natural; h : in Vector ) is
! 147: begin
! 148: if rps = null
! 149: then return;
! 150: else declare
! 151: eqi : Equation_List := rps(i).first;
! 152: lv : Link_To_Vector;
! 153: count : natural := 1;
! 154: begin
! 155: while not Is_Null(eqi) loop
! 156: if count = j
! 157: then lv := Head_Of(eqi);
! 158: for k in h'range loop
! 159: lv(k) := h(k);
! 160: end loop;
! 161: return;
! 162: else count := count + 1;
! 163: eqi := Tail_Of(eqi);
! 164: end if;
! 165: end loop;
! 166: end;
! 167: end if;
! 168: end Change_Hyperplane;
! 169:
! 170: procedure Solve ( i,n : in natural; sols,sols_last : in out Solution_List;
! 171: a : in out Matrix; b : in out Vector;
! 172: nb : in out natural ) is
! 173: begin
! 174: if i > n
! 175: then declare
! 176: aa : Matrix(a'range(1),a'range(2));
! 177: bb : Vector(b'range);
! 178: rcond : double_float;
! 179: ipvt : Standard_Natural_Vectors.Vector(b'range);
! 180: begin
! 181: for k in a'range(1) loop
! 182: for l in a'range(2) loop
! 183: aa(k,l) := a(k,l);
! 184: end loop;
! 185: bb(k) := b(k);
! 186: end loop;
! 187: lufco(aa,n,ipvt,rcond);
! 188: nb := nb + 1;
! 189: if abs(rcond) > Bad_Condition
! 190: then lusolve(aa,n,ipvt,bb);
! 191: declare
! 192: s : Solution(n);
! 193: begin
! 194: s.t := Create(0.0);
! 195: s.m := 1;
! 196: s.v := bb;
! 197: s.err := 0.0; s.rco := rcond; s.res := 0.0;
! 198: Append(sols,sols_last,s);
! 199: end;
! 200: end if;
! 201: exception
! 202: when numeric_error => return;
! 203: end;
! 204: else declare
! 205: eqi : Equation_List := rps(i).first;
! 206: h : Vector(0..n);
! 207: count : natural := 0;
! 208: begin
! 209: while not Is_Null(eqi) loop
! 210: count := count + 1;
! 211: h := Head_Of(eqi).all;
! 212: b(i) := -h(0);
! 213: for j in 1..n loop
! 214: a(i,j) := h(j);
! 215: end loop;
! 216: Solve(i+1,n,sols,sols_last,a,b,nb);
! 217: eqi := Tail_Of(eqi);
! 218: end loop;
! 219: end;
! 220: end if;
! 221: end Solve;
! 222:
! 223: procedure Solve ( sols : in out Solution_List; nl : out natural ) is
! 224:
! 225: n : natural := rps'last;
! 226: m : Matrix(1..n,1..n);
! 227: v : Vector(1..n);
! 228: num : natural := 0;
! 229: last : Solution_List;
! 230:
! 231: begin
! 232: for i in 1..n loop
! 233: for j in 1..n loop
! 234: m(i,j) := Create(0.0);
! 235: end loop;
! 236: v(i) := Create(0.0);
! 237: end loop;
! 238: Solve(1,n,sols,last,m,v,num);
! 239: nl := num;
! 240: end Solve;
! 241:
! 242: procedure Solve ( sols : in out Solution_List; nl : out natural;
! 243: l : in List ) is
! 244:
! 245: n : natural := rps'last;
! 246: m : Matrix(1..n,1..n);
! 247: v : Vector(1..n);
! 248: num : natural := 0;
! 249: temp : List := l;
! 250: pos : Standard_Integer_Vectors.Vector(1..n);
! 251: stop : boolean := false;
! 252: last : Solution_List;
! 253:
! 254: procedure PSolve ( i,n : in natural; sols,sols_last : in out Solution_List;
! 255: a : in out Matrix; b : in out Vector;
! 256: nb : in out natural ) is
! 257: begin
! 258: if i > n
! 259: then declare
! 260: aa : Matrix(a'range(1),a'range(2));
! 261: bb : Vector(b'range);
! 262: rcond : double_float;
! 263: ipvt : Standard_Natural_Vectors.Vector(b'range);
! 264: begin
! 265: for k in a'range(1) loop
! 266: for l in a'range(2) loop
! 267: aa(k,l) := a(k,l);
! 268: end loop;
! 269: bb(k) := b(k);
! 270: end loop;
! 271: lufco(aa,n,ipvt,rcond);
! 272: nb := nb + 1;
! 273: if abs(rcond) > Bad_Condition
! 274: then lusolve(aa,n,ipvt,bb);
! 275: declare
! 276: s : Solution(n);
! 277: begin
! 278: s.t := Create(0.0);
! 279: s.m := 1;
! 280: s.v := bb;
! 281: s.err := 0.0; s.rco := rcond; s.res := 0.0;
! 282: Append(sols,sols_last,s);
! 283: end;
! 284: end if;
! 285: if Is_Null(temp)
! 286: then stop := true;
! 287: else pos := Head_Of(temp).all;
! 288: temp := Tail_Of(temp);
! 289: end if;
! 290: exception
! 291: when numeric_error => return;
! 292: end;
! 293: else declare
! 294: eqi : Equation_List := rps(i).first;
! 295: h : Vector(0..n);
! 296: count : natural := 0;
! 297: begin
! 298: while not Is_Null(eqi) loop
! 299: count := count + 1;
! 300: if count = pos(i)
! 301: then h := Head_Of(eqi).all;
! 302: b(i) := -h(0);
! 303: for j in 1..n loop
! 304: a(i,j) := h(j);
! 305: end loop;
! 306: --put("eq"); put(i,1); put(count,1); put(" ");
! 307: PSolve(i+1,n,sols,sols_last,a,b,nb);
! 308: end if;
! 309: exit when stop;
! 310: eqi := Tail_Of(eqi);
! 311: end loop;
! 312: end;
! 313: end if;
! 314: end PSolve;
! 315:
! 316: begin
! 317: if not Is_Null(temp)
! 318: then pos := Head_Of(temp).all;
! 319: temp := Tail_Of(temp);
! 320: for i in 1..n loop
! 321: for j in 1..n loop
! 322: m(i,j) := Create(0.0);
! 323: end loop;
! 324: v(i) := Create(0.0);
! 325: end loop;
! 326: PSolve(1,n,sols,last,m,v,num);
! 327: nl := num;
! 328: end if;
! 329: end Solve;
! 330:
! 331: function Polynomial ( h : Vector ) return Poly is
! 332:
! 333: res : Poly;
! 334: t : Term;
! 335: n : natural := h'last;
! 336:
! 337: begin
! 338: for j in 0..n loop
! 339: if h(j) /= Create(0.0)
! 340: then t.dg := new Standard_Natural_Vectors.Vector'(1..n => 0);
! 341: t.cf := h(j);
! 342: if j /= 0
! 343: then t.dg(j) := 1;
! 344: end if;
! 345: Add(res,t);
! 346: Clear(t.dg);
! 347: end if;
! 348: end loop;
! 349: return res;
! 350: end Polynomial;
! 351:
! 352: function Create ( i : in natural ) return Poly is
! 353:
! 354: eql : Equation_List := rps(i).first;
! 355: hyp,res : Poly := Null_Poly;
! 356:
! 357: begin
! 358: while not Is_Null(eql) loop
! 359: hyp := Polynomial(Head_Of(eql).all);
! 360: if res = Null_Poly
! 361: then Copy(hyp,res);
! 362: else Mul(res,hyp);
! 363: end if;
! 364: Clear(hyp);
! 365: eql := Tail_Of(eql);
! 366: end loop;
! 367: return res;
! 368: end Create;
! 369:
! 370: function Polynomial_System return Poly_Sys is
! 371:
! 372: res : Poly_Sys(rps'range);
! 373:
! 374: begin
! 375: for i in rps'range loop
! 376: res(i) := Create(i);
! 377: end loop;
! 378: return res;
! 379: end Polynomial_System;
! 380:
! 381: end Random_Product_System;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to random_product_system.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