Changeset e653c9 in git


Ignore:
Timestamp:
Jul 7, 2017, 2:46:14 PM (7 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
fb98c8dd692771969a7b893425029e42cd04a51d
Parents:
67f9e8e1fc9856b995aae2d27c9a3848419e3d13
Message:
doc: cring
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/gitfan.lib

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