source: git/Singular/LIB/gitfan.lib @ ea6e16e

spielwiese
Last change on this file since ea6e16e was ea6e16e, checked in by Sachin <sachinkm308@…>, 4 years ago
replacing execute with create_ring(16)
  • Property mode set to 100644
File size: 88.0 KB
Line 
1////////////////////////////////////////////////////////////////////
2version="version gitfan.lib 4.1.2.0 Feb_2019 ";
3category="Algebraic Geometry";
4info="
5LIBRARY:   gitfan.lib       Compute GIT-fans.
6
7AUTHORS: Janko Boehm, boehm at mathematik.uni-kl.de @*
8         Simon Keicher, keicher at mail.mathematik.uni-tuebingen.de @*
9         Yue Ren, ren at mathematik.uni-kl.de @*
10
11OVERVIEW:
12This library allows you to calculate GIT-fans, torus orbits and GKZ-fans.
13
14In provides features to make use of symmetries of the torus action under consideration.
15
16The main procedure is GITfan which can be directly applied to an ideal and a grading matrix
17encoding the torus action, and returns a fan, the associated GIT-fan.
18We also provide various procedures implementing substeps of the algorithm
19to deal with large computations.
20
21The library uses the package 'gfanlib' by Anders N. Jensen.
22
23For notation, background, and algorithms see [BKR16].
24
25Functions produce debug output if printlevel is positive.
26
27Elements of the symmetric group Sn of type permutation can be created by
28the function permutationFromIntvec.
29The images of 1,...,n can be obtained by permutationToIntvec.
30Composition of permutations can be done by the *-Operator,
31also powers can be computed in the usual way.
32
33
34REFERENCES:
35
36[BKR16] J. Boehm, S. Keicher, Y. Ren: Computing GIT-Fans with Symmetry and the Mori Chamber Decomposition of M06bar,  https://arxiv.org/abs/1603.09241
37
38TYPES:
39  permutation; Permutation in map representation.
40
41PROCEDURES:
42  isAface(ideal,intvec); Checks whether the given face is an a-face.
43  afaces(ideal); Returns a list of intvecs that correspond to the set of all a-faces, optionally for given list of simplex faces.
44  fullDimImages(list,intmat); Finds the afaces which have a full-dimensional projection.
45  minimalAfaces(list); compute the minimal a-faces among the a-faces with full dimensional projection.
46  orbitCones(list,intmat); Returns the list of all orbit cones.
47  GITcone(list,bigintmat); Returns the GIT-cone containing the given weight vector.
48  GITfan(ideal,intmat); Compute GIT-fan.
49  GITfanFromOrbitCones(list,intmat,cone); Compute GIT-fan from orbit cones.
50  GITfanParallel(list,intmat,cone); Compute GIT-fan in parallel from orbit cones.
51  GKZfan(intmat); Returns the GKZ-fan of the matrix Q.
52  movingCone(intmat); Compute the moving cone.
53  computeAfaceOrbits(list,list); Compute orbits of a-faces under a permutation group action.
54  minimalAfaceOrbits(list); Find the minimal a-face orbits.
55  orbitConeOrbits(list,intmat); Project the a-face orbits to orbit cone orbits.
56  minimalOrbitConeOrbits(list); Find the minimal orbit cone orbits.
57  intersectOrbitsWithMovingCone(list,cone); Intersect orbit cone orbits with moving cone.
58  groupActionOnQImage(list,intmat); Determine the induced group action in the target of the grading matrix.
59  groupActionOnHashes(list,list); Determine the induced group action on the set of orbit cones.
60  storeActionOnOrbitConeIndices(list, string); Write the group action on the set of orbit cones to a file.
61  permutationFromIntvec(intvec); Create a permutation from an intvec of images.
62  permutationToIntvec(permutation); Return the intvec of images.
63  evaluateProduct(list,list); Evaluate a list of products of group elements in terms of a given representation of the elements as permutations.
64  GITfanSymmetric(list,intmat,cone,list); Compute GIT-fan from orbit cones by determining a minimal representing set for the orbits of maximal dimensional GIT-cones.
65  GITfanParallelSymmetric(list, intmat, cone, list); Compute GIT-fan in parallel from orbit cones by determining a minimal representing set for the orbits of maximal dimensional GIT-cones.
66  bigintToBinary(bigint,int); Convert bigint into a sparse binary represenation specifying the indices of the one-entries
67  binaryToBigint(intvec); Convert sparse binary represenation specifying the indices of the one-entries to bigint
68  applyPermutationToIntvec(intvec,permutation); Apply permutation to a set of integers represented as an intvec
69  hashToCone(bigint,list); Convert a bigint hash to a GIT-cone
70  hashesToFan(list hashes,list OC)
71gitCone(ideal,bigintmat,bigintmat); Returns the GIT-cone around the given weight vector w
72
73KEYWORDS: library, gitfan, GIT, geometric invariant theory, quotients
74";
75
76
77LIB "customstd.lib";
78LIB "linalg.lib";
79LIB "multigrading.lib";
80LIB "parallel.lib";
81
82static proc mod_init()
83{
84  LIB "gfanlib.so";
85  LIB "gitfan.so";
86  //LIB "subsets.so";
87  newstruct("permutation","intvec image");
88  system("install","permutation","*",composePermutations,2);
89  system("install","permutation","^",permutationPower,2);
90  system("install","permutation","print",printPermutation,1);
91}
92
93static proc emptyString(int n)
94{
95string st;
96for (int i = 1; i<=n;i++){
97  st=st+" ";
98}
99return(st);}
100
101proc printPermutation(permutation sigma)
102{
103 intvec v = permutationToIntvec(sigma);
104 string vsrc,vimg;
105 for (int i = 1; i<=size(v);i++){
106   vsrc=vsrc+emptyString(1+size(string(size(v)))-size(string(i)))+string(i);
107 }
108 for (i = 1; i<=size(v);i++){
109   vimg=vimg+emptyString(1+size(string(size(v)))-size(string(v[i])))+string(v[i]);
110 }
111 print("|"+vsrc+"|");
112 print("|"+vimg+"|");
113}
114
115
116////////////////////////////////////////////////////
117// converts e.g. n=5 to its binary representation, i.e. 0,1,0,1
118// if r = 3.
119// and stores it in an intvec.
120// r gives the bound for n <= 2^r:
121static proc int2face(int n, int r)
122{
123  int k = r-1;
124  intvec v;
125  int n0 = n;
126
127  while(n0 > 0){
128    while(2^k > n0){
129      k--;
130      //v[size(v)+1] = 0;
131    }
132
133    v = k+1,v;
134    n0 = n0 - 2^k;
135    k--;
136  }
137  v = v[1..size(v)-1];
138  return(v);
139}
140example
141{
142  echo = 2;
143  int n = 5;
144  int r = 4;
145  int2face(n, r);
146
147  n = 1;
148  r = 1;
149  int2face(n,r);
150}
151
152////////
153
154
155
156////////////////////////////////////////////////////
157
158proc afaces(ideal a, list #)
159"USAGE: afaces(a [,L]); a: ideal, L: list of intvecs
160PURPOSE: Returns a list of all a-faces (considered as intvecs of 0 and 1, where the i-th entry is 1 if the cone has the i-th unit basis vector as a generator), if L is specified only the faces of the simplex listed in L are considered (e.g. representatives with respect to a group action).
161RETURN: a list of intvecs
162EXAMPLE: example afaces; shows an example
163"
164{
165  list AF;
166  if ((size(#)>0) and (typeof(#[1])=="intvec")){
167    list L = #;
168    for(int i = 1; i<=size(L); i++ ){
169      dbprint("representative "+string(i)+" of "+string(size(L)));
170      if (isAface(a,L[i])){
171        AF[size(AF) + 1] = L[i];
172      }
173    }
174    return(AF);
175  }
176
177  intvec gam0;
178  int r = nvars(basering);
179
180  // check if 0 is an a-face:
181  int bd;
182  if (size(#)>0){bd=#[1];}
183  gam0 = 0;
184  if (size(gam0)>=bd){
185    if (isAface(a,gam0)){
186      AF[size(AF) + 1] = gam0;
187    }
188  }
189  // check for other a-faces:
190  for(int i = 1; i < 2^r; i++ ){
191    gam0 = int2face(i,r);
192    if (size(gam0)>=bd){
193      if (isAface(a,gam0)){
194        AF[size(AF) + 1] = gam0;
195      }
196    }
197  }
198
199  //"done checking a-faces!";
200  return(AF);
201}
202example
203{
204
205  echo = 2;
206  ring R = 0,T(1..3),dp;
207  ideal a = T(1)+T(2)+T(3);
208
209  list F = afaces(a);
210  print(F);
211  print(size(F));
212
213  // 2nd ex //
214  ring R2 = 0,T(1..3),dp;
215  ideal a2 = T(2)^2*T(3)^2+T(1)*T(3);
216
217  list F2 = afaces(a2);
218  print(F2);
219  print(size(F2));
220
221  // 3rd ex //
222  ring R3 = 0,T(1..3),dp;
223  ideal a3 = 0;
224
225  list F3 = afaces(a3);
226  print(F3);
227  print(size(F3));
228
229
230  // 4th ex //
231  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
232  ideal J =
233    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
234    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
235    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
236    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
237    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
238  list F4 = afaces(J);
239  print(size(F4));
240}
241
242
243
244static proc saturateWithRespectToVariable(ideal I, int k)
245{
246  // "saturating with respect to variable "+string(k);
247  ASSUME(1,k>=1);
248  ASSUME(1,k<=nvars(basering));
249
250  def origin = basering;
251  int n = nvars(basering);
252  intvec weightVector = ringlist(origin)[3][1][2];
253
254  string newVars;
255  intvec newWeightVector;
256  for (int i=1; i<k; i++)
257  {
258    newVars = newVars+string(var(i))+",";
259    newWeightVector[i]=weightVector[i];
260  }
261  for (i=k+1; i<=n; i++)
262  {
263    newVars = newVars+string(var(i))+",";
264    newWeightVector[i-1]=weightVector[i];
265  }
266  newVars = newVars+string(var(k));
267  newWeightVector[n]=weightVector[k];
268  ring ringForSaturation = create_ring(ringlist(origin)[1], "("+newVars+")",
269      list(list("wp", newWeightVector)),"no_minpoly");
270  ideal I = satstd(imap(origin,I));
271  I = simplify(I,2+4+32);
272
273  // "finished saturating with respect to variable "+string(k);
274  setring origin;
275  return (imap(ringForSaturation,I));
276}
277
278static proc stepwiseSaturation(ideal I)
279{
280  if (I!=1)
281  {
282    list variablesToBeSaturated;
283    int l = nvars(basering);
284    for (int i=1; i<=l; i++)
285    { variablesToBeSaturated[i]=l-i+1; }
286
287    while (size(variablesToBeSaturated)>0)
288    {
289      I = saturateWithRespectToVariable(I,variablesToBeSaturated[1]);
290      variablesToBeSaturated = delete(variablesToBeSaturated,1);
291      if ((I==1) || (I==-1))
292      {
293        break;
294      }
295    }
296  }
297
298  return (I);
299}
300
301
302
303
304
305
306proc isAface(ideal a, intvec gam0)
307  "USAGE: isAface(a,gam0); a: ideal, gam0:intvec
308PURPOSE: Checks whether gam0 is an a-face w.r.t. the ideal a.
309RETURN: int
310EXAMPLE: example isaface; shows an example
311"
312{
313  // special case: gam0 is the zero-cone:
314  if (size(gam0) == 1 and gam0[1] == 0){
315    poly pz;
316    ideal G;
317    int i;
318    for (int k = 1; k <= size(a); k++) {
319      pz = subst(a[k], var(1), 0);
320      for (i = 2; i <= nvars(basering); i++) {
321        pz = subst(pz, var(i), 0);
322      }
323      G = G, pz;
324    }
325    G = std(G);
326    // monomial inside?:
327    if(G == 1){
328      return(0);
329    }
330    return(1);
331  }
332
333  string initNewRing = "(";
334
335  intvec w = ringlist(basering)[3][1][2];
336  intvec w0;
337  for (int i=1; i<size(gam0); i++)
338  {
339    // take only entries i of w0 with i in gam0
340    initNewRing = initNewRing + string(var(gam0[i])) + ",";
341    w0[i] = w[gam0[i]];
342  }
343  w0[size(gam0)] = w[gam0[size(gam0)]];
344
345  initNewRing = initNewRing + string(var(gam0[size(gam0)])) + ")";
346  def R = basering;
347//  execute(initNewRing);
348  initNewRing;
349  ring Rgam0 = create_ring(0, initNewRing, "wp("+string(w0)+")");
350  ringlist(Rgam0);
351
352  ideal agam0 = imap(R,a);
353
354  ideal G = stepwiseSaturation(agam0);
355
356  if(G == 1)
357  {
358    return(0); // true
359  }
360  return(1); // false
361}
362example
363{
364  echo = 2;
365
366  ring R = 0,(T(1..4)),dp;
367  ideal I = T(1)*T(2)-T(4);
368
369  intvec w = 1,4;
370  intvec v = 1,2,4;
371
372  isAface(I,w); // should be 0
373  isAface(I,v); // should be 1
374}
375
376
377proc applyPermutationToIntvec(intvec v, permutation g)
378"USAGE: applyPermutationToIntvec(v,g); v: intvec, g:permutation
379PURPOSE: Apply g to the set of entries of v. The result is a sorted intvec. We assume that the entries of v are valid arguments of g. We do not require that the input is sorted.
380RETURN: intvec.
381EXAMPLE: example applyPermutationToIntvec; shows an example
382"
383{
384  intvec gv = composeIntvecs(permutationToIntvec(g),v);
385  //for (int i=1;i<=size(v);i++){
386  //  gv[i]=g[v[i]];
387  //}
388  return(sort(gv)[1]);
389}
390example
391{
392permutation g = permutationFromIntvec(intvec(10, 9, 7, 4, 8, 6, 3, 5, 2, 1));
393applyPermutationToIntvec(intvec(1, 3, 4, 6, 8),g);
394}
395
396static proc isContained(intvec v, list orbit){
397  for (int i=1;i<=size(orbit);i++){
398    if (v==orbit[i]){return(1);}
399  }
400  return(0);
401}
402
403
404static proc computeSimplexOrbit(intvec v, list G){
405  list orbit;
406  intvec gv;
407  for (int i=1;i<=size(G);i++){
408    gv=applyPermutationToIntvec(v,G[i]);
409    if (isContained(gv,orbit)==0){orbit[size(orbit)+1]=gv;}
410  }
411  return(orbit);
412}
413
414
415proc computeAfaceOrbits(list AF, list G)
416  "USAGE: computeAfaceOrbits(AF,G); AF list of intvecs, G: list of permutations
417PURPOSE: Computes the orbits of the afaces in the list AF under the group action in G, where G is a list of permutations. We assume that the elements of G form a group and the first entry corresponds to the neutral element.
418RETURN: a list of lists of intvecs
419EXAMPLE: example computeAfaceOrbits; shows an example
420"
421{
422  list listorbits;
423  for (int i=1;i<=size(AF);i++){
424    dbprint(i);
425    listorbits[i]=computeSimplexOrbit(AF[i],G);
426  }
427  return(listorbits);}
428example
429{
430  echo = 2;
431  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
432  ideal J =
433    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
434    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
435    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
436    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
437    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
438  intmat Q[5][10] =
439    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
440    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
441    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
442    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
443    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
444  list simplexSymmetryGroup = G25Action();
445  list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ),
446  intvec( 1, 2, 3, 5, 6 ),
447  intvec( 1, 2, 3, 5, 7 ),
448  intvec( 1, 2, 3, 5, 10 ),
449  intvec( 1, 2, 3, 7, 9 ),
450  intvec( 1, 2, 6, 9, 10 ),
451  intvec( 1, 2, 3, 4, 5, 6 ),
452  intvec( 1, 2, 3, 4, 5, 10 ),
453  intvec( 1, 2, 3, 5, 6, 8 ),
454  intvec( 1, 2, 3, 5, 6, 9 ),
455  intvec( 1, 2, 3, 5, 7, 10 ),
456  intvec( 1, 2, 3, 7, 9, 10 ),
457  intvec( 1, 2, 3, 4, 5, 6, 7 ),
458  intvec( 1, 2, 3, 4, 5, 6, 8 ),
459  intvec( 1, 2, 3, 4, 5, 6, 9 ),
460  intvec( 1, 2, 3, 5, 6, 9, 10 ),
461  intvec( 1, 2, 3, 4, 5, 6, 7, 8 ),
462  intvec( 1, 2, 3, 4, 5, 6, 9, 10 ),
463  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ),
464  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 );
465  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
466  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
467  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
468  apply(afaceOrbits,size);}
469
470
471
472
473static proc isSubset(intvec v, intvec w)
474{
475  int i,j;
476  int jStart=1;
477  for (i=1; i<=size(v); i++)
478  {
479    for (j=jStart; j<=size(w); j++)
480    {
481      if (v[i]==w[j])
482      {
483        break;
484      }
485    }
486    if (j<=size(w))
487    {
488      jStart = j+1;
489    }
490    else
491    {
492      return (0);
493    }
494  }
495  return (1);
496}
497
498
499
500static proc containsOrbitd(list A, list B)
501{
502  intvec a = A[1];
503  for (int i=1; i<=size(B); i++)
504  {
505    if (isSubset(a,B[i]))
506    {
507      return (1);
508    }
509  }
510  return (0);
511}
512
513
514proc minimalAfaceOrbits(list listOfOrbits)
515"USAGE: minimalAfaceOrbits(afaceOrbits); afaceOrbits: list
516PURPOSE: Returns a list of all minimal a-face orbits.
517RETURN: a list of intvecs
518EXAMPLE: example minimalAfaceOrbits; shows an example
519"
520{
521  int i,j;
522  list L = listOfOrbits;
523  for (i=1; i<=size(L); i++)
524  {
525    dbprint(string(i)+" of "+string(size(L)));
526    for (j=1; j<=size(L); j++)
527    {
528      if (i!=j)
529      {
530        if (containsOrbitd(L[i],L[j]))
531        {
532          dbprint("werfe raus (Laenge):  " + string(size(L[j])));
533          L = delete(L,j);
534          if (j<i)
535          {
536            i = i-1;
537          }
538          j = j-1;
539        }
540      }
541    }
542  }
543  return(L);
544}
545example
546{
547  echo = 2;
548  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
549  ideal J =
550    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
551    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
552    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
553    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
554    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
555  intmat Q[5][10] =
556    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
557    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
558    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
559    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
560    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
561  list simplexSymmetryGroup = G25Action();
562  list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ),
563  intvec( 1, 2, 3, 5, 6 ),
564  intvec( 1, 2, 3, 5, 7 ),
565  intvec( 1, 2, 3, 5, 10 ),
566  intvec( 1, 2, 3, 7, 9 ),
567  intvec( 1, 2, 6, 9, 10 ),
568  intvec( 1, 2, 3, 4, 5, 6 ),
569  intvec( 1, 2, 3, 4, 5, 10 ),
570  intvec( 1, 2, 3, 5, 6, 8 ),
571  intvec( 1, 2, 3, 5, 6, 9 ),
572  intvec( 1, 2, 3, 5, 7, 10 ),
573  intvec( 1, 2, 3, 7, 9, 10 ),
574  intvec( 1, 2, 3, 4, 5, 6, 7 ),
575  intvec( 1, 2, 3, 4, 5, 6, 8 ),
576  intvec( 1, 2, 3, 4, 5, 6, 9 ),
577  intvec( 1, 2, 3, 5, 6, 9, 10 ),
578  intvec( 1, 2, 3, 4, 5, 6, 7, 8 ),
579  intvec( 1, 2, 3, 4, 5, 6, 9, 10 ),
580  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ),
581  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 );
582  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
583  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
584  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
585  apply(afaceOrbits,size);
586  list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits);
587  apply(minAfaceOrbits,size);
588}
589
590
591static proc list2string(list L){
592  string s = "";
593
594  for(int i = 1; i <=size(L); i++){
595    s = s + string(size(L[i])) + ", ";
596  }
597
598  return(s);
599}
600
601
602
603proc fullDimImages(list afaces,intmat Q)
604"USAGE: fullDimImages(afaces, Q); afaces: list, Q: intmat
605PURPOSE: Determines the a-faces (represented as intvecs) from the list afaces which have a full-dimensional projection with respect to Q.
606RETURN: a list of intvecs
607EXAMPLE: example fullDimImages; shows an example
608"
609{
610  list L;
611  for (int i=1; i<=size(afaces); i++)
612  {
613          intvec gam0 = afaces[i];
614          intmat Qgam0[nrows(Q)][size(gam0)] = Q[intvec(1..nrows(Q)),gam0];
615          //cone c = coneViaPoints(Qgam0);
616          if(rank(Qgam0) == nrows(Q)){
617      L[size(L)+1]=afaces[i];
618    }
619    kill gam0,Qgam0;
620  }
621  return(L);
622}
623example
624{
625  echo=2;
626  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
627  ideal J =
628    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
629    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
630    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
631    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
632    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
633  intmat Q[5][10] =
634    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
635    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
636    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
637    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
638    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
639
640  list AF= afaces(J,nrows(Q));
641  size(AF);
642  size(fullDimImages(AF,Q));
643}
644
645
646
647
648proc orbitConeOrbits(list F, intmat Q)
649"USAGE: orbitConeOrbits(F, Q); F: list, Q: intmat
650PURPOSE: Projects a list F of a-face orbits to the orbit cones with respect to Q. The function checks whether the projections are of full dimension and returns an error otherwise.
651RETURN: a list of lists of cones
652EXAMPLE: example orbitConeOrbits; shows an example
653"
654{
655  list OC;
656  int j;
657  intmat Qt = transpose(Q);
658  for(int i = 1; i <= size(F); i++){
659    list Orbit = F[i];
660    list U;
661    for(j = 1; j <= size(Orbit); j++){
662      intvec aface = Orbit[j];
663      intmat Qgam0[size(aface)][ncols(Qt)] = Qt[aface,intvec(1..ncols(Qt))];
664      cone c = coneViaPoints(Qgam0);
665      if (dimension(c)!=ncols(Qt)){ERROR("cone should have full dimension");}
666      // insert it even if it is already inside (mincones will take care of this)
667      U[size(U)+1] = c;
668      kill aface;
669      kill Qgam0;
670      kill c;
671    }
672    OC[size(OC)+1] = U;
673    dbprint("current size of OC:");
674    dbprint(size(OC));
675    kill U;
676    kill Orbit;
677  }
678  return(OC);
679}
680example
681{
682  echo = 2;
683  // Note that simplexOrbitRepresentatives and simplexSymmetryGroup are subsets of the actual sets for G25. For the full example see the examples in the documentation
684  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
685  ideal J =
686    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
687    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
688    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
689    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
690    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
691  intmat Q[5][10] =
692    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
693    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
694    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
695    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
696    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
697  list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ),
698    intvec( 1, 2, 3, 5, 6 ),
699    intvec( 1, 2, 3, 5, 7 ),
700    intvec( 1, 2, 3, 5, 10 ),
701    intvec( 1, 2, 3, 7, 9 ),
702    intvec( 1, 2, 3, 4, 5, 6, 9, 10 ),
703    intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ),
704    intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 );
705  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
706  list simplexSymmetryGroup = permutationFromIntvec(intvec( 1 .. 10 )),
707    permutationFromIntvec(intvec( 1, 2, 4, 3, 5, 7, 6, 9, 8, 10 )),
708    permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )),
709    permutationFromIntvec(intvec( 1, 3, 4, 2, 6, 7, 5, 10, 8, 9 )),
710    permutationFromIntvec(intvec( 1, 4, 2, 3, 7, 5, 6, 9, 10, 8 )),
711    permutationFromIntvec(intvec( 1, 4, 3, 2, 7, 6, 5, 10, 9, 8 ));
712  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
713  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
714  apply(afaceOrbits,size);
715  list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits);
716  apply(minAfaceOrbits,size);
717  list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q);
718}
719
720
721proc minimalOrbitConeOrbits(list listOfOrbits)
722"USAGE: minimalOrbitConeOrbits(listOfOrbits); listOfOrbits: list of lists of cones
723PURPOSE: Minimizes a list of orbit cone orbits.
724RETURN: a list of lists of cones
725EXAMPLE: example minimalOrbitConeOrbits; shows an example
726"
727{
728  int i,j;
729  list L = listOfOrbits;
730
731
732  for (i=1; i<=size(L); i++)
733  {
734    dbprint("::::: " + string(i)+" of "+string(size(L)));
735
736    for (j=1; j<=size(L); j++)
737    {
738      dbprint("outer: "+string(i)+"  inner: " + string(j));
739      if (i!=j)
740      {
741        if (containsOrbitdCone(L[i],L[j]))
742        {
743          dbprint("werfe raus, i="+string(i)+",  j="+string(j)+" (Laenge):  " + string(size(L[j])));
744          L = delete(L,j);
745          if (j<i)
746          {
747            i = i-1;
748          }
749          j = j-1;
750        }
751      }
752    }
753  }
754  return(L);
755}
756example
757{
758  echo = 2;
759  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
760  ideal J =
761    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
762    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
763    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
764    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
765    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
766  intmat Q[5][10] =
767    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
768    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
769    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
770    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
771    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
772  list simplexSymmetryGroup = G25Action();
773  list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ),
774  intvec( 1, 2, 3, 5, 6 ),
775  intvec( 1, 2, 3, 5, 7 ),
776  intvec( 1, 2, 3, 5, 10 ),
777  intvec( 1, 2, 3, 7, 9 ),
778  intvec( 1, 2, 6, 9, 10 ),
779  intvec( 1, 2, 3, 4, 5, 6 ),
780  intvec( 1, 2, 3, 4, 5, 10 ),
781  intvec( 1, 2, 3, 5, 6, 8 ),
782  intvec( 1, 2, 3, 5, 6, 9 ),
783  intvec( 1, 2, 3, 5, 7, 10 ),
784  intvec( 1, 2, 3, 7, 9, 10 ),
785  intvec( 1, 2, 3, 4, 5, 6, 7 ),
786  intvec( 1, 2, 3, 4, 5, 6, 8 ),
787  intvec( 1, 2, 3, 4, 5, 6, 9 ),
788  intvec( 1, 2, 3, 5, 6, 9, 10 ),
789  intvec( 1, 2, 3, 4, 5, 6, 7, 8 ),
790  intvec( 1, 2, 3, 4, 5, 6, 9, 10 ),
791  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ),
792  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 );
793  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
794  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
795  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
796  apply(afaceOrbits,size);
797  list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits);
798  apply(minAfaceOrbits,size);
799  list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q);
800  apply(listOfOrbitConeOrbits,size);
801  list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits);
802  size(listOfMinimalOrbitConeOrbits);
803}
804
805
806static proc containsOrbitdCone(list A, list B)
807{
808  cone a = A[1];
809  for (int i=1; i<=size(B); i++)
810  {
811    if (isSubcone(a,B[i]))
812    {
813      dbprint("found subcone i="+string(i));
814      return (1);
815    }
816  }
817  return (0);
818}
819
820static proc isSubcone(cone c1, cone c2)
821{
822  return(containsInSupport(c2, c1));
823}
824
825
826
827
828
829
830proc intersectOrbitsWithMovingCone(list OCmin,cone mov)
831"USAGE: intersectOrbitsWithMovingCone(OCmin, mov); OCmin: list of lists of cones, mov: cone
832PURPOSE: Intersects all cones in the orbits in OCmin with mov discarting all orbits of cones which are not of full dimension. The resulting orbits are duplicate free.
833RETURN: a list of lists of cones
834EXAMPLE: example intersectOrbitsWithMovingCone; shows an example
835"
836{
837  list OCmov;
838  cone ocmov;
839  int i,j;
840  int fulldim=dimension(mov);
841  for (i=1; i<=size(OCmin); i++)
842  {
843    dbprint("intersecting orbit "+string(i)+" with moving cone");
844    ocmov = convexIntersection(OCmin[i][1],mov);
845    if (dimension(ocmov)==fulldim)
846    {
847      list currentOrbit;
848      for (j=1; j<=size(OCmin[i]); j++)
849      {
850        dbprint("checking cone "+string(j)+" of "+string(size(OCmin[i])));
851        ocmov = convexIntersection(OCmin[i][j],mov);
852        ocmov = canonicalizeCone(ocmov);
853        if (!listContainsCone(currentOrbit,ocmov))
854        {
855          dbprint("inserting cone");
856          currentOrbit[size(currentOrbit)+1] = ocmov;
857        }
858      }
859      OCmov[size(OCmov)+1] = currentOrbit;
860      kill currentOrbit;
861    }
862    else
863    {
864      dbprint("intersection with moving cone lower-dimensional, abandoning orbit");
865    }
866  }
867  return(OCmov);
868}
869example {
870  echo=2;
871  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
872  ideal J =
873    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
874    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
875    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
876    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
877    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
878  intmat Q[5][10] =
879    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
880    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
881    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
882    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
883    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
884  list simplexSymmetryGroup = G25Action();
885  list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ),
886  intvec( 1, 2, 3, 5, 6 ),
887  intvec( 1, 2, 3, 5, 7 ),
888  intvec( 1, 2, 3, 5, 10 ),
889  intvec( 1, 2, 3, 7, 9 ),
890  intvec( 1, 2, 6, 9, 10 ),
891  intvec( 1, 2, 3, 4, 5, 6 ),
892  intvec( 1, 2, 3, 4, 5, 10 ),
893  intvec( 1, 2, 3, 5, 6, 8 ),
894  intvec( 1, 2, 3, 5, 6, 9 ),
895  intvec( 1, 2, 3, 5, 7, 10 ),
896  intvec( 1, 2, 3, 7, 9, 10 ),
897  intvec( 1, 2, 3, 4, 5, 6, 7 ),
898  intvec( 1, 2, 3, 4, 5, 6, 8 ),
899  intvec( 1, 2, 3, 4, 5, 6, 9 ),
900  intvec( 1, 2, 3, 5, 6, 9, 10 ),
901  intvec( 1, 2, 3, 4, 5, 6, 7, 8 ),
902  intvec( 1, 2, 3, 4, 5, 6, 9, 10 ),
903  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ),
904  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 );
905  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
906  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
907  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
908  apply(afaceOrbits,size);
909  list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits);
910  apply(minAfaceOrbits,size);
911  list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q);
912  cone mov = movingCone(Q);
913  intersectOrbitsWithMovingCone(listOfOrbitConeOrbits,mov);
914}
915
916
917static proc coneToString(cone C){
918  bigintmat F = facets(C);
919  ring R=0,(x),dp;
920  matrix FF[nrows(F)][ncols(F)];
921  int i,j;
922  for (i=1; i<=nrows(F);i++){
923    for (j=1; j<=ncols(F);j++){
924      FF[i,j]=number(F[i,j]);
925    }
926  }
927  return("bigintmat F["+string(nrows(F))+"]["+string(ncols(F))+"]="+string(FF)+";cone C=canonicalizeCone(coneViaInequalities(F));");
928}
929
930proc storeOrbitConeOrbits(list OC, string fn)
931"USAGE: storeOrbitConeOrbits(OC, fn); OC: list of lists of cones, fn: string
932PURPOSE: Writes OC to file fn in Singular readable format. Can be importet using <.
933RETURN: nothing
934"
935{
936  int i,j;
937  write(":w "+fn,"list OC;");
938  for (i=1;i<=size(OC);i++){
939    write(":a "+fn,"list OCj;");
940    for (j=1;j<=size(OC[i]);j++){
941      write(":a "+fn,coneToString(OC[i][j]));
942      write(":a "+fn,"OCj[size(OCj)+1]=C;");
943      write(":a "+fn,"kill C;kill F;");
944    }
945    write(":a "+fn,"OC[size(OC)+1]=OCj;");
946    write(":a "+fn,"kill OCj;");
947  }
948}
949
950static proc storeCone(cone C,string fn){
951  intmat H = intmat(inequalities(C));
952  write(":w "+fn,"bigintmat H["+string(nrows(H))+"]["+string(ncols(H))+"] =");
953  write(":w "+fn,string(H)+";");
954  write(":a "+fn,"cone C= coneViaInequalities(H); kill H;");
955}
956
957static proc listToIntvec(list L){
958  intvec v;
959  for (int i = 1;i<=size(L);i++){
960    v[i]=L[i];}
961  return(v);
962}
963
964static proc intvecToList(intvec L){
965  list v;
966  for (int i = 1;i<=size(L);i++){
967    v[i]=L[i];}
968  return(v);
969}
970
971static proc stackMatrices(bigintmat M, bigintmat N){
972  if (ncols(M)!=ncols(N)){ERROR("matrices should have the same number of columns");}
973  bigintmat MN[nrows(M)+nrows(N)][ncols(N)];
974  MN[1..nrows(M),1..ncols(M)]=M[1..nrows(M),1..ncols(M)];
975  MN[(nrows(M)+1)..(nrows(M)+nrows(N)),1..ncols(M)]=N[1..nrows(N),1..ncols(N)];
976  return(MN);
977}
978
979
980proc movingCone(intmat Q)
981"USAGE: movingCone(Q); Q: intmat
982PURPOSE: Computes the moving cone from the grading matrix, with the degrees in the columns of Q.
983RETURN: a cone
984EXAMPLE: example movingCone; shows an example
985"
986{
987  cone Qgammadual = coneViaInequalities(transpose(Q));
988  bigintmat R = facets(Qgammadual);
989  list idx1=intvecToList(intvec(1..nrows(R)));
990  intvec idx2=1..ncols(R);
991  intvec idx1iv;
992  bigintmat Ri[nrows(R)-1][ncols(R)];
993  list C;
994  for (int i = 1; i<=nrows(R);i++){
995    idx1iv = listToIntvec(delete(idx1,i));
996    Ri=R[idx1iv,idx2];
997    C[i]=coneViaPoints(Ri);
998    dbprint(string(i)+"   "+string(nrows(facets(C[i]))));
999  }
1000  cone mov = convexIntersection(C);
1001  return(mov);
1002}
1003example
1004{
1005  echo=2;
1006  intmat Q[5][10] =
1007    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
1008    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
1009    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
1010    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
1011    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
1012  cone mov = movingCone(Q);
1013  mov;
1014  rays(mov);
1015}
1016
1017
1018static proc pivotIndices(matrix H)
1019{
1020  intvec p;
1021  int pp;
1022  int i,j;
1023  int l=1;
1024  for (i=1; i<=ncols(H); i++)
1025  {
1026    for (j=nrows(H); j>=0; j--)
1027    {
1028      if (H[j,i]!=0)
1029      {
1030        if (j>pp)
1031        {
1032          p[l] = i;
1033          l++;
1034          pp = j;
1035        }
1036        break;
1037      }
1038    }
1039  }
1040  return (p);
1041}
1042
1043
1044
1045proc groupActionOnQImage(list G,intmat Q)
1046"USAGE: groupActionOnQImage(G,Q); G: list of permutations, Q: intmat
1047PURPOSE: Given the group G of permutations acting on the simplex on ncols(Q) objects, computes the corresponding group action on the image of Q. We assume that the basering has characteristic 0.
1048RETURN: list of matrices
1049EXAMPLE: example groupActionOnQImage; shows an example
1050"
1051{
1052  matrix Qmat = transpose(matrix(Q));
1053  matrix H = gauss_nf(Qmat);
1054  intvec indices = pivotIndices(H);
1055  intmat Qbasis[nrows(Q)][size(indices)]=Q[1..nrows(Q),indices];
1056  matrix QbasisInv = inverse(Qbasis);
1057  list L;
1058  intmat Qt = transpose(Q);
1059  for(int i = 1; i <= size(G); i++){
1060    intvec sig = permutationToIntvec(G[i]);
1061    intmat Bsig = perm2mat(sig, indices,Qt);
1062    matrix Asig = Bsig * QbasisInv;
1063    L[size(L)+1] = matrixToIntmat(Asig);
1064    kill sig;
1065    kill Bsig;
1066    kill Asig;
1067  }
1068  return(L);
1069}
1070example {
1071  echo=2;
1072  ring R = 0,(x),dp;
1073  intmat Q[5][10] =
1074    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
1075    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
1076    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
1077    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
1078    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
1079list generatorsG = permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )),
1080permutationFromIntvec(intvec( 5, 7, 1, 6, 9, 2, 8, 4, 10, 3 ));
1081groupActionOnQImage(generatorsG,Q);
1082}
1083
1084
1085
1086static proc perm2mat(intvec sig, intvec indices, intmat Q){
1087  intvec sigind = 1:size(indices);
1088  for(int i = 1; i <= size(indices); i++){
1089    if(indices[i] > size(sig)){
1090      sigind[i] = indices[i];
1091    } else {
1092      sigind[i] = sig[indices[i]];
1093    }
1094  }
1095  intmat Asig[size(indices)][ncols(Q)];
1096  for(i = 1; i <= size(sigind); i++){
1097    Asig[i,1..ncols(Q)] = Q[sigind[i], 1..ncols(Q)];
1098  }
1099  // the new basis must be in the cols:
1100  return(transpose(Asig));
1101}
1102
1103///////
1104
1105
1106
1107
1108static proc imageCone(cone c, intmat A){
1109  bigintmat ineqs = inequalities(c);
1110  cone cc = coneViaInequalities (ineqs * A);
1111  return(cc);
1112}
1113
1114
1115static proc matrixToIntmat(matrix A){
1116  int i,j;
1117  intmat Aint[nrows(A)][ncols(A)];
1118  for(i = 1; i<=nrows(A); i++){
1119    for(j = 1; j<=ncols(A); j++){
1120      if (deg(A[i,j])>0){ERROR("entries should be constants");}
1121      if (denominator(number(A[i,j]))!=1){ERROR("entries should be integers");}
1122      Aint[i,j]=int(number(A[i,j]));
1123    }
1124  }
1125  return(Aint);
1126}
1127
1128
1129proc groupActionOnHashes(list Asigma, list OCmov)
1130"USAGE: groupActionOnHashes(Asigma,OCmov); Asigma: list, OCmov: list of list of cones
1131PURPOSE: From the list of orbits of orbitcones, and the symmetry group represenation given by the matrices in Asigma, compute the corresponding permutation representation of the symmetry group on the orbit cones. The permutations are specified in a map representation of length the sum of the size of the orbits of OCmov.
1132RETURN: list of permutations
1133EXAMPLE: example groupActionOnHashes; shows an example
1134"
1135{
1136// for each A in S6:
1137// for each c in OC90:
1138// store index of A*c
1139// --> intvec v_A
1140// store this in the list Ind
1141  list Ind;
1142  int i,j,b,k,sizepreviousorbits;
1143  list remaining;
1144  intmat A;
1145  cone c,Ac;
1146  for(i = 1; i<=size(Asigma); i++){
1147    intvec vA;
1148    A = intmat(Asigma[i]);
1149    dbprint("element "+string(i)+" of symmetry group");
1150    sizepreviousorbits=0;
1151    for(b = 1; b <= size(OCmov); b++){
1152      remaining = intvecToList(intvec(1..size(OCmov[b])));
1153      for(j = 1; j <= size(OCmov[b]); j++){
1154        Ac = imageCone(OCmov[b][j], A);
1155
1156        // find out index:
1157        for(k= 1; k <= size(remaining); k++){
1158          if(OCmov[b][remaining[k]] == Ac){
1159            vA[j+sizepreviousorbits] = remaining[k]+sizepreviousorbits;
1160            remaining = delete(remaining,k);
1161            break;
1162          }
1163        }
1164      }
1165
1166      sizepreviousorbits=size(vA);
1167    }
1168    Ind[size(Ind)+1] = permutationFromIntvec(vA);
1169    dbprint("vA: "+string(vA));
1170    kill vA;
1171  }
1172  return(Ind);
1173}
1174example
1175{
1176  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
1177  ideal J =
1178    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
1179    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
1180    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
1181    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
1182    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
1183  intmat Q[5][10] =
1184    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
1185    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
1186    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
1187    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
1188    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
1189  list AF= afaces(J);
1190  list OC = orbitCones(AF,Q);
1191  list generatorsG = permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )),
1192                     permutationFromIntvec(intvec( 5, 7, 1, 6, 9, 2, 8, 4, 10, 3 ));
1193  list Asigmagens = groupActionOnQImage(generatorsG,Q);
1194groupActionOnHashes(Asigmagens,list(OC));
1195
1196
1197  list simplexSymmetryGroup = G25Action();
1198  list orb = findOrbits(simplexSymmetryGroup,nrows(Q));
1199  list simplexOrbitRepresentatives;
1200  for (int i=1;i<=size(orb);i++){simplexOrbitRepresentatives[i]=orb[i][1];}
1201  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
1202  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
1203  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
1204  list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits);
1205  list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q);
1206  list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits);
1207  list Asigma = groupActionOnQImage(simplexSymmetryGroup,Q);
1208  groupActionOnHashes(Asigma,listOfOrbitConeOrbits);
1209
1210}
1211
1212
1213
1214
1215
1216
1217
1218
1219static proc composePermutationsGAP(permutation sigma, permutation tau){
1220  intvec sigmaTauImage;
1221
1222  for (int i=1;i<=size(sigma.image);i++){
1223    sigmaTauImage[i]=tau.image[sigma.image[i]];
1224  }
1225
1226  permutation sigmaTau;
1227  sigmaTau.image = sigmaTauImage;
1228  return(sigmaTau);
1229}
1230
1231static proc composePermutations(permutation sigma, permutation tau){
1232  intvec sigmaTauImage;
1233
1234  for (int i=1;i<=size(sigma.image);i++){
1235    sigmaTauImage[i]=sigma.image[tau.image[i]];
1236  }
1237
1238  permutation sigmaTau;
1239  sigmaTau.image = sigmaTauImage;
1240  return(sigmaTau);
1241}
1242
1243
1244static proc permutationPower(permutation sigma, int k){
1245  int i;
1246  if (k==0)
1247  {
1248    permutation identity;
1249    identity.image = intvec(1..size(sigma.image));
1250    return (identity);
1251  }
1252  if (k<0)
1253  {
1254    // if k<0 invert sigma and replace k with -k
1255    intvec sigmaImage = sigma.image;
1256    for (i=1; i<=size(sigmaImage); i++)
1257    {
1258      sigma.image[sigmaImage[i]] = i;
1259    }
1260    k = -k;
1261  }
1262  permutation sigmaToPowerK = sigma;
1263  for (i=2; i<=k; i++)
1264  {
1265    sigmaToPowerK = composePermutations(sigmaToPowerK, sigma);
1266  }
1267  return (sigmaToPowerK);
1268}
1269
1270proc permutationFromIntvec(intvec sigmaImage)
1271"USAGE: permutationFromIntvec(sigmaImage); sigmaImage: intvec
1272PURPOSE: Create a permutation from an intvec of images.
1273RETURN: permutation
1274EXAMPLE: example permutationFromIntvec; shows an example
1275"
1276{
1277  permutation sigma;
1278  sigma.image = sigmaImage;
1279  return (sigma);
1280}
1281
1282
1283proc evaluateProduct(list generatorsGperm, string st)
1284"USAGE: evaluateProduct(generatorsGperm,st); generatorsGperm: list, st: string
1285PURPOSE: Evaluates a formal product of variables xi in st, where xi corresponds to the permutation generatorsGperm[i].
1286RETURN: permutation
1287EXAMPLE: example evaluateProduct; shows an example
1288"
1289{
1290  for (int i=1;i<=size(generatorsGperm);i++){
1291    execute("permutation x"+string(i)+"= generatorsGperm[i];");
1292  }
1293  if (st==""){
1294    st="x1^0";
1295  }
1296  execute("permutation sigma = "+st);
1297  return(sigma);}
1298example
1299{
1300  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
1301  ideal J =
1302    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
1303    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
1304    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
1305    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
1306    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
1307  intmat Q[5][10] =
1308    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
1309    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
1310    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
1311    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
1312    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
1313  list AF= afaces(J);
1314  list OC = orbitCones(AF,Q);
1315  list generatorsG = permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )),
1316                     permutationFromIntvec(intvec( 5, 7, 1, 6, 9, 2, 8, 4, 10, 3 ));
1317  list Asigmagens = groupActionOnQImage(generatorsG,Q);
1318//list actionOnOrbitconeIndicesForGenerators = groupActionOnHashes(Asigmagens,OC);
1319string elementInTermsOfGenerators =
1320"(x2^-1*x1^-1)^3*x1^-1";
1321//evaluateProduct(actionOnOrbitconeIndicesForGenerators, elementInTermsOfGenerators);
1322}
1323
1324
1325proc permutationToIntvec(permutation sigma)
1326"USAGE: permutationToIntvec(sigma); sigma: permutation
1327PURPOSE: Convert a permutation to an intvec of images.
1328RETURN: intvec
1329EXAMPLE: example permutationToIntvec; shows an example
1330"
1331{return(sigma.image);}
1332example
1333{
1334permutation sigma = permutationFromIntvec(intvec( 1, 2, 4, 3, 5, 7, 6, 9, 8, 10 ));
1335sigma;
1336permutationToIntvec(sigma);
1337}
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350proc storeActionOnOrbitConeIndices(list Ind,string fn)
1351"USAGE: storeActionOnOrbitConeIndices(generatorsGperm,st); generatorsGperm: list, fn: string
1352PURPOSE: Write the action on the set of orbit cones to the file fn in Singular readable format.
1353RETURN: nothing
1354"
1355{
1356  string s = "list actionOnOrbitconeIndices;";
1357  for(int i =1; i <= size(Ind); i++){
1358    s = s + "intvec v = " + string(Ind[i]) + ";" + newline;
1359    s = s + "actionOnOrbitconeIndices[size(actionOnOrbitconeIndices)+1] = permutationFromIntvec(v);" + newline + "kill v;" + newline;
1360  }
1361  write(":w "+fn, s);
1362}
1363
1364
1365
1366
1367
1368static proc getNeighborHash(list OC, bigintmat w, bigintmat v, int mu)
1369{
1370  int success = 0;
1371  int zz;
1372  intvec Jtmp;
1373  bigintmat wtmp;
1374
1375  while(!success)
1376  {
1377    mu = mu*2;
1378    wtmp = mu*w - v;
1379    success = 1;
1380    for(zz = 1; zz <= size(OC); zz++)
1381    {
1382      if(containsInSupport(OC[zz], wtmp))
1383      {
1384        if(!containsInSupport(OC[zz], w))
1385        {
1386          success = 0;
1387          Jtmp = 0;
1388          break;
1389        }
1390        // insert index zz:
1391        if(size(Jtmp) ==1 && Jtmp[1] == 0)
1392        {
1393          Jtmp[1] = zz;
1394        }
1395        else
1396        {
1397          Jtmp[size(Jtmp)+1] = zz;
1398        }
1399      }
1400    }
1401  }
1402
1403  return(Jtmp);
1404}
1405
1406
1407
1408
1409///////////////////////////////////////
1410
1411proc GITfanSymmetric(list OC, bigintmat Q, cone Qgamma, list actiononorbitcones, list #)
1412  "USAGE: GITfanSymmetric(OC, Q, Qgamma, actiononorbitcones [, file1, file2]); OC:list, Q:bigintmat, Qgamma:cone, actiononorbitcones: list of intvec, file1:string, file2:string
1413PURPOSE: Returns the common refinement of the cones given in
1414the list OC which is supposed to contain the orbit cones intersected with Qgamma. The list actiononorbitcones is supposed to contain the symmetry group acting as permutations of on the list of orbit cones in OC. The optional argument can be used to specify one or two strings with file names, where the first file will contain the hashes of the GIT-cones and the second argument the actual cones in their H-representation.
1415To obtain the whole GIT-fan Qgamma has to be take the cone generated by the columns of Q.
1416RETURN: a list containing the bigint hashes of the GIT cones.
1417EXAMPLE: example GITfanSymmetric; shows an example
1418"
1419{
1420  actiononorbitcones=apply(actiononorbitcones,permutationToIntvec);
1421  /**
1422   * stores the hashes of all maximal GIT cones computed
1423   */
1424  list hashesOfCones;
1425  /**
1426   * stores to the corresponding maximal GIT cone in hashesOfCones
1427   * - 0, if guaranteed that all adjacent GIT cones are known
1428   *        or are to be computed in the next iteration       (*)
1429   * - hash as intvec, otherwise
1430   */
1431  list workingList;
1432
1433
1434  /**
1435   * compute starting cone
1436   */
1437  bigintmat w,v;
1438  cone lambda;
1439  intvec lambdaHash;
1440  while(dimension(lambda) <> nrows(Q)){
1441    w = randConeEl(transpose(Q),100);
1442    dbprint("testing "+string(w));
1443    if (containsRelatively(Qgamma,w)) {
1444      lambda,lambdaHash = GITcone(OC,w);
1445      dbprint("computed cone of dimension "+string(dimension(lambda)));
1446    }
1447  }
1448  int nCones = 1;     // essentially size(hashesOfCones)
1449  int nConesOpen = 1; // essentially size(workingList)+1, see (*) above
1450
1451
1452  /**
1453   * initialize lists
1454   */
1455  int posToInsert = 1;
1456  bigint hashToInsert = binaryToBigint(lambdaHash);
1457  intvec sigLambdaHash;
1458  for(int i = 2; i <= size(actiononorbitcones); i++){
1459    sigLambdaHash = composeIntvecs(actiononorbitcones[i],lambdaHash);
1460    hashToInsert = min(hashToInsert,binaryToBigint(sigLambdaHash));
1461  }
1462  hashesOfCones[1] = hashToInsert;
1463  workingList[1] = int(0);
1464  if (size(#)>0) {write(":w "+#[1],string(hashesOfCones[1]) + ",");}
1465  if (size(#)>1) {write(":w "+#[2],"");writeGitconeToFile(lambda,#[2]);}
1466
1467
1468  /**
1469   * traverse fan
1470   */
1471  int t,tt;
1472  list FL;
1473  intvec neighbourHash;
1474  int mu = 1024;
1475  while (lambdaHash>0)
1476  {
1477    tt=timer;
1478
1479    /**
1480     * compute all facets of lambda
1481     */
1482    t = timer;
1483    FL = listOfFacetsAndInteriorVectors(lambda, Qgamma);
1484    dbprint("time for facets: "+string(timer - t));
1485
1486    /**
1487     * compute all neighbours of lambda
1488     */
1489    for (i=size(FL[1]); i>0; i--)
1490    {
1491      v = FL[1][i][1]; // interior facet normal
1492      w = FL[1][i][2]; // interior facet point
1493      neighbourHash = getNeighborHash(OC,w,v,mu);
1494      posToInsert,hashToInsert = findPosToInsertSymmetric(hashesOfCones,neighbourHash,actiononorbitcones);
1495      if(posToInsert > 0)
1496      {
1497        if (size(#)>0){write(":a "+#[1],string(binaryToBigint(neighbourHash)) + ",");}
1498        hashesOfCones = insertToList(hashesOfCones,hashToInsert,posToInsert);
1499        workingList = insertToList(workingList,neighbourHash,posToInsert);
1500        nConesOpen++;
1501        nCones++;
1502      }
1503    }
1504
1505    /**
1506     * pick lambdaHash and lambda for next iteration,
1507     * set respective entry in workingList to 0
1508     */
1509    lambdaHash = 0;
1510    for (i=size(workingList); i>0; i--)
1511    {
1512      if (typeof(workingList[i])=="intvec")
1513      {
1514        lambdaHash = workingList[i];
1515        lambda = gitConeFromHash(OC,lambdaHash);
1516        if (size(#)>1) {writeGitconeToFile(lambda,#[2]);}
1517        workingList[i] = int(0);
1518        break;
1519      }
1520    }
1521    nConesOpen--;
1522
1523    dbprint("overall: "+string(nCones)+"   open: "+string(nConesOpen)+"   time for loop: "+string(timer-tt));
1524  }
1525  return(hashesOfCones);
1526}
1527example
1528{
1529  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
1530  ideal J =
1531    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
1532    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
1533    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
1534    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
1535    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
1536  intmat Q[5][10] =
1537    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
1538    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
1539    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
1540    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
1541    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
1542  list simplexSymmetryGroup = G25Action();
1543  list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ),
1544  intvec( 1, 2, 3, 5, 6 ),
1545  intvec( 1, 2, 3, 5, 7 ),
1546  intvec( 1, 2, 3, 5, 10 ),
1547  intvec( 1, 2, 3, 7, 9 ),
1548  intvec( 1, 2, 6, 9, 10 ),
1549  intvec( 1, 2, 3, 4, 5, 6 ),
1550  intvec( 1, 2, 3, 4, 5, 10 ),
1551  intvec( 1, 2, 3, 5, 6, 8 ),
1552  intvec( 1, 2, 3, 5, 6, 9 ),
1553  intvec( 1, 2, 3, 5, 7, 10 ),
1554  intvec( 1, 2, 3, 7, 9, 10 ),
1555  intvec( 1, 2, 3, 4, 5, 6, 7 ),
1556  intvec( 1, 2, 3, 4, 5, 6, 8 ),
1557  intvec( 1, 2, 3, 4, 5, 6, 9 ),
1558  intvec( 1, 2, 3, 5, 6, 9, 10 ),
1559  intvec( 1, 2, 3, 4, 5, 6, 7, 8 ),
1560  intvec( 1, 2, 3, 4, 5, 6, 9, 10 ),
1561  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ),
1562  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 );
1563  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
1564  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
1565  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
1566  apply(afaceOrbits,size);
1567  list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits);
1568  apply(minAfaceOrbits,size);
1569  list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q);
1570  apply(listOfOrbitConeOrbits,size);
1571  list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits);
1572  size(listOfMinimalOrbitConeOrbits);
1573  list Asigma = groupActionOnQImage(simplexSymmetryGroup,Q);
1574  list actionOnOrbitconeIndices = groupActionOnHashes(Asigma,listOfOrbitConeOrbits);
1575  list OClist = listOfOrbitConeOrbits[1];
1576  for (int i =2;i<=size(listOfOrbitConeOrbits);i++){
1577    OClist = OClist + listOfOrbitConeOrbits[i];
1578  }
1579cone mov = coneViaPoints(transpose(Q));
1580mov = canonicalizeCone(mov);
1581printlevel = 3;
1582list Sigma = GITfanSymmetric(OClist, Q, mov, actionOnOrbitconeIndices);
1583Sigma;
1584}
1585
1586
1587proc GITfanParallelSymmetric(list OC, bigintmat Q, cone Qgamma, list actiononorbitcones, list #)
1588  "USAGE: GITfanParallelSymmetric(OC, Q, Qgamma, actiononorbitcones [, file1]); OC:list, Q:bigintmat, Qgamma:cone, actiononorbitcones: list of intvec, file1:string
1589PURPOSE: Returns the common refinement of the cones given in
1590the list OC which is supposed to contain the orbit cones intersected with Qgamma. The list actiononorbitcones is supposed to contain the symmetry group acting as permutations of on the list of orbit cones in OC. The optional argument can be used to specify a name for a file which will contain the hashes of the GIT-cones.
1591To obtain the whole GIT-fan Qgamma has to be take the cone generated by the columns of Q.
1592RETURN: a list containing the bigint hashes of the GIT cones.
1593NOTE: The proceduce uses parallel computation for the construction of the GIT-cones.
1594EXAMPLE: example GITfanParallelSymmetric; shows an example
1595"
1596{
1597  actiononorbitcones=apply(actiononorbitcones,permutationToIntvec);
1598  /**
1599   * stores the hashes of all maximal GIT cones computed
1600   */
1601  list hashesOfCones;
1602  /**
1603   * stores to the corresponding maximal GIT cone in hashesOfCones
1604   * - 0, if guaranteed that all adjacent GIT cones are known
1605   *        or are to be computed in the next iteration       (*)
1606   * - hash as intvec, otherwise
1607   */
1608  list workingList;
1609
1610
1611  /**
1612   * compute starting cone
1613   */
1614  bigintmat w,v;
1615  cone lambda;
1616  intvec lambdaHash;
1617  while(dimension(lambda) <> nrows(Q)){
1618    w = randConeEl(transpose(Q),100);
1619    dbprint("testing "+string(w));
1620    if (containsRelatively(Qgamma,w)) {
1621      lambda,lambdaHash = GITcone(OC,w);
1622      dbprint("computed cone of dimension "+string(dimension(lambda)));
1623    }
1624  }
1625  int nCones = 1;     // essentially size(hashesOfCones)
1626  int nConesOpen = 1; // essentially size(workingList)+1, see (*) above
1627
1628
1629  /**
1630   * initialize lists
1631   */
1632  bigint hashToInsert = binaryToBigint(lambdaHash);
1633  int posToInsert = 1;
1634  intvec sigLambdaHash;
1635  for(int i = 2; i <= size(actiononorbitcones); i++){
1636    sigLambdaHash = composeIntvecs(actiononorbitcones[i],lambdaHash);
1637    hashToInsert = min(hashToInsert,binaryToBigint(sigLambdaHash));
1638  }
1639  hashesOfCones[1] = hashToInsert;
1640  workingList[1] = int(0);
1641  if (size(#)>0) {write(":w "+#[1],string(binaryToBigint(lambdaHash)) + ",");}
1642  //if (size(#)>1) {write(":w "+#[2],"list listOfMaximalCones;");writeGitconeToFile(lambda,#[2]);}
1643
1644
1645  list iterationArgs = list(list(lambdaHash,OC,Qgamma,actiononorbitcones));
1646  list iterationRes;
1647
1648  /**
1649   * traverse fan
1650   */
1651  int j,t,tloop;
1652  list FL;
1653  list neighbourHashes;
1654  intvec neighbourHash;
1655  while (size(iterationArgs)>0)
1656  {
1657    tloop=rtimer;
1658
1659    /**
1660     * compute all neighbours of lambda
1661     */
1662    t = rtimer;
1663    iterationRes = parallelWaitAll("computeNeighbourMinimalHashes",iterationArgs);
1664    dbprint("time neighbours: "+string(rtimer - t));
1665
1666    /**
1667     * central book keeping
1668     */
1669    t = rtimer;
1670    for (i=1; i<=size(iterationRes); i++)
1671    {
1672      neighbourHashes = iterationRes[i];
1673      for (j=1; j<=size(neighbourHashes); j++)
1674      {
1675        neighbourHash = neighbourHashes[j];
1676        hashToInsert = binaryToBigint(neighbourHash);
1677        posToInsert = findPlaceToInsert(hashesOfCones,hashToInsert);
1678        if(posToInsert > 0)
1679        {
1680          if (size(#)>0){write(":a "+#[1],string(binaryToBigint(neighbourHash)) + ",");}
1681          hashesOfCones = insertToList(hashesOfCones,hashToInsert,posToInsert);
1682          workingList = insertToList(workingList,neighbourHash,posToInsert);
1683          nConesOpen++;
1684          nCones++;
1685        }
1686      }
1687    }
1688    nConesOpen = nConesOpen - size(iterationArgs);
1689    dbprint("time bookkeeping: "+string(rtimer - t));
1690
1691    /**
1692     * pick arguments for next iteration,
1693     * set respective entry in workingList to 0
1694     */
1695    t = rtimer;
1696    iterationArgs = list();
1697    for (i=size(workingList); i>0; i--)
1698    {
1699      if (typeof(workingList[i])=="intvec")
1700      {
1701        iterationArgs[size(iterationArgs)+1] = list(workingList[i],OC,Qgamma,actiononorbitcones);
1702        workingList[i] = int(0);
1703        if (size(iterationArgs) >= getcores())
1704        {
1705          break;
1706        }
1707      }
1708    }
1709    dbprint("time preparation: "+string(rtimer - t));
1710
1711    dbprint("overall: "+string(nCones)+"   open: "+string(nConesOpen)+"   time loop: "+string(rtimer-tloop));
1712  }
1713  return(hashesOfCones);
1714}
1715example
1716{
1717  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
1718  ideal J =
1719    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
1720    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
1721    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
1722    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
1723    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
1724  intmat Q[5][10] =
1725    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
1726    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
1727    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
1728    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
1729    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
1730  list simplexSymmetryGroup = G25Action();
1731  list simplexOrbitRepresentatives = intvec( 1, 2, 3, 4, 5 ),
1732  intvec( 1, 2, 3, 5, 6 ),
1733  intvec( 1, 2, 3, 5, 7 ),
1734  intvec( 1, 2, 3, 5, 10 ),
1735  intvec( 1, 2, 3, 7, 9 ),
1736  intvec( 1, 2, 6, 9, 10 ),
1737  intvec( 1, 2, 3, 4, 5, 6 ),
1738  intvec( 1, 2, 3, 4, 5, 10 ),
1739  intvec( 1, 2, 3, 5, 6, 8 ),
1740  intvec( 1, 2, 3, 5, 6, 9 ),
1741  intvec( 1, 2, 3, 5, 7, 10 ),
1742  intvec( 1, 2, 3, 7, 9, 10 ),
1743  intvec( 1, 2, 3, 4, 5, 6, 7 ),
1744  intvec( 1, 2, 3, 4, 5, 6, 8 ),
1745  intvec( 1, 2, 3, 4, 5, 6, 9 ),
1746  intvec( 1, 2, 3, 5, 6, 9, 10 ),
1747  intvec( 1, 2, 3, 4, 5, 6, 7, 8 ),
1748  intvec( 1, 2, 3, 4, 5, 6, 9, 10 ),
1749  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9 ),
1750  intvec( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 );
1751  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
1752  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
1753  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
1754  apply(afaceOrbits,size);
1755  list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits);
1756  apply(minAfaceOrbits,size);
1757  list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q);
1758  apply(listOfOrbitConeOrbits,size);
1759  list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits);
1760  size(listOfMinimalOrbitConeOrbits);
1761  list Asigma = groupActionOnQImage(simplexSymmetryGroup,Q);
1762  list actionOnOrbitconeIndices = groupActionOnHashes(Asigma,listOfOrbitConeOrbits);
1763  list OClist = listOfOrbitConeOrbits[1];
1764  for (int i =2;i<=size(listOfOrbitConeOrbits);i++){
1765    OClist = OClist + listOfOrbitConeOrbits[i];
1766  }
1767  cone mov = coneViaPoints(transpose(Q));
1768  mov = canonicalizeCone(mov);
1769  list Sigma = GITfanParallelSymmetric(OClist, Q, mov, actionOnOrbitconeIndices);
1770  Sigma;
1771}
1772
1773
1774// not static to be used in parallel.lib
1775proc computeNeighbourMinimalHashes(intvec lambdaHash, list OC, cone Qgamma, list actiononorbitcones)
1776{
1777  /**
1778   * compute all facets of lambda
1779   */
1780  cone lambda = gitConeFromHash(OC,lambdaHash);
1781  list FL = listOfFacetsAndInteriorVectors(lambda, Qgamma);
1782
1783  /**
1784   * compute all minimal hashes of neighbours of lambda
1785   */
1786  int i,j;
1787  bigintmat v,w;
1788  intvec neighbourHash;
1789  intvec neighbourHashPerm;
1790  bigint nPerm;
1791  intvec neighbourHashMin;
1792  bigint nMin;
1793  list neighbourHashes;
1794  for (i=size(FL[1]); i>0; i--)
1795  {
1796    v = FL[1][i][1]; // interior facet normal
1797    w = FL[1][i][2]; // interior facet point
1798    neighbourHash = getNeighborHash(OC,w,v,1024);
1799    neighbourHashMin = neighbourHash;
1800    nMin = binaryToBigint(neighbourHash);
1801    for (j=size(actiononorbitcones); j>1; j--)
1802    {
1803      neighbourHashPerm = composeIntvecs(actiononorbitcones[j],neighbourHash);
1804      nPerm = binaryToBigint(neighbourHashPerm);
1805      if (nPerm < nMin)
1806      {
1807        nMin = nPerm;
1808        neighbourHashMin = neighbourHashPerm;
1809      }
1810    }
1811    neighbourHashes[i] = neighbourHashMin;
1812  }
1813
1814  return (neighbourHashes);
1815}
1816
1817
1818
1819static proc writeGitconeToFile(cone lambda,string fn)
1820{
1821  bigintmat H = facets(lambda);
1822
1823  int rows = nrows(H);
1824  int cols = ncols(H);
1825  string toBeWritten = "bigintmat H["+string(rows)+"]["+string(cols)+"]=";
1826  int i,j;
1827  for (i=1; i<=rows; i++)
1828  {
1829    for (j=1; j<=cols; j++)
1830    {
1831      toBeWritten = toBeWritten + string(H[i,j]) + ",";
1832    }
1833  }
1834  toBeWritten = toBeWritten[1..size(toBeWritten)-1];
1835  toBeWritten = toBeWritten + ";";
1836  write(":a "+fn,toBeWritten);
1837
1838  toBeWritten = "listOfMaximalCones[size(listOfMaximalCones)+1] = coneViaInequalities(V);";
1839  write(":a "+fn,toBeWritten);
1840
1841  toBeWritten = "kill V;";
1842  write(":a "+fn,toBeWritten);
1843
1844  toBeWritten = "";
1845  write(":a "+fn,toBeWritten);
1846}
1847
1848
1849
1850static proc listOfFacetsAndInteriorVectors(cone lambda, cone Qgamma, list #){
1851  list FL;
1852  int numboundary;
1853
1854  bigintmat FL0 = facets(lambda);
1855  bigintmat H[1][ncols(FL0)];
1856  bigintmat FL1[nrows(FL0)-1][ncols(FL0)]; // delete row of H from FL0
1857
1858  bigintmat w;
1859  cone facetCone;
1860  for(int i = 1; i <= nrows(FL0); i++){
1861    H = FL0[i,1..ncols(FL0)];
1862
1863    if(i > 1 and i < nrows(FL0)){
1864      FL1 = FL0[1..i-1, 1..ncols(FL0)], FL0[i+1..nrows(FL0), 1..ncols(FL0)];
1865    } else {
1866      if(i == nrows(FL0)){
1867        FL1 = FL0[1..i-1, 1..ncols(FL0)];
1868      } else { // i = 1:
1869        FL1 = FL0[2..nrows(FL0), 1..ncols(FL0)];
1870      }
1871    }
1872
1873    facetCone = coneViaInequalities(FL1, H);
1874
1875    if (size(#)>0)
1876    {
1877      if (containsInSupport(facetCone,#[1]))
1878      {
1879        i++;
1880        continue;
1881      }
1882    }
1883
1884    w = relativeInteriorPoint(facetCone);
1885
1886    if(containsRelatively(Qgamma,w)){
1887      FL[size(FL) + 1] = list(H,w);
1888    } else {
1889      numboundary=numboundary+1;
1890    }
1891  }
1892
1893  return(list(FL,numboundary));
1894}
1895
1896
1897
1898
1899
1900////////////////////////////////
1901
1902// CAREFUL: we assume that actiononorbitcones[1] = id.
1903static proc findPosToInsertSymmetric(list hashesOfCones, intvec J, list actiononorbitcones){
1904
1905  // 1.: compute minimal hash sigJ
1906  bigint n = binaryToBigint(J);
1907  intvec sigJ;
1908  for(int i = 2; i <= size(actiononorbitcones); i++){
1909    sigJ = composeIntvecs(actiononorbitcones[i],J);
1910    n = min(n,binaryToBigint(sigJ));
1911  }
1912
1913  // 2.:check if minimal hash is already in list
1914  return(findPlaceToInsert(hashesOfCones,n),n);
1915}
1916
1917//////////////////////
1918
1919// insert n at position pos and move bigger elements by one index
1920proc insertToList(list L, def elementToInsert, int posToInsert){
1921  if(posToInsert == 1){
1922    return(list(elementToInsert)+L);
1923  }
1924  if(posToInsert == size(L)+1){
1925    return(L+list(elementToInsert));
1926  }
1927  return(list(L[1..(posToInsert-1)],elementToInsert,L[posToInsert..size(L)]));
1928}
1929
1930proc GITcone(list OCcones, bigintmat w)
1931"USAGE: GITcone(OCcones, w); OCcones: list of orbit cones, w: bigintmat with one row
1932PURPOSE: Returns the intersection of all orbit cones containing w.
1933RETURN: cone,intvec with the GIT cone containing w, and the hash of this cone (the indices of the orbit cones contributing to the intersection)
1934EXAMPLE: example GITcone; shows an example
1935"
1936{
1937  list OCrelcones;
1938  intvec Hash;
1939  for(int i = 1; i <= size(OCcones); i++){
1940    if(containsInSupport(OCcones[i], w)){
1941      OCrelcones[size(OCrelcones)+1]=OCcones[i];
1942      Hash[size(Hash)+1] = i;
1943    }
1944  }
1945  Hash = intvec(Hash[2..size(Hash)]);
1946  return(convexIntersection(OCrelcones),Hash);
1947}
1948example {
1949  echo=2;
1950  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
1951  ideal J =
1952    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
1953    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
1954    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
1955    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
1956    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
1957  intmat Q[5][10] =
1958    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
1959    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
1960    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
1961    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
1962    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
1963
1964  list AF= afaces(J,nrows(Q));
1965  AF=fullDimImages(AF,Q);
1966  AF = minimalAfaces(AF);
1967  list OC = orbitCones(AF,Q);
1968  bigintmat w[1][nrows(Q)];
1969  int j;
1970  for(int i = 1; i <= nrows(Q); i++)
1971  {
1972    for(j=1;j<=nrows(Q);j++)
1973    {
1974      w[1,j]= w[1,j] + Q[1,i];
1975    }
1976  }
1977  GITcone(OC,w);
1978}
1979
1980
1981static proc gitConeFromHash(list OC, intvec Hash)
1982{
1983  list OCtoIntersect = OC[Hash];
1984  return (convexIntersection(OCtoIntersect));
1985}
1986
1987
1988
1989
1990
1991static proc randConeEl(bigintmat Q, int bound){
1992  bigintmat w[1][ncols(Q)];
1993
1994  for(int i = 1; i <= nrows(Q); i++){
1995    bigintmat v[1][ncols(Q)] = Q[i,1..ncols(Q)];
1996
1997    w = w + random(1,bound) * v;
1998
1999    kill v;
2000  }
2001
2002  return(w);
2003}
2004
2005
2006
2007static proc subdividelist(list OC,int ncores){
2008  if (ncores>size(OC)){ncores=size(OC);}
2009  int percore = size(OC) div (ncores);
2010  int lastcore = size(OC) mod (ncores);
2011  list OCLL;
2012  int j=1;
2013  int starti=1;
2014  int endi;
2015  for (int i=1;i<=ncores;i++){
2016    endi = starti+percore-1;
2017    if (i<=lastcore) {endi=endi+1;}
2018    list OCLLi=OC[starti..endi];
2019    starti=endi+1;
2020    OCLL[i]=OCLLi;
2021    kill OCLLi;
2022  }
2023  return(OCLL);
2024}
2025//list OC = 1,2,3,4,5,6,7,8,9,10,11,12;
2026//subdividelist(OC,4);
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039proc orbitCones(list AF, intmat Q, list #)
2040"USAGE: orbitCones(AF, Q[, d]); AF: list of intvecs, Q: intmat, d: int
2041PURPOSE: Returns the list consisting of all cones Q(gam0) where gam0 in AF. If the optional argument d is given then the function returns only the orbit cones of dimension at least d
2042RETURN: a list of cones
2043EXAMPLE: example orbitCones; shows an example
2044"
2045{
2046  list OC;
2047  intvec gam0;
2048  int j;
2049  cone c;
2050  int d = 0;
2051  if (size(#)>0) {
2052    d = #[1];
2053  }
2054
2055  for(int i = 1; i <= size(AF); i++){
2056    gam0 = AF[i];
2057
2058    if(gam0 == 0){
2059      bigintmat M[1][nrows(Q)];
2060    } else {
2061      bigintmat M[size(gam0)][nrows(Q)];
2062
2063      for (j = 1; j <= size(gam0); j++){
2064        M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]];
2065      }
2066    }
2067
2068    c = coneViaPoints(M);
2069
2070    if(listContainsCone(OC, c) == 0 && dimension(c)>=d){
2071      OC[size(OC)+1] = c;
2072    }
2073
2074    kill M;
2075  }
2076
2077  return(OC);
2078}
2079example
2080{
2081  echo=2;
2082  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
2083  ideal J =
2084    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
2085    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
2086    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
2087    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
2088    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
2089  intmat Q[5][10] =
2090    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
2091    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
2092    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
2093    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
2094    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
2095
2096  list AF= afaces(J);
2097  print(size(AF));
2098  list OC = orbitCones(AF,Q);
2099  size(OC);
2100}
2101
2102
2103
2104///////////////////////////////////////
2105
2106proc GITfanFromOrbitCones(list OC, bigintmat Q, cone Qgamma, list #)
2107  "USAGE: GITfanFromOrbitCones(OC, Q, Qgamma [, file1, file2]); OC:list, Q:bigintmat, Qgamma:cone, file1:string, file2:string
2108PURPOSE: Returns the common refinement of the cones given in
2109the list OC which is supposed to contain the orbit cones intersected with Qgamma. The optional argument can be used
2110to specify one or two strings with file names, where the first file will contain the hashes of the
2111GIT-cones and the second argument the actual cones in their H-representation.
2112To obtain the whole GIT-fan Qgamma has to be take the cone generated by the columns of Q.
2113RETURN: a list containing the bigint hashes of the GIT cones.
2114EXAMPLE: example GITfanFromOrbitCones; shows an example
2115"
2116{
2117  /**
2118   * stores the hashes of all maximal GIT cones computed
2119   */
2120  list hashesOfCones;
2121  /**
2122   * stores to the corresponding maximal GIT cone in hashesOfCones
2123   * - 0, if guaranteed that all adjacent GIT cones are known
2124   *        or are to be computed in the next iteration       (*)
2125   * - hash as intvec, otherwise
2126   */
2127  list workingList;
2128
2129
2130  /**
2131   * compute starting cone
2132   */
2133  bigintmat w,v;
2134  cone lambda;
2135  intvec lambdaHash;
2136  while(dimension(lambda) <> nrows(Q)){
2137    w = randConeEl(transpose(Q),100);
2138    dbprint("testing "+string(w));
2139    if (containsRelatively(Qgamma,w)) {
2140      lambda,lambdaHash = GITcone(OC,w);
2141      dbprint("computed cone of dimension "+string(dimension(lambda)));
2142    }
2143  }
2144  int nCones = 1;     // essentially size(hashesOfCones)
2145  int nConesOpen = 1; // essentially size(workingList)+1, see (*) above
2146
2147
2148  /**
2149   * initialize lists
2150   */
2151  int posToInsert = 1;
2152  bigint hashToInsert = binaryToBigint(lambdaHash);
2153  hashesOfCones[1] = hashToInsert;
2154  workingList[1] = int(0);
2155  if (size(#)>0) {write(":w "+#[1],string(hashesOfCones[1]) + ",");}
2156  if (size(#)>1) {write(":w "+#[2],"list listOfMaximalCones;");writeGitconeToFile(lambda,#[2]);}
2157
2158
2159  /**
2160   * traverse fan
2161   */
2162  int i,t,tt;
2163  list FL;
2164  intvec neighbourHash;
2165  int mu = 1024;
2166  while (lambdaHash>0)
2167  {
2168    tt=timer;
2169
2170    /**
2171     * compute all facets of lambda
2172     */
2173    t = timer;
2174    FL = listOfFacetsAndInteriorVectors(lambda, Qgamma);
2175    dbprint("time for facets: "+string(timer - t));
2176
2177    /**
2178     * compute all neighbours of lambda
2179     */
2180    for (i=size(FL[1]); i>0; i--)
2181    {
2182      v = FL[1][i][1]; // interior facet normal
2183      w = FL[1][i][2]; // interior facet point
2184      neighbourHash = getNeighborHash(OC,w,v,mu);
2185      hashToInsert = binaryToBigint(neighbourHash);
2186      posToInsert = findPlaceToInsert(hashesOfCones,hashToInsert);
2187
2188      if(posToInsert > 0)
2189      {
2190        if (size(#)>0){write(":a "+#[1],string(binaryToBigint(neighbourHash)) + ",");}
2191        hashesOfCones = insertToList(hashesOfCones,hashToInsert,posToInsert);
2192        workingList = insertToList(workingList,neighbourHash,posToInsert);
2193        nConesOpen++;
2194        nCones++;
2195      }
2196    }
2197
2198    /**
2199     * pick lambdaHash and lambda for next iteration,
2200     * set respective entry in workingList to 0
2201     */
2202    lambdaHash = 0;
2203    for (i=size(workingList); i>0; i--)
2204    {
2205      if (typeof(workingList[i])=="intvec")
2206      {
2207        lambdaHash = workingList[i];
2208        lambda = gitConeFromHash(OC,lambdaHash);
2209        if (size(#)>1) {writeGitconeToFile(lambda,#[2]);}
2210        workingList[i] = int(0);
2211        break;
2212      }
2213    }
2214    nConesOpen--;
2215
2216    dbprint("overall: "+string(nCones)+"   open: "+string(nConesOpen)+"   time for loop: "+string(timer-tt));
2217  }
2218  return(hashesOfCones);
2219}
2220example
2221{
2222  echo=2;
2223  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
2224  ideal J =
2225    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
2226    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
2227    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
2228    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
2229    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
2230  intmat Q[5][10] =
2231    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
2232    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
2233    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
2234    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
2235    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
2236
2237  list AF= afaces(J,nrows(Q));
2238  AF=fullDimImages(AF,Q);
2239  AF = minimalAfaces(AF);
2240  list OC = orbitCones(AF,Q);
2241  cone Qgamma = coneViaPoints(transpose(Q));
2242  list GIT = GITfanFromOrbitCones(OC,Q,Qgamma);
2243  size(GIT);
2244}
2245
2246
2247
2248
2249
2250proc GITfanParallel(list OC, intmat Q, cone Qgamma, list #)
2251  "USAGE: GITfanParallel(OC, Q, Qgamma [, file1]); OC:list, Q:intmat, Qgamma:cone, file1:string
2252PURPOSE: Returns the common refinement of the cones given in
2253the list OC which is supposed to contain the orbit cones intersected with Qgamma. The optional argument can be used to specify a name for a file which will contain the hashes of the GIT-cones.
2254To obtain the whole GIT-fan Qgamma has to be take the cone generated by the columns of Q.
2255RETURN: a list containing the bigint hashes of the GIT cones.
2256NOTE: The proceduce uses parallel computation for the construction of the GIT-cones.
2257EXAMPLE: example GITfanParallel; shows an example
2258"
2259{
2260  /**
2261   * stores the hashes of all maximal GIT cones computed
2262   */
2263  list hashesOfCones;
2264  /**
2265   * stores to the corresponding maximal GIT cone in hashesOfCones
2266   * - 0, if guaranteed that all adjacent GIT cones are known
2267   *        or are to be computed in the next iteration       (*)
2268   * - hash as intvec, otherwise
2269   */
2270  list workingList;
2271
2272
2273  /**
2274   * compute starting cone
2275   */
2276  bigintmat w,v;
2277  cone lambda;
2278  intvec lambdaHash;
2279  while(dimension(lambda) <> nrows(Q)){
2280    w = randConeEl(transpose(Q),100);
2281    dbprint("testing "+string(w));
2282    if (containsRelatively(Qgamma,w)) {
2283      lambda,lambdaHash = GITcone(OC,w);
2284      dbprint("computed cone of dimension "+string(dimension(lambda)));
2285    }
2286  }
2287  int nCones = 1;     // essentially size(hashesOfCones)
2288  int nConesOpen = 1; // essentially size(workingList)+1, see (*) above
2289
2290
2291  /**
2292   * initialize lists
2293   */
2294  bigint hashToInsert = binaryToBigint(lambdaHash);
2295  int posToInsert = 1;
2296  hashesOfCones[1] = hashToInsert;
2297  workingList[1] = int(0);
2298  if (size(#)>0) {write(":w "+#[1],string(binaryToBigint(lambdaHash)) + ",");}
2299  //if (size(#)>1) {write(":w "+#[2],"list listOfMaximalCones;");writeGitconeToFile(lambda,#[2]);}
2300
2301
2302  list iterationArgs = list(list(lambdaHash,OC,Qgamma));
2303  list iterationRes;
2304
2305  /**
2306   * traverse fan
2307   */
2308  int i,j,t,tloop;
2309  list FL;
2310  list neighbourHashes;
2311  intvec neighbourHash;
2312  while (size(iterationArgs)>0)
2313  {
2314    tloop=rtimer;
2315
2316    /**
2317     * compute all neighbours of lambda
2318     */
2319    t = rtimer;
2320    iterationRes = parallelWaitAll("computeNeighbourHashes",iterationArgs);
2321    dbprint("time neighbours: "+string(rtimer - t));
2322
2323    /**
2324     * central book keeping
2325     */
2326    t = rtimer;
2327    for (i=1; i<=size(iterationRes); i++)
2328    {
2329      neighbourHashes = iterationRes[i];
2330      for (j=1; j<=size(neighbourHashes); j++)
2331      {
2332        neighbourHash = neighbourHashes[j];
2333        hashToInsert = binaryToBigint(neighbourHash);
2334        posToInsert = findPlaceToInsert(hashesOfCones,hashToInsert);
2335        if(posToInsert > 0)
2336        {
2337          if (size(#)>0){write(":a "+#[1],string(binaryToBigint(neighbourHash)) + ",");}
2338          hashesOfCones = insertToList(hashesOfCones,hashToInsert,posToInsert);
2339          workingList = insertToList(workingList,neighbourHash,posToInsert);
2340          nConesOpen++;
2341          nCones++;
2342        }
2343      }
2344    }
2345    nConesOpen = nConesOpen - size(iterationArgs);
2346    dbprint("time bookkeeping: "+string(rtimer - t));
2347
2348    /**
2349     * pick arguments for next iteration,
2350     * set respective entry in workingList to 0
2351     */
2352    t = rtimer;
2353    iterationArgs = list();
2354    for (i=size(workingList); i>0; i--)
2355    {
2356      if (typeof(workingList[i])=="intvec")
2357      {
2358        iterationArgs[size(iterationArgs)+1] = list(workingList[i],OC,Qgamma);
2359        workingList[i] = int(0);
2360        if (size(iterationArgs) >= getcores())
2361        {
2362          break;
2363        }
2364      }
2365    }
2366    dbprint("time preparation: "+string(rtimer - t));
2367
2368    dbprint("overall: "+string(nCones)+"   open: "+string(nConesOpen)+"   time loop: "+string(rtimer-tloop));
2369  }
2370
2371  if (size(#)==0)
2372  {
2373    return(hashesOfCones);
2374  }
2375}
2376example
2377{
2378  echo=2;
2379  setcores(4);
2380  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
2381  ideal J =
2382    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
2383    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
2384    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
2385    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
2386    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
2387  intmat Q[5][10] =
2388    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
2389    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
2390    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
2391    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
2392    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
2393
2394  list AF= afaces(J);
2395  print(size(AF));
2396  list OC = orbitCones(AF,Q);
2397  cone Qgamma = coneViaPoints(transpose(Q));
2398  list GIT = GITfanParallel(OC,Q,Qgamma);
2399  size(GIT);
2400}
2401
2402// not static to be used in parallel.lib
2403proc computeNeighbourHashes(intvec lambdaHash, list OC, cone Qgamma)
2404{
2405  /**
2406   * compute all facets of lambda
2407   */
2408  cone lambda = gitConeFromHash(OC,lambdaHash);
2409  list FL = listOfFacetsAndInteriorVectors(lambda, Qgamma);
2410
2411  /**
2412   * compute all minimal hashes of neighbours of lambda
2413   */
2414  bigintmat v,w;
2415  intvec neighbourHash;
2416  list neighbourHashes;
2417  for (int i=size(FL[1]); i>0; i--)
2418  {
2419    v = FL[1][i][1]; // interior facet normal
2420    w = FL[1][i][2]; // interior facet point
2421    neighbourHash = getNeighborHash(OC,w,v,1024);
2422    neighbourHashes[i] = neighbourHash;
2423  }
2424
2425  return (neighbourHashes);
2426}
2427
2428
2429
2430proc minimalAfaces(list listOfAfaces)
2431"USAGE: minimalAfaces(listOfAfaces); listOfAfaces: list
2432PURPOSE: Returns a list of all minimal a-faces. Note that listOfAfaces must only contain
2433afaces which project to full dimension.
2434RETURN: a list of intvecs
2435EXAMPLE: example minimalAfaces; shows an example
2436"
2437{
2438  int i,j;
2439  for (i=1; i<=size(listOfAfaces); i++)
2440  {
2441    for (j=1; j<=size(listOfAfaces); j++)
2442    {
2443      if (i!=j)
2444      {
2445        if (isSubset(listOfAfaces[i],listOfAfaces[j]))
2446        {
2447          listOfAfaces = delete(listOfAfaces,j);
2448          if (j<i)
2449          {
2450            i = i-1;
2451          }
2452          j = j-1;
2453        }
2454      }
2455    }
2456  }
2457  return(listOfAfaces);
2458}
2459example
2460{
2461  echo=2;
2462  setcores(4);
2463  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
2464  ideal J =
2465    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
2466    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
2467    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
2468    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
2469    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
2470  intmat Q[5][10] =
2471    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
2472    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
2473    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
2474    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
2475    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
2476  list AF= afaces(J,nrows(Q));
2477  size(AF);
2478  AF=fullDimImages(AF,Q);
2479  size(AF);
2480  AF=minimalAfaces(AF);
2481  size(AF);
2482}
2483
2484
2485
2486proc GITfan(ideal J, intmat Q, list #)
2487"USAGE: GITfan(J,Q [, G]); J:ideal, Q:intmat, G:list
2488PURPOSE: Computes the GIT fan associated to J and Q. Optionally a symmetry group action on the column space of Q can be specified.
2489RETURN: a fan, the GIT fan.
2490NOTE: The proceduce uses parallel computation for the construction of the GIT-cones. The a-faces are not computed in parallel. This can be done by calling the aface procedure specifying a list of simplex faces. If used with the optional argument G, the orbit decomposition of the simplex of columns of Q is computed. Refer to the Singular documentation on how to do this more efficiently using GAP.
2491EXAMPLE: example GITfan; shows an example
2492"
2493{
2494  list GIT;
2495  list OC;
2496  if (size(#)>0){
2497      (GIT,OC) = GITfanWrapperWithSymmetry(J,Q,#);
2498  } else {
2499  list AF= afaces(J,nrows(Q));
2500  AF=fullDimImages(AF,Q);
2501  AF = minimalAfaces(AF);
2502  OC = orbitCones(AF,Q);
2503  cone Qgamma = coneViaPoints(transpose(Q));
2504  GIT = GITfanParallel(OC,Q,Qgamma);
2505  }
2506  fan Sigma = hashesToFan(GIT,OC);
2507  return(Sigma);
2508}
2509example
2510{
2511  echo=2;
2512  setcores(4);
2513  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
2514  ideal J =
2515    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
2516    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
2517    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
2518    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
2519    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
2520  intmat Q[5][10] =
2521    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
2522    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
2523    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
2524    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
2525    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
2526  fan GIT = GITfan(J,Q);
2527
2528  intmat Q[5][10] =
2529    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
2530    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
2531    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
2532    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
2533    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
2534  list simplexSymmetryGroup = G25Action();
2535  fan GIT2 = GITfan(J,Q,simplexSymmetryGroup);
2536
2537}
2538
2539
2540proc hashToCone(bigint v, list OC)
2541"USAGE: hashToCone(v, OC): v bigint, OC list of cones.
2542ASSUME: the elements of OC are the orbit cones used in the hash representation of the GIT cones.
2543RETURN: a cone, the intersection of the cones in OC according to the binary representation of the hash v.
2544EXAMPLE: example hashToCone; shows an example
2545"
2546{
2547  intvec J = bigintToBinary(v, size(OC));
2548  return(convexIntersection(list(OC[J])));
2549}
2550example
2551{
2552  echo = 2;
2553  setcores(4);
2554  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
2555  ideal J =
2556    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
2557    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
2558    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
2559    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
2560    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
2561  intmat Q[5][10] =
2562    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
2563    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
2564    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
2565    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
2566    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
2567  list AF= afaces(J,nrows(Q));
2568  AF=fullDimImages(AF,Q);
2569  AF = minimalAfaces(AF);
2570  list OC = orbitCones(AF,Q);
2571  bigint v = 21300544;
2572  hashToCone(v, OC);
2573}
2574
2575
2576proc bigintToBinary(bigint n, int r)
2577"
2578USAGE: bigintToBinary(n, r): n bigint, r int.
2579ASSUME: n is smaller then 2^r.
2580RETURN: an intvec, with entries the positions of 1 in the binary representation of n with r bits.
2581EXAMPLE: example bigintToBinary; shows an example
2582"
2583{
2584  int k = r-1;
2585
2586  intvec v;
2587  bigint n0 = n;
2588
2589  while(n0 > 0){
2590    bigint tmp = bigint(2)^k;
2591
2592    while(tmp > n0){
2593      k--;
2594      tmp = bigint(2)^k;
2595    }
2596
2597    v = k+1,v;
2598    n0 = n0 - tmp;
2599    k--;
2600
2601    kill tmp;
2602  }
2603
2604  v = v[1..size(v)-1];
2605  return(v);
2606}
2607example
2608{
2609  echo = 2;
2610  bigintToBinary(bigint(2)^90-1, 90);
2611}
2612
2613
2614
2615
2616proc hashesToFan(list hashes, list OC)
2617"USAGE: hashesToFan(hashes, OC): hashes list of bigint, OC list of cones.
2618ASSUME: the elements of OC are the orbit cones used in the hash representation of the GIT cones.
2619RETURN: a fan, with maximal cones the intersections of the cones in OC according to the binary representation of the hashes.
2620EXAMPLE: example hashesToFan; shows an example
2621"
2622{
2623  fan Sigma = emptyFan(ambientDimension(OC[1]));
2624  for (int i=1;i<=size(hashes);i++){
2625    insertCone(Sigma,hashToCone(hashes[i],OC),0);
2626  }
2627  return(Sigma);
2628}
2629example
2630{
2631  echo = 2;
2632  setcores(4);
2633  ring R = 0,T(1..10),wp(1,1,1,1,1,1,1,1,1,1);
2634  ideal J =
2635    T(5)*T(10)-T(6)*T(9)+T(7)*T(8),
2636    T(1)*T(9)-T(2)*T(7)+T(4)*T(5),
2637    T(1)*T(8)-T(2)*T(6)+T(3)*T(5),
2638    T(1)*T(10)-T(3)*T(7)+T(4)*T(6),
2639    T(2)*T(10)-T(3)*T(9)+T(4)*T(8);
2640  intmat Q[5][10] =
2641    1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
2642    1, 0, 0, 0, 1, 1, 1, 0, 0, 0,
2643    0, 1, 1, 0, 0, 0, -1, 1, 0, 0,
2644    0, 1, 0, 1, 0, -1, 0, 0, 1, 0,
2645    0, 0, 1, 1, -1, 0, 0, 0, 0, 1;
2646  list AF= afaces(J,nrows(Q));
2647  AF=fullDimImages(AF,Q);
2648  AF = minimalAfaces(AF);
2649  list OC = orbitCones(AF,Q);
2650  cone Qgamma = coneViaPoints(transpose(Q));
2651  list GIT = GITfanParallel(OC,Q,Qgamma);
2652  fan Sigma = hashesToFan(GIT,OC);
2653}
2654
2655
2656
2657
2658/////////////////////////////////////
2659
2660proc gkzFan(intmat Q)
2661"USAGE: gkzFan(Q); a: ideal, Q:intmat
2662PURPOSE: Returns the GKZ-fan of the matrix Q.
2663RETURN: a fan.
2664EXAMPLE: example gkzFan; shows an example
2665"
2666{
2667  // only difference to gitFan:
2668  // it suffices to consider all faces
2669  // that are simplicial:
2670  list OC = simplicialToricOrbitCones(Q);
2671  print(size(OC));
2672  cone Qgamma = coneViaPoints(transpose(Q));
2673  list GIT = GITfanParallel(OC,Q,Qgamma);
2674  fan Sigma = hashesToFan(GIT,OC);
2675  return(Sigma);
2676}
2677example
2678{
2679  echo=2;
2680  intmat Q[3][4] =
2681    1,0,1,0,
2682    0,1,0,1,
2683    0,0,1,1;
2684
2685  gkzFan(Q);
2686}
2687
2688
2689
2690/////////////////////////////////////
2691
2692// Computes all simplicial orbit cones
2693// w.r.t. the 0-ideal:
2694static proc simplicialToricOrbitCones(bigintmat Q){
2695  intvec gam0;
2696  list OC;
2697  cone c;
2698  int r = ncols(Q);
2699  int j;
2700
2701  for(int i = 1; i < 2^r; i++ ){
2702    gam0 = int2face(i,r);
2703
2704    // each simplicial cone is generated by
2705    // exactly nrows(Q) many columns of Q:
2706    if(size(gam0) == nrows(Q)){
2707      bigintmat M[size(gam0)][nrows(Q)];
2708
2709      for(j = 1; j <= size(gam0); j++){
2710        M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]];
2711      }
2712
2713      c = coneViaPoints(M);
2714
2715      if((dimension(c) == nrows(Q)) and (!(listContainsCone(OC, c)))){
2716        OC[size(OC)+1] = c;
2717      }
2718
2719      kill M;
2720    }
2721  }
2722
2723  return(OC);
2724}
2725example
2726{
2727  echo = 2;
2728  bigintmat Q[3][4] =
2729  1,0,1,0,
2730  0,1,0,1,
2731  0,0,1,1;
2732
2733  list OC = simplicialToricOrbitCones(Q);
2734  print(OC);
2735
2736  bigintmat Q[5][15] =
2737    1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,
2738    0,1,0,0,0,1,0,0,0,1,1,1,0,0,0,
2739    0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,
2740    0,0,0,1,0,0,0,1,0,0,1,0,1,0,1,
2741    0,0,0,0,1,0,0,0,1,0,0,1,0,1,1;
2742
2743  list OC = simplicialToricOrbitCones(Q);
2744  print(size(OC));
2745}
2746
2747
2748proc G25Action()
2749"USAGE: G25Action(Q);
2750PURPOSE: Returns a representation of S5 as a subgroup of S10 with the action on the Grassmannian G25.
2751RETURN: list of with elements of type permutation.
2752EXAMPLE: example G25Action; shows an example
2753"
2754{
2755list simplexSymmetryGroup = permutationFromIntvec(intvec( 1 .. 10 )),
2756permutationFromIntvec(intvec( 1, 2, 4, 3, 5, 7, 6, 9, 8, 10 )),
2757permutationFromIntvec(intvec( 1, 3, 2, 4, 6, 5, 7, 8, 10, 9 )),
2758permutationFromIntvec(intvec( 1, 3, 4, 2, 6, 7, 5, 10, 8, 9 )),
2759permutationFromIntvec(intvec( 1, 4, 2, 3, 7, 5, 6, 9, 10, 8 )),
2760permutationFromIntvec(intvec( 1, 4, 3, 2, 7, 6, 5, 10, 9, 8 )),
2761permutationFromIntvec(intvec( 1, 5, 6, 7, 2, 3, 4, 8, 9, 10 )),
2762permutationFromIntvec(intvec( 1, 5, 7, 6, 2, 4, 3, 9, 8, 10 )),
2763permutationFromIntvec(intvec( 1, 6, 5, 7, 3, 2, 4, 8, 10, 9 )),
2764permutationFromIntvec(intvec( 1, 6, 7, 5, 3, 4, 2, 10, 8, 9 )),
2765permutationFromIntvec(intvec( 1, 7, 5, 6, 4, 2, 3, 9, 10, 8 )),
2766permutationFromIntvec(intvec( 1, 7, 6, 5, 4, 3, 2, 10, 9, 8 )),
2767permutationFromIntvec(intvec( 2, 1, 3, 4, 5, 8, 9, 6, 7, 10 )),
2768permutationFromIntvec(intvec( 2, 1, 4, 3, 5, 9, 8, 7, 6, 10 )),
2769permutationFromIntvec(intvec( 2, 3, 1, 4, 8, 5, 9, 6, 10, 7 )),
2770permutationFromIntvec(intvec( 2, 3, 4, 1, 8, 9, 5, 10, 6, 7 )),
2771permutationFromIntvec(intvec( 2, 4, 1, 3, 9, 5, 8, 7, 10, 6 )),
2772permutationFromIntvec(intvec( 2, 4, 3, 1, 9, 8, 5, 10, 7, 6 )),
2773permutationFromIntvec(intvec( 2, 5, 8, 9, 1, 3, 4, 6, 7, 10 )),
2774permutationFromIntvec(intvec( 2, 5, 9, 8, 1, 4, 3, 7, 6, 10 )),
2775permutationFromIntvec(intvec( 2, 8, 5, 9, 3, 1, 4, 6, 10, 7 )),
2776permutationFromIntvec(intvec( 2, 8, 9, 5, 3, 4, 1, 10, 6, 7 )),
2777permutationFromIntvec(intvec( 2, 9, 5, 8, 4, 1, 3, 7, 10, 6 )),
2778permutationFromIntvec(intvec( 2, 9, 8, 5, 4, 3, 1, 10, 7, 6 )),
2779permutationFromIntvec(intvec( 3, 1, 2, 4, 6, 8, 10, 5, 7, 9 )),
2780permutationFromIntvec(intvec( 3, 1, 4, 2, 6, 10, 8, 7, 5, 9 )),
2781permutationFromIntvec(intvec( 3, 2, 1, 4, 8, 6, 10, 5, 9, 7 )),
2782permutationFromIntvec(intvec( 3, 2, 4, 1, 8, 10, 6, 9, 5, 7 )),
2783permutationFromIntvec(intvec( 3, 4, 1, 2, 10, 6, 8, 7, 9, 5 )),
2784permutationFromIntvec(intvec( 3, 4, 2, 1, 10, 8, 6, 9, 7, 5 )),
2785permutationFromIntvec(intvec( 3, 6, 8, 10, 1, 2, 4, 5, 7, 9 )),
2786permutationFromIntvec(intvec( 3, 6, 10, 8, 1, 4, 2, 7, 5, 9 )),
2787permutationFromIntvec(intvec( 3, 8, 6, 10, 2, 1, 4, 5, 9, 7 )),
2788permutationFromIntvec(intvec( 3, 8, 10, 6, 2, 4, 1, 9, 5, 7 )),
2789permutationFromIntvec(intvec( 3, 10, 6, 8, 4, 1, 2, 7, 9, 5 )),
2790permutationFromIntvec(intvec( 3, 10, 8, 6, 4, 2, 1, 9, 7, 5 )),
2791permutationFromIntvec(intvec( 4, 1, 2, 3, 7, 9, 10, 5, 6, 8 )),
2792permutationFromIntvec(intvec( 4, 1, 3, 2, 7, 10, 9, 6, 5, 8 )),
2793permutationFromIntvec(intvec( 4, 2, 1, 3, 9, 7, 10, 5, 8, 6 )),
2794permutationFromIntvec(intvec( 4, 2, 3, 1, 9, 10, 7, 8, 5, 6 )),
2795permutationFromIntvec(intvec( 4, 3, 1, 2, 10, 7, 9, 6, 8, 5 )),
2796permutationFromIntvec(intvec( 4, 3, 2, 1, 10, 9, 7, 8, 6, 5 )),
2797permutationFromIntvec(intvec( 4, 7, 9, 10, 1, 2, 3, 5, 6, 8 )),
2798permutationFromIntvec(intvec( 4, 7, 10, 9, 1, 3, 2, 6, 5, 8 )),
2799permutationFromIntvec(intvec( 4, 9, 7, 10, 2, 1, 3, 5, 8, 6 )),
2800permutationFromIntvec(intvec( 4, 9, 10, 7, 2, 3, 1, 8, 5, 6 )),
2801permutationFromIntvec(intvec( 4, 10, 7, 9, 3, 1, 2, 6, 8, 5 )),
2802permutationFromIntvec(intvec( 4, 10, 9, 7, 3, 2, 1, 8, 6, 5 )),
2803permutationFromIntvec(intvec( 5, 1, 6, 7, 2, 8, 9, 3, 4, 10 )),
2804permutationFromIntvec(intvec( 5, 1, 7, 6, 2, 9, 8, 4, 3, 10 )),
2805permutationFromIntvec(intvec( 5, 2, 8, 9, 1, 6, 7, 3, 4, 10 )),
2806permutationFromIntvec(intvec( 5, 2, 9, 8, 1, 7, 6, 4, 3, 10 )),
2807permutationFromIntvec(intvec( 5, 6, 1, 7, 8, 2, 9, 3, 10, 4 )),
2808permutationFromIntvec(intvec( 5, 6, 7, 1, 8, 9, 2, 10, 3, 4 )),
2809permutationFromIntvec(intvec( 5, 7, 1, 6, 9, 2, 8, 4, 10, 3 )),
2810permutationFromIntvec(intvec( 5, 7, 6, 1, 9, 8, 2, 10, 4, 3 )),
2811permutationFromIntvec(intvec( 5, 8, 2, 9, 6, 1, 7, 3, 10, 4 )),
2812permutationFromIntvec(intvec( 5, 8, 9, 2, 6, 7, 1, 10, 3, 4 )),
2813permutationFromIntvec(intvec( 5, 9, 2, 8, 7, 1, 6, 4, 10, 3 )),
2814permutationFromIntvec(intvec( 5, 9, 8, 2, 7, 6, 1, 10, 4, 3 )),
2815permutationFromIntvec(intvec( 6, 1, 5, 7, 3, 8, 10, 2, 4, 9 )),
2816permutationFromIntvec(intvec( 6, 1, 7, 5, 3, 10, 8, 4, 2, 9 )),
2817permutationFromIntvec(intvec( 6, 3, 8, 10, 1, 5, 7, 2, 4, 9 )),
2818permutationFromIntvec(intvec( 6, 3, 10, 8, 1, 7, 5, 4, 2, 9 )),
2819permutationFromIntvec(intvec( 6, 5, 1, 7, 8, 3, 10, 2, 9, 4 )),
2820permutationFromIntvec(intvec( 6, 5, 7, 1, 8, 10, 3, 9, 2, 4 )),
2821permutationFromIntvec(intvec( 6, 7, 1, 5, 10, 3, 8, 4, 9, 2 )),
2822permutationFromIntvec(intvec( 6, 7, 5, 1, 10, 8, 3, 9, 4, 2 )),
2823permutationFromIntvec(intvec( 6, 8, 3, 10, 5, 1, 7, 2, 9, 4 )),
2824permutationFromIntvec(intvec( 6, 8, 10, 3, 5, 7, 1, 9, 2, 4 )),
2825permutationFromIntvec(intvec( 6, 10, 3, 8, 7, 1, 5, 4, 9, 2 )),
2826permutationFromIntvec(intvec( 6, 10, 8, 3, 7, 5, 1, 9, 4, 2 )),
2827permutationFromIntvec(intvec( 7, 1, 5, 6, 4, 9, 10, 2, 3, 8 )),
2828permutationFromIntvec(intvec( 7, 1, 6, 5, 4, 10, 9, 3, 2, 8 )),
2829permutationFromIntvec(intvec( 7, 4, 9, 10, 1, 5, 6, 2, 3, 8 )),
2830permutationFromIntvec(intvec( 7, 4, 10, 9, 1, 6, 5, 3, 2, 8 )),
2831permutationFromIntvec(intvec( 7, 5, 1, 6, 9, 4, 10, 2, 8, 3 )),
2832permutationFromIntvec(intvec( 7, 5, 6, 1, 9, 10, 4, 8, 2, 3 )),
2833permutationFromIntvec(intvec( 7, 6, 1, 5, 10, 4, 9, 3, 8, 2 )),
2834permutationFromIntvec(intvec( 7, 6, 5, 1, 10, 9, 4, 8, 3, 2 )),
2835permutationFromIntvec(intvec( 7, 9, 4, 10, 5, 1, 6, 2, 8, 3 )),
2836permutationFromIntvec(intvec( 7, 9, 10, 4, 5, 6, 1, 8, 2, 3 )),
2837permutationFromIntvec(intvec( 7, 10, 4, 9, 6, 1, 5, 3, 8, 2 )),
2838permutationFromIntvec(intvec( 7, 10, 9, 4, 6, 5, 1, 8, 3, 2 )),
2839permutationFromIntvec(intvec( 8, 2, 5, 9, 3, 6, 10, 1, 4, 7 )),
2840permutationFromIntvec(intvec( 8, 2, 9, 5, 3, 10, 6, 4, 1, 7 )),
2841permutationFromIntvec(intvec( 8, 3, 6, 10, 2, 5, 9, 1, 4, 7 )),
2842permutationFromIntvec(intvec( 8, 3, 10, 6, 2, 9, 5, 4, 1, 7 )),
2843permutationFromIntvec(intvec( 8, 5, 2, 9, 6, 3, 10, 1, 7, 4 )),
2844permutationFromIntvec(intvec( 8, 5, 9, 2, 6, 10, 3, 7, 1, 4 )),
2845permutationFromIntvec(intvec( 8, 6, 3, 10, 5, 2, 9, 1, 7, 4 )),
2846permutationFromIntvec(intvec( 8, 6, 10, 3, 5, 9, 2, 7, 1, 4 )),
2847permutationFromIntvec(intvec( 8, 9, 2, 5, 10, 3, 6, 4, 7, 1 )),
2848permutationFromIntvec(intvec( 8, 9, 5, 2, 10, 6, 3, 7, 4, 1 )),
2849permutationFromIntvec(intvec( 8, 10, 3, 6, 9, 2, 5, 4, 7, 1 )),
2850permutationFromIntvec(intvec( 8, 10, 6, 3, 9, 5, 2, 7, 4, 1 )),
2851permutationFromIntvec(intvec( 9, 2, 5, 8, 4, 7, 10, 1, 3, 6 )),
2852permutationFromIntvec(intvec( 9, 2, 8, 5, 4, 10, 7, 3, 1, 6 )),
2853permutationFromIntvec(intvec( 9, 4, 7, 10, 2, 5, 8, 1, 3, 6 )),
2854permutationFromIntvec(intvec( 9, 4, 10, 7, 2, 8, 5, 3, 1, 6 )),
2855permutationFromIntvec(intvec( 9, 5, 2, 8, 7, 4, 10, 1, 6, 3 )),
2856permutationFromIntvec(intvec( 9, 5, 8, 2, 7, 10, 4, 6, 1, 3 )),
2857permutationFromIntvec(intvec( 9, 7, 4, 10, 5, 2, 8, 1, 6, 3 )),
2858permutationFromIntvec(intvec( 9, 7, 10, 4, 5, 8, 2, 6, 1, 3 )),
2859permutationFromIntvec(intvec( 9, 8, 2, 5, 10, 4, 7, 3, 6, 1 )),
2860permutationFromIntvec(intvec( 9, 8, 5, 2, 10, 7, 4, 6, 3, 1 )),
2861permutationFromIntvec(intvec( 9, 10, 4, 7, 8, 2, 5, 3, 6, 1 )),
2862permutationFromIntvec(intvec( 9, 10, 7, 4, 8, 5, 2, 6, 3, 1 )),
2863permutationFromIntvec(intvec( 10, 3, 6, 8, 4, 7, 9, 1, 2, 5 )),
2864permutationFromIntvec(intvec( 10, 3, 8, 6, 4, 9, 7, 2, 1, 5 )),
2865permutationFromIntvec(intvec( 10, 4, 7, 9, 3, 6, 8, 1, 2, 5 )),
2866permutationFromIntvec(intvec( 10, 4, 9, 7, 3, 8, 6, 2, 1, 5 )),
2867permutationFromIntvec(intvec( 10, 6, 3, 8, 7, 4, 9, 1, 5, 2 )),
2868permutationFromIntvec(intvec( 10, 6, 8, 3, 7, 9, 4, 5, 1, 2 )),
2869permutationFromIntvec(intvec( 10, 7, 4, 9, 6, 3, 8, 1, 5, 2 )),
2870permutationFromIntvec(intvec( 10, 7, 9, 4, 6, 8, 3, 5, 1, 2 )),
2871permutationFromIntvec(intvec( 10, 8, 3, 6, 9, 4, 7, 2, 5, 1 )),
2872permutationFromIntvec(intvec( 10, 8, 6, 3, 9, 7, 4, 5, 2, 1 )),
2873permutationFromIntvec(intvec( 10, 9, 4, 7, 8, 3, 6, 2, 5, 1 )),
2874permutationFromIntvec(intvec( 10, 9, 7, 4, 8, 6, 3, 5, 2, 1 ));
2875return(simplexSymmetryGroup);
2876}
2877example
2878{
2879  echo = 2;
2880  G25Action();
2881}
2882
2883
2884proc findOrbits(list G, int #)
2885"USAGE: findOrbits(G [,d]); G list of permutations in a subgroup of the symmetric group; d int minimum cardinality of simplices to be considered; if d is not specified all orbits are computed.
2886PURPOSE: Computes the orbit decomposition of the action of G.
2887RETURN: list of intvec.
2888EXAMPLE: example findOrbits; shows an example
2889"
2890{
2891int d;
2892if (size(#)>0){d=#;}
2893int n = size(permutationToIntvec(G[1]));
2894list listOrbits;
2895list finished;
2896int tst;
2897bigint startel;
2898intvec startelbin;
2899list neworbit,neworbitint;
2900int i,posToInsert;
2901bigint nn;
2902while (size(finished)<2^n){
2903   startel=0;
2904   if (size(finished)>0){
2905      tst=0;
2906      while ((tst==0) and (startel+1<=size(finished))){
2907        if (finished[int(startel+1)]<>startel){
2908          tst=1;
2909        } else {
2910          startel=startel+1;
2911        }
2912      }
2913   }
2914   if (startel==0){
2915     neworbit[1]= list();
2916     neworbitint[1]=0;
2917   } else {
2918      startelbin=bigintToBinary(startel,n);
2919       neworbitint=list(startel);
2920       for (i=2;i<=size(G);i++){
2921         nn=binaryToBigint(applyPermutationToIntvec(startelbin,G[i]));
2922         //nn;neworbitint;
2923         //"place to insert";
2924         posToInsert = findPlaceToInsert(neworbitint,nn);
2925         //"pos";posToInsert;
2926         if(posToInsert > 0)
2927         {
2928          //"vorher";neworbitint;
2929          neworbitint = insertToList(neworbitint,nn,posToInsert);
2930          //"nachher";neworbitint;
2931         }
2932       }
2933       for (i=1;i<=size(neworbitint);i++){
2934          neworbit[i]=bigintToBinary(neworbitint[i],n);
2935       }
2936   }
2937   if (size(neworbit[1])>=d){
2938     listOrbits[size(listOrbits)+1]=neworbit;
2939   }
2940   finished=sort(finished+neworbitint)[1];
2941}
2942return(listOrbits);}
2943example
2944{
2945list G = G25Action();
2946list orb = findOrbits(G);
2947for (int i=1;i<=size(orb);i++){orb[i][1];}
2948}
2949
2950
2951
2952static proc GITfanWrapperWithSymmetry(ideal J, intmat Q, list simplexSymmetryGroup){
2953  list orb = findOrbits(simplexSymmetryGroup,nrows(Q));
2954  list simplexOrbitRepresentatives;
2955  for (int i=1;i<=size(orb);i++){simplexOrbitRepresentatives[i]=orb[i][1];}
2956  list afaceOrbitRepresentatives=afaces(J,simplexOrbitRepresentatives);
2957  list fulldimAfaceOrbitRepresentatives=fullDimImages(afaceOrbitRepresentatives,Q);
2958  list afaceOrbits=computeAfaceOrbits(fulldimAfaceOrbitRepresentatives,simplexSymmetryGroup);
2959  list minAfaceOrbits = minimalAfaceOrbits(afaceOrbits);
2960  list listOfOrbitConeOrbits = orbitConeOrbits(minAfaceOrbits,Q);
2961  list listOfMinimalOrbitConeOrbits = minimalOrbitConeOrbits(listOfOrbitConeOrbits);
2962  list Asigma = groupActionOnQImage(simplexSymmetryGroup,Q);
2963  list actionOnOrbitconeIndices = groupActionOnHashes(Asigma,listOfOrbitConeOrbits);
2964  list OClist = listOfOrbitConeOrbits[1];
2965  for (i =2;i<=size(listOfOrbitConeOrbits);i++){
2966    OClist = OClist + listOfOrbitConeOrbits[i];
2967  }
2968  cone mov = coneViaPoints(transpose(Q));
2969  mov = canonicalizeCone(mov);
2970  list Sigma = GITfanParallelSymmetric(OClist, Q, mov, actionOnOrbitconeIndices);
2971  return(Sigma, OClist);
2972}
2973///////////////////////////////////////
2974
2975proc gitCone(ideal a, bigintmat Q, bigintmat w)
2976"USAGE: gitCone(a, Q, w); a: ideal, Q:bigintmat, w:bigintmat
2977PURPOSE: Returns the GIT-cone lambda(w), i.e. the intersection of all
2978orbit cones containing the vector w.
2979NOTE: call this only if you are interested in a single GIT-cone.
2980RETURN: a cone.
2981EXAMPLE: example gitCone; shows an example
2982"
2983{
2984  list OC =  orbitConesOld(a, intmat(Q));
2985  cone lambda = nrows(Q);
2986
2987  for(int i = 1; i <= size(OC); i++)
2988  {
2989    cone c = OC[i];
2990
2991    if(containsInSupport(c, w))
2992    {
2993      lambda = convexIntersectionOld(lambda, c);
2994    }
2995
2996    kill c;
2997  }
2998
2999  return(lambda);
3000}
3001example
3002{
3003  echo=2;
3004  intmat Q[3][4] =
3005    1,0,1,0,
3006    0,1,0,1,
3007    0,0,1,1;
3008
3009  ring R = 0,T(1..4),dp;
3010  ideal a = 0;
3011
3012  bigintmat w[1][3] = 3,3,1;
3013  cone lambda = gitCone(a, Q, w);
3014  rays(lambda);
3015
3016  bigintmat w2[1][3] = 1,1,1;
3017  cone lambda2 = gitCone(a, Q, w2);
3018  rays(lambda2);
3019}
3020
3021proc orbitConesOld(ideal a, bigintmat Q, list #)
3022"USAGE:  orbitConesOld(a, Q, b, c); a: ideal, Q: bigintmat, b: int, c: int
3023PURPOSE: Returns a list consisting of all cones Q(gam0) where gam0 is an a-face.
3024         Moreover, it is possible to specify a dimensional bound b,
3025         upon which only cones of that dimension and above are returned.
3026         Lastly, as the computation is parallizable, one can specify c,
3027         the number of cores to be used by the computation.
3028RETURN:  a list of cones
3029EXAMPLE: example orbitConesOld; shows an example
3030"
3031{
3032  list AF;
3033
3034  if((size(#) > 1) and (typeof(#[2]) == "int"))
3035  {
3036    AF = afaces(a, nrows(Q), #[2]);
3037  }
3038  else
3039  {
3040    AF = afaces(a, nrows(Q));
3041  }
3042
3043  int dimensionBound = 0;
3044  if((size(#) > 0) and (typeof(#[1]) == "int"))
3045  {
3046    dimensionBound = #[1];
3047  }
3048
3049  list OC;
3050  intvec gam0;
3051  int j;
3052
3053  for(int i = 1; i <= size(AF); i++)
3054  {
3055    gam0 = AF[i];
3056
3057    if(gam0 == 0)
3058    {
3059      bigintmat M[1][nrows(Q)];
3060    }
3061    else
3062    {
3063      bigintmat M[size(gam0)][nrows(Q)];
3064      for (j = 1; j <= size(gam0); j++)
3065      {
3066        M[j,1..ncols(M)] = Q[1..nrows(Q),gam0[j]];
3067      }
3068    }
3069    cone c = coneViaPoints(M);
3070
3071    if((dimension(c) >= dimensionBound) and (!(listContainsCone(OC, c))))
3072    {
3073      OC[size(OC)+1] = c;
3074    }
3075
3076    kill M, c;
3077  }
3078
3079  return(OC);
3080}
3081example
3082{
3083  echo=2;
3084  intmat Q[3][4] =
3085    1,0,1,0,
3086    0,1,0,1,
3087    0,0,1,1;
3088
3089  ring R = 0,T(1..4),dp;
3090  ideal a = 0;
3091
3092  orbitConesOld(a, Q);
3093}
3094
Note: See TracBrowser for help on using the repository browser.