![[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|
Return to volumes.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift |
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;