Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Implift/volumes.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Integer_Vectors; use Standard_Integer_Vectors;
! 2: with Standard_Integer_Norms; use Standard_Integer_Norms;
! 3: with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
! 4: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
! 5: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
! 6: with Transformations; use Transformations;
! 7: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
! 8: with Integer_Vectors_Utilities; use Integer_Vectors_Utilities;
! 9: with Integer_Support_Functions; use Integer_Support_Functions;
! 10: with Integer_Face_Enumerators; use Integer_Face_Enumerators;
! 11: with Lists_of_Vectors_Utilities; use Lists_of_Vectors_Utilities;
! 12: with Arrays_of_Lists_Utilities; use Arrays_of_Lists_Utilities;
! 13: --with Face_Enumerator_of_Sum; use Face_Enumerator_of_Sum;
! 14:
! 15: package body Volumes is
! 16:
! 17: -- INVARIANT CONDITION :
! 18: -- In order to get p(v) > 0, the zero vector must be in the first list,
! 19: -- so shifting is necessary.
! 20: -- The shift vector must be equal to the maximal element in the list,
! 21: -- w.r.t. the graded lexicographic ordening.
! 22: -- In this way, the shift vector is unique, which allows to `mirror'
! 23: -- the same operations for the mixed homotopy continuation.
! 24:
! 25: -- AUXILIARIES :
! 26:
! 27: function Create_Facet ( n : natural; facet : Vector; pts : VecVec )
! 28: return VecVec is
! 29:
! 30: res : VecVec(1..n);
! 31: cnt : natural := 0;
! 32:
! 33: begin
! 34: for i in facet'range loop
! 35: cnt := cnt+1;
! 36: res(cnt) := pts(facet(i));
! 37: end loop;
! 38: return res;
! 39: end Create_Facet;
! 40:
! 41: function Determinant ( vv : VecVec ) return integer is
! 42:
! 43: a : Matrix(vv'range,vv'range);
! 44:
! 45: begin
! 46: for k in a'range(1) loop
! 47: for l in a'range(2) loop
! 48: a(k,l) := vv(k)(l);
! 49: end loop;
! 50: end loop;
! 51: -- Upper_Triangulate(a);
! 52: return Det(a);
! 53: end Determinant;
! 54:
! 55: -- VOLUME COMPUTATIONS :
! 56:
! 57: function Vol ( n : natural; l : List ) return natural is
! 58:
! 59: -- DESCRIPTION :
! 60: -- Computes the volume of the simplex spanned by the list
! 61: -- and the origin.
! 62:
! 63: res : integer;
! 64: vv : VecVec(1..n);
! 65: tmp : List;
! 66:
! 67: begin
! 68: tmp := l;
! 69: for i in vv'range loop
! 70: vv(i) := Head_Of(tmp);
! 71: tmp := Tail_Of(tmp);
! 72: end loop;
! 73: res := Determinant(vv);
! 74: if res < 0
! 75: then return -res;
! 76: else return res;
! 77: end if;
! 78: end Vol;
! 79:
! 80: function Volume ( n,i : natural; l : List; v : Vector; pv : integer )
! 81: return natural is
! 82:
! 83: -- DESCRIPTION :
! 84: -- The points in l all belong to the same hyperplane:
! 85: -- < x , v > - pv = 0;
! 86: -- this procedure computes the volume of the polytope generated by
! 87: -- the points in l and the origin.
! 88:
! 89: ll : natural := Length_Of(l);
! 90:
! 91: begin
! 92: if ll < n
! 93: then return 0;
! 94: elsif ll = n
! 95: then return Vol(n,l);
! 96: else declare
! 97: t : Transfo;
! 98: tl,rl : List;
! 99: vl : integer;
! 100: begin
! 101: if pv > 0
! 102: then t := Build_Transfo(v,i);
! 103: tl := t*l;
! 104: rl := Reduce(tl,i);
! 105: Clear(t); Deep_Clear(tl);
! 106: -- Remove_Duplicates(rl);
! 107: if Length_Of(rl) >= n-1
! 108: then vl := pv*Volume(n-1,rl);
! 109: else vl := 0;
! 110: end if;
! 111: Deep_Clear(rl);
! 112: return vl;
! 113: else return 0;
! 114: end if;
! 115: end;
! 116: end if;
! 117: end Volume;
! 118:
! 119: function Sum ( n,m : natural; l : List ) return natural is
! 120:
! 121: -- DESCRIPTION :
! 122: -- Computes the volume of the polytope generated by the points in l;
! 123: -- where n > 1 and n < m = Length_Of(l).
! 124:
! 125: res : natural;
! 126: fix : Link_to_Vector;
! 127: pts : VecVec(1..m-1);
! 128: nulvec : Vector(1..n) := (1..n => 0);
! 129: vl,wl,vl_last : List;
! 130: sh : boolean;
! 131:
! 132: procedure Compute_Facet ( facet : in vector; cont : out boolean ) is
! 133:
! 134: vv : constant VecVec := Create_Facet(n,facet,pts);
! 135: pl : List;
! 136: v : Link_to_Vector;
! 137: i,pv,sup : integer;
! 138:
! 139: begin
! 140: v := Compute_Normal(vv);
! 141: i := Pivot(v);
! 142: if i <= v'last
! 143: then pv := vv(vv'first)*v;
! 144: if pv < 0
! 145: then for j in v'range loop
! 146: v(j) := -v(j);
! 147: end loop;
! 148: pv := -pv;
! 149: end if;
! 150: if (pv > 0) and then not Is_In(vl,v)
! 151: then sup := Maximal_Support(wl,v.all);
! 152: if sup = pv
! 153: then Append(vl,vl_last,v);
! 154: pl := Face(wl,v.all,pv);
! 155: res := res + Volume(n,i,pl,v.all,pv);
! 156: Deep_Clear(pl);
! 157: end if;
! 158: end if;
! 159: end if;
! 160: cont := true;
! 161: end Compute_Facet;
! 162:
! 163: procedure Compute_Facets is new Enumerate_Faces(Compute_Facet);
! 164:
! 165: begin
! 166: res := 0;
! 167: if not Is_In(l,nulvec)
! 168: then fix := Graded_Max(l);
! 169: wl := Shift(l,fix); sh := true;
! 170: else wl := l; sh := false;
! 171: end if;
! 172: Move_to_Front(wl,nulvec);
! 173: pts := Shallow_Create(Tail_Of(wl));
! 174: Compute_Facets(n-1,pts);
! 175: Deep_Clear(vl);
! 176: if sh
! 177: then Deep_Clear(wl); Clear(fix);
! 178: end if;
! 179: return res;
! 180: end Sum;
! 181:
! 182: function Volume ( n : natural; l : List ) return natural is
! 183: m : natural;
! 184: begin
! 185: if n = 1
! 186: then declare
! 187: min,max : integer := 0;
! 188: d : Link_to_Vector := Head_Of(l);
! 189: begin
! 190: Min_Max(l,d'first,min,max);
! 191: return (max - min);
! 192: end;
! 193: else m := Length_Of(l);
! 194: if m <= n
! 195: then return 0;
! 196: else return Sum(n,m,l);
! 197: end if;
! 198: end if;
! 199: end Volume;
! 200:
! 201: function Volume ( n,i : natural; l : List; v : Vector;
! 202: pv : integer; tv : Tree_of_Vectors ) return natural is
! 203:
! 204: -- DESCRIPTION :
! 205: -- The points in l all belong to the same hyperplane:
! 206: -- < x , v > - pv = 0;
! 207: -- this procedure computes the volume of the polytope generated by
! 208: -- the points in l and the origin.
! 209:
! 210: ll : natural := Length_Of(l);
! 211:
! 212: begin
! 213: if ll < n
! 214: then return 0;
! 215: elsif ll = n
! 216: then return Vol(n,l);
! 217: else declare
! 218: t : Transfo;
! 219: tl,rl : List;
! 220: vl : integer;
! 221: begin
! 222: if pv > 0
! 223: then t := Build_Transfo(v,i);
! 224: tl := t*l;
! 225: rl := Reduce(tl,i);
! 226: Clear(t); Deep_Clear(tl);
! 227: -- Remove_Duplicates(rl);
! 228: if Length_Of(rl) >= n-1
! 229: then vl := pv*Volume(n-1,rl,tv);
! 230: else vl := 0;
! 231: end if;
! 232: Deep_Clear(rl);
! 233: return vl;
! 234: else return 0;
! 235: end if;
! 236: end;
! 237: end if;
! 238: end Volume;
! 239:
! 240: function Sum ( n,m : natural; l : List; tv : Tree_of_Vectors )
! 241: return natural is
! 242:
! 243: -- DESCRIPTION :
! 244: -- Computes the volume of the polytope generated by the points in l;
! 245: -- where n > 1 and n < m = Length_Of(l).
! 246: -- The tree of degrees tv is not empty.
! 247:
! 248: res : natural;
! 249: fix : Link_to_Vector;
! 250: nulvec : Vector(1..n) := (1..n => 0);
! 251: wl : List;
! 252: tmp : Tree_of_Vectors;
! 253: sh : boolean;
! 254:
! 255: begin
! 256: res := 0;
! 257: if not Is_In(l,nulvec)
! 258: then fix := Graded_Max(l);
! 259: wl := Shift(l,fix); sh := true;
! 260: else wl := l; sh := false;
! 261: end if;
! 262: Move_to_Front(wl,nulvec);
! 263: tmp := tv;
! 264: while not Is_Null(tmp) loop
! 265: declare
! 266: nd : node := Head_Of(tmp);
! 267: v : Link_to_Vector := nd.d;
! 268: pv,i : integer;
! 269: pl : List;
! 270: begin
! 271: i := Pivot(v);
! 272: pv := Maximal_Support(wl,v.all);
! 273: pl := Face(wl,v.all,pv);
! 274: if (nd.ltv = null) or else Is_Null(nd.ltv.all)
! 275: then res := res + Volume(n,i,pl,v.all,pv);
! 276: else res := res + Volume(n,i,pl,v.all,pv,nd.ltv.all);
! 277: end if;
! 278: Deep_Clear(pl);
! 279: end;
! 280: tmp := Tail_Of(tmp);
! 281: end loop;
! 282: if sh
! 283: then Deep_Clear(wl); Clear(fix);
! 284: end if;
! 285: return res;
! 286: end Sum;
! 287:
! 288: function Volume ( n : natural; l : List; tv : Tree_of_Vectors )
! 289: return natural is
! 290: m : natural;
! 291: begin
! 292: if n = 1
! 293: then declare
! 294: min,max : integer := 0;
! 295: d : Link_to_Vector := Head_Of(l);
! 296: begin
! 297: Min_Max(l,d'first,min,max);
! 298: return (max - min);
! 299: end;
! 300: else m := Length_Of(l);
! 301: if m <= n
! 302: then return 0;
! 303: elsif not Is_Null(tv)
! 304: then return Sum(n,m,l,tv);
! 305: else return Sum(n,m,l);
! 306: end if;
! 307: end if;
! 308: end Volume;
! 309:
! 310: procedure Volume ( n,i : in natural; l : in List; v : in Vector;
! 311: pv : in integer; tv : in out Tree_of_Vectors;
! 312: vlm : out natural ) is
! 313:
! 314: -- DESCRIPTION :
! 315: -- The points in l all belong to the same hyperplane:
! 316: -- < x , v > - pv = 0;
! 317: -- this procedure computes the volume of the polytope generated by
! 318: -- the points in l and the origin.
! 319:
! 320: ll : natural := Length_Of(l);
! 321: vl : natural;
! 322:
! 323: begin
! 324: if ll < n
! 325: then vlm := 0;
! 326: elsif ll = n
! 327: then vl := Vol(n,l);
! 328: if vl > 0
! 329: then declare
! 330: nd : node;
! 331: begin
! 332: nd.d := new Vector'(v);
! 333: Construct(nd,tv);
! 334: end;
! 335: end if;
! 336: vlm := vl;
! 337: else declare
! 338: t : Transfo;
! 339: tl,rl : List;
! 340: begin
! 341: if pv > 0
! 342: then t := Build_Transfo(v,i);
! 343: tl := t*l;
! 344: rl := Reduce(tl,i);
! 345: Clear(t); Deep_Clear(tl);
! 346: -- Remove_Duplicates(rl);
! 347: if Length_Of(rl) >= n-1
! 348: then declare
! 349: tmp : Tree_of_Vectors;
! 350: begin
! 351: Volume(n-1,rl,tmp,vl);
! 352: if vl = 0
! 353: then Clear(tmp);
! 354: else
! 355: declare
! 356: nd : node;
! 357: begin
! 358: nd.d := new Vector'(v);
! 359: if not Is_Null(tmp)
! 360: then nd.ltv := new Tree_of_Vectors'(tmp);
! 361: end if;
! 362: Construct(nd,tv);
! 363: end;
! 364: end if;
! 365: end;
! 366: vlm := pv*vl;
! 367: else vlm := 0;
! 368: end if;
! 369: Deep_Clear(rl);
! 370: else vlm := 0;
! 371: end if;
! 372: end;
! 373: end if;
! 374: end Volume;
! 375:
! 376: procedure Sum ( n,m : in natural; l : in List;
! 377: tv : in out Tree_of_Vectors; vlm : out natural ) is
! 378:
! 379: -- DESCRIPTION :
! 380: -- Computes the volume of the polytope generated by the points in l;
! 381: -- where n > 1 and n < m = Length_Of(l).
! 382:
! 383: res : natural;
! 384: pts : VecVec(1..m-1);
! 385: fix : Link_to_Vector;
! 386: nulvec : Vector(1..n) := (1..n => 0);
! 387: wl : List;
! 388: sh : boolean;
! 389:
! 390: procedure Compute_Facet ( facet : in vector; cont : out boolean ) is
! 391:
! 392: vv : constant VecVec := Create_Facet(n,facet,pts);
! 393: pl : List;
! 394: v : Link_to_Vector;
! 395: i,pv,sup,pvol : integer;
! 396:
! 397: begin
! 398: v := Compute_Normal(vv);
! 399: i := Pivot(v);
! 400: if i <= v'last
! 401: then pv := vv(vv'first)*v;
! 402: if pv < 0
! 403: then for j in v'range loop
! 404: v(j) := -v(j);
! 405: end loop;
! 406: pv := -pv;
! 407: end if;
! 408: if (pv > 0) and then not Is_In(tv,v)
! 409: then sup := Maximal_Support(wl,v.all);
! 410: if sup = pv
! 411: then pl := Face(wl,v.all,pv);
! 412: Volume(n,i,pl,v.all,pv,tv,pvol);
! 413: res := res + pvol;
! 414: Deep_Clear(pl);
! 415: end if;
! 416: end if;
! 417: end if;
! 418: cont := true;
! 419: end Compute_Facet;
! 420:
! 421: procedure Compute_Facets is new Enumerate_Faces(Compute_Facet);
! 422:
! 423: begin
! 424: res := 0;
! 425: if not Is_In(l,nulvec)
! 426: then fix := Graded_Max(l);
! 427: wl := Shift(l,fix); sh := true;
! 428: else wl := l; sh := false;
! 429: end if;
! 430: Move_to_Front(wl,nulvec);
! 431: pts := Shallow_Create(Tail_Of(wl));
! 432: Compute_Facets(n-1,pts);
! 433: if sh
! 434: then Deep_Clear(wl); Clear(fix);
! 435: end if;
! 436: vlm := res;
! 437: end Sum;
! 438:
! 439: procedure Volume ( n : in natural; l : in List;
! 440: tv : in out Tree_of_Vectors; vlm : out natural ) is
! 441: m : natural;
! 442: begin
! 443: if n = 1
! 444: then declare
! 445: min,max : integer := 0;
! 446: d : Link_to_Vector := Head_Of(l);
! 447: begin
! 448: Min_Max(l,d'first,min,max);
! 449: vlm := max - min;
! 450: end;
! 451: else m := Length_Of(l);
! 452: if m <= n
! 453: then vlm := 0;
! 454: else Sum(n,m,l,tv,vlm);
! 455: end if;
! 456: end if;
! 457: end Volume;
! 458:
! 459: -- MIXED VOLUME COMPUTATIONS WITHOUT TREE OF USEFUL DIRECTIONS :
! 460:
! 461: function Two_Terms_Mixed_Vol ( n : natural; al : Array_of_Lists )
! 462: return natural is
! 463:
! 464: -- DESCRIPTION :
! 465: -- returns the mixed volume of the polytopes generated by the
! 466: -- points in al, where the first polytope is generated by only
! 467: -- two points.
! 468:
! 469: first,second : Link_to_Vector;
! 470: l : List renames al(al'first);
! 471: res : natural;
! 472:
! 473: begin
! 474: first := Head_Of(l);
! 475: second := Head_Of(Tail_Of(l));
! 476: declare
! 477: d : Vector(first'range);
! 478: piv : integer := 0;
! 479: begin
! 480: for i in d'range loop
! 481: d(i) := first(i) - second(i);
! 482: if (piv = 0) and then (d(i) /= 0)
! 483: then piv := i;
! 484: end if;
! 485: end loop;
! 486: if piv = 0
! 487: then return 0;
! 488: else if d(piv) < 0
! 489: then Min(d);
! 490: end if;
! 491: declare
! 492: t : Transfo := Rotate(d,piv);
! 493: tr_al : Array_of_Lists(al'first..(al'last-1));
! 494: degen : boolean := false;
! 495: begin
! 496: Apply(t,d);
! 497: for i in tr_al'range loop
! 498: tr_al(i) := Transform_and_Reduce(t,piv,al(i+1));
! 499: Remove_Duplicates(tr_al(i));
! 500: degen := (Length_Of(tr_al(i)) <= 1);
! 501: exit when degen;
! 502: end loop;
! 503: Clear(t);
! 504: if not degen
! 505: then res := d(piv)*Mixed_Volume(n-1,tr_al);
! 506: else res := 0;
! 507: end if;
! 508: Deep_Clear(tr_al);
! 509: end;
! 510: end if;
! 511: end;
! 512: return res;
! 513: end Two_Terms_Mixed_Vol;
! 514:
! 515: function Facet_Normal ( n : natural; facet,pts : VecVec ) return Vector is
! 516:
! 517: res,pt1,pt2 : Vector(1..n);
! 518: im : Matrix(1..n,1..n);
! 519: cnt : natural := 0;
! 520:
! 521: begin
! 522: for i in facet'range loop
! 523: if facet(i)'length > 1
! 524: then pt1 := pts(facet(i)(facet(i)'first)).all;
! 525: for j in facet(i)'first+1..facet(i)'last loop
! 526: pt2 := pts(facet(i)(j)).all;
! 527: cnt := cnt + 1;
! 528: for k in pt1'range loop
! 529: im(cnt,k) := pt2(k) - pt1(k);
! 530: end loop;
! 531: end loop;
! 532: end if;
! 533: end loop;
! 534: for j in 1..n loop
! 535: im(n,j) := 0;
! 536: end loop;
! 537: Upper_Triangulate(im);
! 538: Scale(im);
! 539: res := (1..n => 0);
! 540: Solve0(im,res);
! 541: Normalize(res);
! 542: -- put("The facet normal : "); put(res); new_line;
! 543: return res;
! 544: end Facet_Normal;
! 545:
! 546: function Minkowski_Sum ( n : natural; al : Array_of_Lists )
! 547: return natural is
! 548:
! 549: -- DESCRIPTION :
! 550: -- Computes the mixed volume of the polytope generated
! 551: -- by the points in al, where n > 1.
! 552:
! 553: res,m,ptslen : natural;
! 554: vl,vl_last,al1 : List;
! 555: typ,ind : Vector(1..n);
! 556: perm,mix : Link_to_Vector;
! 557: wal : Array_of_Lists(al'range) := Interchange2(al);
! 558:
! 559: procedure Update ( v : in Vector; i : in integer;
! 560: added : in out boolean ) is
! 561:
! 562: -- DESCRIPTION :
! 563: -- This procedure computes the support of the first list
! 564: -- in the direction v; if this support is not zero,
! 565: -- the projection will be computed.
! 566: -- The parameter added becomes true if v has been added to vl.
! 567:
! 568: pal : Array_of_Lists(al'first..(al'last-1));
! 569: pv : integer;
! 570: degen : boolean;
! 571:
! 572: begin
! 573: if not Is_In(vl,v)
! 574: then pv := Maximal_Support(al1,v);
! 575: if pv > 0
! 576: then Projection(wal,v,i,pal,degen);
! 577: if not degen
! 578: then declare
! 579: mv : integer := Mixed_Volume(n-1,pal);
! 580: begin
! 581: if mv > 0
! 582: then res := res + pv*mv;
! 583: Append(vl,vl_last,v);
! 584: added := true;
! 585: end if;
! 586: end;
! 587: Deep_Clear(pal);
! 588: end if;
! 589: end if;
! 590: end if;
! 591: end Update;
! 592:
! 593: procedure Enumerate_Facets ( lpts : in Array_of_Lists; len : in natural ) is
! 594:
! 595: pts : VecVec(1..len);
! 596: cnt : integer;
! 597:
! 598: procedure Compute_Facet ( facet : in VecVec; cont : out boolean ) is
! 599:
! 600: v,w : Vector(1..n);
! 601: i : integer;
! 602: added : boolean;
! 603:
! 604: begin
! 605: v := Facet_Normal(n,facet,pts);
! 606: i := Pivot(v);
! 607: if i <= v'last
! 608: then added := false;
! 609: Update(v,i,added);
! 610: if not added
! 611: then Min(v); w := v;
! 612: else w := -v; added := false;
! 613: end if;
! 614: Update(w,i,added);
! 615: end if;
! 616: cont := true;
! 617: end Compute_Facet;
! 618:
! 619: procedure Compute_Facets is new Enumerate_Faces_of_Sum(Compute_Facet);
! 620:
! 621: begin
! 622: pts(ind(1)..ind(2)-1) := Shallow_Create(Tail_Of(al1));
! 623: cnt := lpts'first + mix(mix'first);
! 624: for i in mix'first+1..mix'last loop
! 625: if i < mix'last
! 626: then pts(ind(i)..ind(i+1)-1) := Shallow_Create(lpts(cnt));
! 627: cnt := cnt + mix(i);
! 628: else pts(ind(i)..len) := Shallow_Create(lpts(cnt));
! 629: end if;
! 630: end loop;
! 631: Compute_Facets(ind(mix'range),typ(mix'range),n-1,pts);
! 632: end Enumerate_Facets;
! 633:
! 634: begin
! 635: m := Length_Of(wal(wal'first));
! 636: if m = 2
! 637: then return Two_Terms_Mixed_Vol(n,wal);
! 638: elsif m > 2
! 639: then
! 640: Mixture(al,perm,mix);
! 641: -- put("Type of mixture : "); put(mix); new_line;
! 642: -- put(" with permutation : "); put(perm); new_line;
! 643: wal := Permute(perm.all,al);
! 644: declare
! 645: shiftvec : Link_to_Vector;
! 646: nulvec : Vector(1..n) := (1..n => 0);
! 647: shifted : boolean;
! 648: cnt : integer;
! 649: begin
! 650: -- SHIFT OF THE FIRST LIST ( then all pv >= 0) :
! 651: if not Is_In(wal(wal'first),nulvec)
! 652: then shiftvec := Graded_Max(wal(wal'first));
! 653: al1 := Shift(wal(wal'first),shiftvec); shifted := true;
! 654: else al1 := wal(wal'first); shifted := false;
! 655: end if;
! 656: -- ENUMERATE FACES OF SUM OF THE RIGHT TYPE :
! 657: Move_to_Front(al1,nulvec);
! 658: wal(wal'first) := al1;
! 659: res := 0;
! 660: typ(1) := mix(mix'first)-1;
! 661: typ(2) := mix(mix'first+1); ind(1) := 1;
! 662: ind(2) := ind(1) + Length_Of(al1) - 1; -- skip null vector
! 663: cnt := wal'first + mix(mix'first);
! 664: for i in mix'first+2..mix'last loop
! 665: typ(i) := mix(i);
! 666: ind(i) := ind(i-1) + Length_Of(wal(cnt));
! 667: cnt := cnt + mix(i-1);
! 668: end loop;
! 669: ptslen := ind(mix'last) + Length_Of(wal(cnt)) - 1;
! 670: Enumerate_Facets(wal,ptslen);
! 671: -- CLEANING UP :
! 672: Deep_Clear(vl); Clear(perm); Clear(mix);
! 673: if shifted
! 674: then Clear(shiftvec);
! 675: Deep_Clear(al1);
! 676: end if;
! 677: return res;
! 678: end;
! 679: else -- m < 2
! 680: return 0;
! 681: end if;
! 682: end Minkowski_Sum;
! 683:
! 684: function Mixed_Volume ( n : natural; al : Array_of_Lists )
! 685: return natural is
! 686: begin
! 687: if (n = 0) or Is_Null(al(al'first))
! 688: then return 0;
! 689: elsif n = 1
! 690: then declare
! 691: min,max : integer := 0;
! 692: d : Link_to_Vector := Head_Of(al(al'first));
! 693: begin
! 694: Min_Max(al(al'first),d'first,min,max);
! 695: return (max - min);
! 696: end;
! 697: elsif All_Equal(al)
! 698: then return Volume(n,al(al'first));
! 699: else return Minkowski_Sum(n,al);
! 700: end if;
! 701: end Mixed_Volume;
! 702:
! 703: -- MIXED VOLUME COMPUTATIONS WITH TREE OF USEFUL DIRECTIONS :
! 704:
! 705: function Two_Terms_Mixed_Volume ( n : natural; al : Array_of_Lists;
! 706: tv : Tree_of_Vectors )
! 707: return natural is
! 708:
! 709: -- DESCRIPTION :
! 710: -- returns the mixed volume of the polytopes generated by the
! 711: -- points in al, where the first polytope is generated by only
! 712: -- two points.
! 713:
! 714: first,second : Link_to_Vector;
! 715: l : List renames al(al'first);
! 716: res : natural;
! 717:
! 718: begin
! 719: first := Head_Of(l);
! 720: second := Head_Of(Tail_Of(l));
! 721: declare
! 722: d : Vector(first'range);
! 723: piv : integer := 0;
! 724: begin
! 725: for i in d'range loop
! 726: d(i) := first(i) - second(i);
! 727: if (piv = 0) and then (d(i) /= 0)
! 728: then piv := i;
! 729: end if;
! 730: end loop;
! 731: if piv = 0
! 732: then return 0;
! 733: else if d(piv) < 0
! 734: then Min(d);
! 735: end if;
! 736: declare
! 737: t : Transfo := Rotate(d,piv);
! 738: tr_al : Array_of_Lists(al'first..(al'last-1));
! 739: degen : boolean := false;
! 740: begin
! 741: Apply(t,d);
! 742: for i in tr_al'range loop
! 743: tr_al(i) := Transform_and_Reduce(t,piv,al(i+1));
! 744: Remove_Duplicates(tr_al(i));
! 745: degen := (Length_Of(tr_al(i)) <= 1);
! 746: exit when degen;
! 747: end loop;
! 748: Clear(t);
! 749: if not degen
! 750: then res := d(piv)*Mixed_Volume(n-1,tr_al,tv);
! 751: else res := 0;
! 752: end if;
! 753: Deep_Clear(tr_al);
! 754: end;
! 755: end if;
! 756: end;
! 757: return res;
! 758: end Two_Terms_Mixed_Volume;
! 759:
! 760: function Minkowski_Sum ( n : natural; al : Array_of_Lists;
! 761: tv : Tree_of_Vectors ) return natural is
! 762:
! 763: -- DESCRIPTION :
! 764: -- Computes the mixed volume of the polytope generated
! 765: -- by the points in al, where n > 1.
! 766: -- The tree of degrees is not empty.
! 767:
! 768: res,m : natural;
! 769: al1 : List;
! 770: wal : Array_of_Lists(al'range) := Interchange2(al);
! 771: tmp : Tree_of_Vectors;
! 772: perm,mix : Link_to_Vector;
! 773:
! 774: begin
! 775: m := Length_Of(wal(wal'first));
! 776: if m = 2
! 777: then return Two_Terms_Mixed_Volume(n,wal,tv);
! 778: elsif m > 2
! 779: then
! 780: Mixture(al,perm,mix);
! 781: wal := Permute(perm.all,al);
! 782: declare
! 783: shiftvec : Link_to_Vector;
! 784: nulvec : Vector(1..n) := (1..n => 0);
! 785: shifted : boolean;
! 786: begin
! 787: -- SHIFT OF THE FIRST LIST ( then all pv >= 0) :
! 788: if not Is_In(wal(wal'first),nulvec)
! 789: then shiftvec := Graded_Max(wal(wal'first));
! 790: al1 := Shift(wal(wal'first),shiftvec); shifted := true;
! 791: else al1 := wal(wal'first); shifted := false;
! 792: end if;
! 793: Move_to_Front(al1,nulvec);
! 794: wal(wal'first) := al1;
! 795: -- COMPUTING THE MIXED VOLUME :
! 796: tmp := tv; res := 0;
! 797: while not Is_Null(tmp) loop
! 798: declare
! 799: nd : node := Head_Of(tmp);
! 800: v : Link_to_Vector := nd.d;
! 801: i : integer := Pivot(v);
! 802: pv : integer := Maximal_Support(al1,v.all);
! 803: pal : Array_of_Lists(al'first..(al'last-1));
! 804: degen : boolean;
! 805: begin
! 806: Projection(wal,v.all,i,pal,degen);
! 807: if (nd.ltv = null) or else Is_Null(nd.ltv.all)
! 808: then res := res + pv*Mixed_Volume(n-1,pal);
! 809: else res := res + pv*Mixed_Volume(n-1,pal,nd.ltv.all);
! 810: end if;
! 811: Deep_Clear(pal);
! 812: end;
! 813: tmp := Tail_Of(tmp);
! 814: end loop;
! 815: -- CLEANING UP :
! 816: Clear(perm); Clear(mix);
! 817: if shifted
! 818: then Clear(shiftvec); Deep_Clear(al1);
! 819: end if;
! 820: return res;
! 821: end;
! 822: else -- m < 2
! 823: return 0;
! 824: end if;
! 825: end Minkowski_Sum;
! 826:
! 827: function Mixed_Volume ( n : natural; al : Array_of_Lists;
! 828: tv : Tree_of_Vectors ) return natural is
! 829: begin
! 830: if (n = 0) or Is_Null(al(al'first))
! 831: then return 0;
! 832: elsif n = 1
! 833: then declare
! 834: min,max : integer := 0;
! 835: d : Link_to_Vector := Head_Of(al(al'first));
! 836: begin
! 837: Min_Max(al(al'first),d'first,min,max);
! 838: return (max - min);
! 839: end;
! 840: elsif All_Equal(al)
! 841: then return Volume(n,al(al'first),tv);
! 842: elsif not Is_Null(tv)
! 843: then return Minkowski_Sum(n,al,tv);
! 844: else return Minkowski_Sum(n,al);
! 845: end if;
! 846: end Mixed_Volume;
! 847:
! 848: -- MIXED VOLUME COMPUTATIONS WITH CREATION OF TREE OF USEFUL DIRECTIONS :
! 849:
! 850: procedure Two_Terms_Mixed_Vol
! 851: ( n : in natural; al : in Array_of_Lists;
! 852: tv : in out Tree_of_Vectors; mv : out natural ) is
! 853:
! 854: -- DESCRIPTION :
! 855: -- returns the mixed volume of the polytopes generated by the
! 856: -- points in al, where the first polytope is generated by only
! 857: -- two points.
! 858:
! 859: first,second : Link_to_Vector;
! 860: l : List renames al(al'first);
! 861:
! 862: begin
! 863: first := Head_Of(l);
! 864: second := Head_Of(Tail_Of(l));
! 865: declare
! 866: d : Vector(first'range);
! 867: piv : integer := 0;
! 868: begin
! 869: for i in d'range loop
! 870: d(i) := first(i) - second(i);
! 871: if (piv = 0) and then (d(i) /= 0)
! 872: then piv := i;
! 873: end if;
! 874: end loop;
! 875: if piv = 0
! 876: then mv := 0;
! 877: else if d(piv) < 0
! 878: then Min(d);
! 879: end if;
! 880: declare
! 881: t : Transfo := Rotate(d,piv);
! 882: tr_al : Array_of_Lists(al'first..(al'last-1));
! 883: mvl : natural;
! 884: degen : boolean := false;
! 885: begin
! 886: Apply(t,d);
! 887: for i in tr_al'range loop
! 888: tr_al(i) := Transform_and_Reduce(t,piv,al(i+1));
! 889: Remove_Duplicates(tr_al(i));
! 890: degen := (Length_Of(tr_al(i)) <= 1);
! 891: exit when degen;
! 892: end loop;
! 893: Clear(t);
! 894: if not degen
! 895: then Mixed_Volume(n-1,tr_al,tv,mvl); mv := d(piv)*mvl;
! 896: else mv := 0;
! 897: end if;
! 898: Deep_Clear(tr_al);
! 899: end;
! 900: end if;
! 901: end;
! 902: end Two_Terms_Mixed_Vol;
! 903:
! 904: procedure Minkowski_Sum ( n : in natural; al : in Array_of_Lists;
! 905: tv : in out Tree_of_Vectors; mv : out natural ) is
! 906:
! 907: -- DESCRIPTION :
! 908: -- Computes the mixed volume of the polytope generated
! 909: -- by the points in al, where n > 1.
! 910:
! 911: res,m,ptslen : natural;
! 912: al1 : List;
! 913: typ,ind : Vector(1..n);
! 914: perm,mix : Link_to_Vector;
! 915: wal : Array_of_Lists(al'range) := Interchange2(al);
! 916:
! 917: procedure Update ( v : in Vector; i : in integer;
! 918: added : in out boolean ) is
! 919:
! 920: -- DESCRIPTION :
! 921: -- This procedure computes the support of the first list
! 922: -- in the direction v; if this support is not zero,
! 923: -- the projection will be computed.
! 924: -- The parameter added becomes true if v has been added to vl.
! 925:
! 926: pal : Array_of_Lists(al'first..(al'last-1));
! 927: pv : integer;
! 928: degen : boolean;
! 929:
! 930: begin
! 931: if not Is_In(tv,v)
! 932: then pv := Maximal_Support(al1,v);
! 933: if pv > 0
! 934: then Projection(wal,v,i,pal,degen);
! 935: if not degen
! 936: then declare
! 937: tmp : Tree_of_Vectors;
! 938: mv : natural;
! 939: begin
! 940: Mixed_Volume(n-1,pal,tmp,mv);
! 941: if mv = 0
! 942: then Clear(tmp);
! 943: else res := res + pv*mv;
! 944: declare
! 945: nd : node;
! 946: begin
! 947: nd.d := new Vector'(v);
! 948: if not Is_Null(tmp)
! 949: then nd.ltv := new Tree_of_Vectors'(tmp);
! 950: end if;
! 951: Construct(nd,tv);
! 952: end;
! 953: added := true;
! 954: end if;
! 955: end;
! 956: Deep_Clear(pal);
! 957: end if;
! 958: end if;
! 959: end if;
! 960: end Update;
! 961:
! 962: procedure Enumerate_Facets ( lpts : in Array_of_Lists; len : in natural ) is
! 963:
! 964: pts : VecVec(1..len);
! 965: cnt : integer;
! 966:
! 967: procedure Compute_Facet ( facet : in VecVec; cont : out boolean ) is
! 968:
! 969: v,w : Vector(1..n);
! 970: i : integer;
! 971: added : boolean;
! 972: begin
! 973: v := Facet_Normal(n,facet,pts);
! 974: i := Pivot(v);
! 975: if i <= v'last
! 976: then added := false;
! 977: Update(v,i,added);
! 978: if not added
! 979: then Min(v); w := v;
! 980: else w := -v; added := false;
! 981: end if;
! 982: Update(w,i,added);
! 983: end if;
! 984: cont := true;
! 985: end Compute_Facet;
! 986:
! 987: procedure Compute_Facets is
! 988: new Enumerate_Faces_of_Sum(Compute_Facet);
! 989:
! 990: begin
! 991: pts(ind(1)..ind(2)-1) := Shallow_Create(Tail_Of(al1));
! 992: cnt := lpts'first + mix(mix'first);
! 993: for i in mix'first+1..mix'last loop
! 994: if i < mix'last
! 995: then pts(ind(i)..ind(i+1)-1) := Shallow_Create(lpts(cnt));
! 996: cnt := cnt + mix(i);
! 997: else pts(ind(i)..len) := Shallow_Create(lpts(cnt));
! 998: end if;
! 999: end loop;
! 1000: Compute_Facets(ind(mix'range),typ(mix'range),n-1,pts);
! 1001: end Enumerate_Facets;
! 1002:
! 1003: begin
! 1004: m := Length_Of(wal(wal'first));
! 1005: if m = 2
! 1006: then Two_Terms_Mixed_Vol(n,wal,tv,mv);
! 1007: elsif m > 2
! 1008: then
! 1009: Mixture(al,perm,mix);
! 1010: wal := Permute(perm.all,al);
! 1011: declare
! 1012: shiftvec : Link_to_Vector;
! 1013: nulvec : Vector(1..n) := (1..n => 0);
! 1014: shifted : boolean;
! 1015: cnt : integer;
! 1016: begin
! 1017: -- SHIFT OF THE FIRST LIST ( then all pv >= 0) :
! 1018: if not Is_In(wal(wal'first),nulvec)
! 1019: then shiftvec := Graded_Max(wal(wal'first));
! 1020: al1 := Shift(wal(wal'first),shiftvec); shifted := true;
! 1021: else al1 := wal(wal'first); shifted := false;
! 1022: end if;
! 1023: -- ENUMERATE FACES OF SUM OF THE RIGHT TYPE :
! 1024: Move_to_Front(al1,nulvec);
! 1025: wal(wal'first) := al1;
! 1026: res := 0;
! 1027: typ(1) := mix(mix'first)-1;
! 1028: typ(2) := mix(mix'first+1); ind(1) := 1;
! 1029: ind(2) := ind(1) + Length_Of(al1) - 1; -- skip null vector
! 1030: cnt := wal'first + mix(mix'first);
! 1031: for i in mix'first+2..mix'last loop
! 1032: typ(i) := mix(i);
! 1033: ind(i) := ind(i-1) + Length_Of(wal(cnt));
! 1034: cnt := cnt + mix(i-1);
! 1035: end loop;
! 1036: ptslen := ind(mix'last) + Length_Of(wal(cnt)) - 1;
! 1037: Enumerate_Facets(wal,ptslen);
! 1038: -- CLEANING UP :
! 1039: Clear(perm); Clear(mix);
! 1040: if shifted
! 1041: then Clear(shiftvec); Deep_Clear(al1);
! 1042: end if;
! 1043: mv := res;
! 1044: end;
! 1045: else -- m < 2
! 1046: mv := 0;
! 1047: end if;
! 1048: end Minkowski_Sum;
! 1049:
! 1050: procedure Mixed_Volume ( n : in natural; al : in Array_of_Lists;
! 1051: tv : in out Tree_of_Vectors; mv : out natural ) is
! 1052: begin
! 1053: if (n = 0) or Is_Null(al(al'first))
! 1054: then mv := 0;
! 1055: elsif n = 1
! 1056: then declare
! 1057: min,max : integer := 0;
! 1058: d : Link_to_Vector := Head_Of(al(al'first));
! 1059: begin
! 1060: Min_Max(al(al'first),d'first,min,max);
! 1061: mv := max - min;
! 1062: end;
! 1063: elsif All_Equal(al)
! 1064: then Volume(n,al(al'first),tv,mv);
! 1065: else Minkowski_Sum(n,al,tv,mv);
! 1066: end if;
! 1067: end Mixed_Volume;
! 1068:
! 1069: end Volumes;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to volumes.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 / Implift