orbit_fkts.c

/* last change: 26.09.00 by Oliver Heidbuechel */

#include <typedef.h>
#include <matrix.h>
#include <bravais.h>
#include <base.h>
#include <graph.h>
#include <sort.h>
#include <zass.h>
#include <longtools.h>


/* --------------------------------------------------------------------------- */
/* Add two cocycles in standard representation C_d1 x ... x C_dn and           */
/* transform the result into standard representation                           */
/* --------------------------------------------------------------------------- */
matrix_TYP *add_mod_D(matrix_TYP *A,
                      matrix_TYP *B,
                      matrix_TYP *D,
                      int first,
                      int diff)
{
   int i;

   matrix_TYP *sol;


   sol = init_mat(diff, 1, "");
   for (i = 0; i < diff; i++){
      sol->array.SZ[i][0] = (A->array.SZ[i][0] + B->array.SZ[i][0]) %
                                       D->array.SZ[i + first][i + first];
      if (sol->array.SZ[i][0] < 0)
         sol->array.SZ[i][0] += D->array.SZ[i + first][i + first];
   }

   return(sol);
}



/* --------------------------------------------------------------------------- */
/* return min(a,b)                                                             */
/* --------------------------------------------------------------------------- */
static int min(int a,
               int b)
{
   if (a < b)
      return(a);
   else
      return(b);
}



/* --------------------------------------------------------------------------- */
/* which affine classes are in the image                                       */
/* --------------------------------------------------------------------------- */
/* names: names for the affine classes                                         */
/* aff_no: total number of affine classes                                      */
/* image_gen: generating system of the image                                   */
/* D: the second return value of cohomology                                    */
/* anz: save the number of affine classes in the image here                    */
/* --------------------------------------------------------------------------- */
int *aff_classes_in_image(MP_INT *names,
                          int aff_no,
                          int coho_size,
                          matrix_TYP **orbit,
                          matrix_TYP **image_gen,
                          int gen_no,
                          matrix_TYP *D,
                          int *anz)
{
   int i, j, k, pos, /* m, */
       first, last, diff,
       insg, flag = 0,
      *list, *aff_list;

   matrix_TYP *tmp;

   MP_INT val;


   mpz_init(&val);
   aff_list = (int *)calloc(aff_no, sizeof(int));
   list = (int *)calloc(coho_size, sizeof(int));
   list[0] = -1;
   anz[0] = 1;
   aff_list[0] = 1;
   for (first = 0; first < D->cols && D->array.SZ[first][first] == 1; first++);
   for (last = first; last < D->cols && D->array.SZ[last][last] != 0; last++);
   diff = last - first;
   if (orbit[0] == NULL)
      orbit[0] = init_mat(diff, 1, "");
   /* m = coho_size; */

   /* orbit algorithm */
   for (i = 0; i < coho_size && anz[0] < aff_no; i++){
      if (list[i] == -1){
         for (j = 0; j < gen_no; j++){
            tmp = add_mod_D(orbit[i], image_gen[j], D, first, diff);
            valuation(tmp, D, &val);
            pos = mpz_get_ui(&val);
            if (list[pos] == 0){
               /* new element in the orbit */
               list[pos] = -1;
               if (orbit[pos] == NULL){
                  orbit[pos] = tmp;
                  tmp = NULL;
               }
               for (k = 0; k < aff_no; k++){
                  if (mpz_cmp(&val, &names[k]) == 0){
                     /* found representative of the k-th affine class */
                     aff_list[k] = 1;
                     anz[0]++;
                     break;
                  }
               }
            }
            if (tmp != NULL)
               free_mat(tmp);
	    /*
            if (list[pos] == -1){
               m = min(pos, m);
            }
	    */
         }
         list[i] = 1;
	 /*
	 if (m == i){
	    m = coho_size;
	 }
	 else{
            i = m - 1;
	 }
	 */
	 i = -1;
      }
   }

   /* clean */
   mpz_clear(&val);
   free(list);

   return(aff_list);
}





/* --------------------------------------------------------------------------- */
/* calculate all elements in the group generated by M                          */
/* in the cohomology group given by the second return value of cohomology      */
/* the elements of M have to be in in this cohomology group                    */
/* return list, which number/element of the cohomology group is in our group   */
/* --------------------------------------------------------------------------- */
/* M: generating system                                                        */
/* D: the second return value of cohomology                                    */
/* anz: save the order of the group here                                       */
/* elements: save the elements of the group (at the correct place in the       */
/*             cohomology group) here                                          */
/* coho_size: size of the cohomology group                                     */
/* --------------------------------------------------------------------------- */
int *aufspannen(int coho_size,
                matrix_TYP **elements,
                matrix_TYP **M,
                int gen_no,
                matrix_TYP *D,
                int *anz)
{
   int i, j, pos, /* m, */
       first, last, diff,
       *list;

   matrix_TYP *tmp;

   MP_INT val;


   mpz_init(&val);
   list = (int *)calloc(coho_size, sizeof(int));
   list[0] = -1;
   anz[0] = 1;
   for (first = 0; first < D->cols && D->array.SZ[first][first] == 1; first++);
   for (last = first; last < D->cols && D->array.SZ[last][last] != 0; last++);
   diff = last - first;
   if (elements[0] == NULL)
      elements[0] = init_mat(diff, 1, "");
   /* m = coho_size; */

   /* orbit algorithm */
   for (i = 0; i < coho_size && anz[0] < coho_size; i++){
      if (list[i] == -1){
         for (j = 0; j < gen_no; j++){
            tmp = add_mod_D(elements[i], M[j], D, first, diff);
            valuation(tmp, D, &val);
            pos = mpz_get_ui(&val);
            if (list[pos] == 0){
               /* new element in the orbit */
               anz[0]++;
               list[pos] = -1;
               if (elements[pos] == NULL){
                  elements[pos] = tmp;
                  tmp = NULL;
               }
            }
            if (tmp != NULL)
               free_mat(tmp);
	    /*
            if (list[pos] == -1){
               m = min(pos, m);
            }
	    */
         }
	 list[i] = 1;
	 /*
	 if (m == i){
	    m = coho_size;
	 }
	 else{
	    i = m - 1;
	 }
	 */
	 i = -1;
      }
   }

   for (i = 0; i < coho_size; i++){
      if (list[i] == -1)
         list[i] = 1;
   }

   /* clean */
   mpz_clear(&val);

   return(list);
}



/* ----------------------------------------------------------------------------- */
static void mod(int *a,
                int b)
{
   a[0] %= b;
   if (a[0] < 0)
      a[0] += b;
}



/* ----------------------------------------------------------------------------- */
static int search_next(int *list,
                       int *ker_list,
                       int coho_size)
{
   int i;

   if (ker_list != NULL){
      for (i = 1; i < coho_size; i++){
         if (ker_list[i] == 1 && list[i] == 0)
            return(i);
      }
   }
   else{
      for (i = 1; i < coho_size; i++){
         if (list[i] == 0)
            return(i);
      }
   }

   return(-1);
}


/* ----------------------------------------------------------------------------- */
/* calculate the representation of S on H^1(...)/Ker(phi)                        */
/* ----------------------------------------------------------------------------- */
matrix_TYP **new_representation(matrix_TYP **S,
                                int S_no,
                                H1_mod_ker_TYP H1_mod_ker,
                                matrix_TYP *A)
{
   int i, j, k, d, s,
       erz_no, A_first;

   matrix_TYP **rep,
               *tmp, *temp;


   erz_no = H1_mod_ker.erz_no;
   rep = (matrix_TYP **)calloc(S_no, sizeof(matrix_TYP *));
   for (A_first = 0; A_first < A->cols && A->array.SZ[A_first][A_first] == 1; A_first++);

   for (i = 0; i < S_no; i++){
      rep[i] = init_mat(erz_no, erz_no, "");
      for (j = 0; j < erz_no; j++){

         /* bilde Basiselement ab mit Erzeuger von S1 */
         tmp = mat_mul(S[i], H1_mod_ker.M[j]);
         for (k = 0; k < tmp->rows; k++){
            tmp->array.SZ[k][0] %= A->array.SZ[k + A_first][k + A_first];
            if (tmp->array.SZ[k][0] < 0)
               tmp->array.SZ[k][0] += A->array.SZ[k + A_first][k + A_first];
         }

	 /* Berechne Standardform als Element von H^1/Ker */
         temp = mat_mul(H1_mod_ker.i, tmp);
         free_mat(tmp);
         for (k = 0; k < temp->rows; k++){
            temp->array.SZ[k][0] %= H1_mod_ker.D->array.SZ[k][k];
            if (temp->array.SZ[k][0] < 0)
               temp->array.SZ[k][0] += H1_mod_ker.D->array.SZ[k][k];
         }

	 /* trage in Matrix ein */
         for (k = 0; k < erz_no; k++){
            rep[i]->array.SZ[k][j] = temp->array.SZ[k + H1_mod_ker.D_first][0];
         }
         free_mat(temp);
      }
   }
   return(rep);
}



/* ----------------------------------------------------------------------------- */
/* from the H^1/Ker representation to the H^1 representation                     */
/* ----------------------------------------------------------------------------- */
static matrix_TYP *back_to_H1(matrix_TYP *elem,
                              H1_mod_ker_TYP Hk)
{
   int i;

   matrix_TYP *mat;

   rational eins, zahl;

   eins.z = eins.n = zahl.n = 1;


   mat = copy_mat(Hk.M[0]);
   for (i = 0; i < mat->rows; i++){
      mat->array.SZ[i][0] *= elem->array.SZ[0][0];
   }

   for (i = 1; i < Hk.erz_no; i++){
      if (elem->array.SZ[i][0] != 0){
         zahl.z = elem->array.SZ[i][0];
         mat = mat_addeq(mat, Hk.M[i], eins, zahl);
      }
   }

   return(mat);
}



/* ----------------------------------------------------------------------------- */
/* orbit algorithm: act with S on H^1(...)/Ker(phi)                              */
/* in orbit 0 is only the 0 element                                              */
/* save WORDS for the stabilizers of ksi + Ker(phi) for a representative in      */
/* each orbit                                                                    */
/* start = NULL: berechne alle Orbits                                            */
/* start != NULL: berechne nur den Orbit von start                               */
/* ----------------------------------------------------------------------------- */
matrix_TYP ***H1_mod_ker_orbit_alg(H1_mod_ker_TYP H1_mod_ker,
                                   matrix_TYP **S,
                                   int S_no,
                                   int *anz,
                                   int **length,
                                   int ****WORDS,
                                   int **WORDS_no,
				   matrix_TYP *start)
{
   matrix_TYP *D,
              *tmp,
             **elem,
            ***rep;

   int i, j, k, /* m, */ i__,
       pos, number, addmem,
       *list, *smallest, **words,
       coho_size, erz_no, counter = 0;

   MP_INT val, cohom_size, MP_i;


   mpz_init(&val);
   erz_no = H1_mod_ker.erz_no;
   D = init_mat(erz_no, erz_no, "");
   for (i = 0; i < erz_no; i++){
      D->array.SZ[i][i] = H1_mod_ker.D->array.SZ[i + H1_mod_ker.D_first][i + H1_mod_ker.D_first];
   }
   cohom_size = cohomology_size(D);
   coho_size = mpz_get_ui(&cohom_size);
   elem = (matrix_TYP **)calloc(coho_size, sizeof(matrix_TYP *));
   list = (int *)calloc(coho_size, sizeof(int));
   length[0] = (int *)calloc(coho_size + 1, sizeof(int));
   smallest = (int *)calloc(coho_size + 1, sizeof(int));
   WORDS[0] = (int ***)calloc(coho_size + 1, sizeof(int **));
   WORDS_no[0] = (int *)calloc(coho_size + 1, sizeof(int));

   if (start == NULL){
      if (coho_size > 1)
         i = 1;
      else
         i = -1;
      anz[0] = 1;
      counter = length[0][0] = 1;
      smallest[0] = 0;
   }
   else{
      anz[0] = 1;
      for (k = 0; k < erz_no; k++)
         mod(start->array.SZ[k], D->array.SZ[k][k]);
      valuation(start, D, &val);
      i = mpz_get_ui(&val);
   }


   while (i != -1){
      WORDS[0][anz[0]] = (int **)calloc(1024, sizeof(int *));
      words = (int **)calloc(coho_size, sizeof(int *));
      list[i] = -anz[0];
      smallest[anz[0]] = i;
      length[0][anz[0]] = 1;
      mpz_init_set_si(&MP_i, i);
      elem[i] = reverse_valuation(&MP_i, D);
      mpz_clear(&MP_i);
      /* m = coho_size; */
      for (; i < coho_size && counter < coho_size; i++){
         if (list[i] == -anz[0]){
            for (j = 0; j < S_no; j++){
               tmp = mat_mul(S[j], elem[i]);
               for (k = 0; k < erz_no; k++)
                  mod(tmp->array.SZ[k], D->array.SZ[k][k]);
               valuation(tmp, D, &val);
               pos = mpz_get_ui(&val);

               /* paranoia test */
               if (list[pos] != 0 && abs(list[pos]) != anz[0]){
                  fprintf(stderr, "ERROR in H1_mod_ker_orbit_alg!\n");
                  exit(8);
               }

               if (list[pos] == 0){
                  /* new element in the orbit */
                  /*words[pos] = (int *)calloc(MIN_SPEICHER, sizeof(int));*/
                  if (i != smallest[anz[0]]){
                     words[pos] = (int *)calloc(words[i][0] + 2, sizeof(int)); /*N*/
                     /*memcpy(words[pos], words[i], (words[i][0] + 1) * sizeof(int));*/
		     for (i__ = 0; i__ < words[i][0] + 1; i__++){
		        words[pos][i__] = words[i][i__];
		     }
                     words[pos][words[i][0] + 1] = -(j+1);
                     words[pos][0]++;
                  }
                  else{
		     words[pos] = (int *)calloc(2, sizeof(int)); /*N*/
                     words[pos][0] = 1;
                     words[pos][1] = -(j+1);
                  }
                  list[pos] = -anz[0];
                  counter++;
                  length[0][anz[0]]++;
                  elem[pos] = tmp;
                  tmp = NULL;
               }
               else{
                  /* we got this element already */
		  if (WORDS_no[0][anz[0]] % 1024 == 0){ /*N*/
		     WORDS[0][anz[0]] = (int **)realloc(WORDS[0][anz[0]],
		                        (1024 + WORDS_no[0][anz[0]]) * sizeof(int *));
                  }

                  /*WORDS[0][anz[0]][WORDS_no[0][anz[0]]] = (int *)calloc(MIN_SPEICHER, sizeof(int));*/
		  if (pos != smallest[anz[0]]){ /*N*/
		     addmem = words[pos][0];
		  }
		  else{
		     addmem = 0;
		  }
                  if (i != smallest[anz[0]]){
                     WORDS[0][anz[0]][WORDS_no[0][anz[0]]] = (int *)calloc(words[i][0] + 2 + addmem, sizeof(int));/*N*/
                     /*memcpy(WORDS[0][anz[0]][WORDS_no[0][anz[0]]], words[i], (words[i][0] + 1) * sizeof(int)); */
		     for (i__ = 0; i__ < words[i][0] + 1; i__++){
                        WORDS[0][anz[0]][WORDS_no[0][anz[0]]][i__] = words[i][i__];
		     }
                     WORDS[0][anz[0]][WORDS_no[0][anz[0]]][words[i][0] + 1] = -(j+1);
                     WORDS[0][anz[0]][WORDS_no[0][anz[0]]][0]++;
                  }
                  else{
                     WORDS[0][anz[0]][WORDS_no[0][anz[0]]] = (int *)calloc(2 + addmem, sizeof(int));/*N*/
                     WORDS[0][anz[0]][WORDS_no[0][anz[0]]][0] = 1;
                     WORDS[0][anz[0]][WORDS_no[0][anz[0]]][1] = -(j+1);
                  }
                  if (pos != smallest[anz[0]]){
                     number = WORDS[0][anz[0]][WORDS_no[0][anz[0]]][0] + words[pos][0];
                     for (k = 0; k < words[pos][0]; k++){
                        WORDS[0][anz[0]][WORDS_no[0][anz[0]]][number - k] = -words[pos][k + 1];
                     }
                     WORDS[0][anz[0]][WORDS_no[0][anz[0]]][0] = number;
                  }
                  WORDS_no[0][anz[0]]++;


                  free_mat(tmp);
               }
	       /*
               if (list[pos] == -anz[0]){
                  m = min(pos, m);
               }
	       */
            }
            list[i] = anz[0];
	    /*
	    if (m == i){
	       m = coho_size;
	    }
	    else{
               i = m - 1;
	    }
	    */
	    i = -1;
         }
      }
      if (start == NULL){
         i = search_next(list, NULL, coho_size);
      }
      else{
         i = -1;
      }
      anz[0]++;
      for (j = 0; j < coho_size; j++)
         if (words[j] != NULL)
            free(words[j]);
      free(words);
   }


   /* representatives */
   rep = (matrix_TYP ***)calloc(anz[0], sizeof(matrix_TYP **));
   for (i = 1; i < anz[0]; i++){
      rep[i] = (matrix_TYP **)calloc(length[0][i], sizeof(matrix_TYP *));
      counter = 0;
      for (j = 0; j < coho_size; j++){
         if ((list[j] == i || list[j] == -i) && counter < length[0][i]){
            rep[i][counter] = back_to_H1(elem[j], H1_mod_ker);
            counter++;
         }
      }
   }

   /* clean */
   mpz_clear(&cohom_size);
   mpz_clear(&val);
   free_mat(D);
   free(list);
   free(smallest);
   for (i = 1; i < coho_size; i++)
      if (elem[i] != NULL)
         free_mat(elem[i]);
   free(elem);


   return(rep);
}


/* ----------------------------------------------------------------------------- */
/* orbitalgorithm on a subgroup of H^1                                           */
/* ker_elements: list of elements in the subgroup                                */
/* H^1 is given by D, the second return value of cohomology                      */
/* coho_size: size of H^1                                                        */
/* N: operating group                                                            */
/* no: number of elements in N                                                   */
/* ker_list: list, which element of the cohomology group is in the subgroup      */
/* ker_order: order of the subgroup                                              */
/* anz: save the number of orbits here                                           */
/* length: save the length of an orbit here                                      */
/* ----------------------------------------------------------------------------- */
/* in orbit 0 is only the 0-element 						 */
/* ----------------------------------------------------------------------------- */
matrix_TYP **orbit_ker(matrix_TYP **ker_elements,
                       int ker_order,
                       matrix_TYP *D,
                       matrix_TYP **N,
                       int no,
                       int coho_size,
                       int *ker_list,
                       int *anz,
                       int **length)
{
   int i, j, k, /* m, */ pos,
       first, last, diff, counter,
       ttt = 0,
      *list, *smallest;

   matrix_TYP **orbit_rep,
               *tmp;

   MP_INT val;


   for (first = 0; first < D->cols && D->array.SZ[first][first] == 1; first++);
   for (last = first; last < D->cols && D->array.SZ[last][last] != 0; last++);
   diff = last - first;
   length[0] = (int *)calloc(coho_size, sizeof(int));
   anz[0] = counter = length[0][0] = 1;
   for (i = 1; i < coho_size; i++){
      if (ker_list[i] == 1)
         break;
   }

   /* only 0-element in the subgroup */
   if (i == coho_size){
      orbit_rep = (matrix_TYP **)calloc(1, sizeof(matrix_TYP *));
      orbit_rep[0] = init_mat(diff, 1, "");
      return(orbit_rep);
   }

   mpz_init(&val);
   list = (int *)calloc(coho_size, sizeof(int));
   smallest = (int *)calloc(coho_size, sizeof(int));

   /* orbit algorithm */
   while (i != -1){
      list[i] = -anz[0];
      smallest[anz[0]] = i;
      length[0][anz[0]] = 1;
      /* m = coho_size; */
      for (; i < coho_size && counter < ker_order; i++){
         if (list[i] == -anz[0]){
            for (j = 0; j < no; j++){
               tmp = mat_mul(N[j], ker_elements[i]);
               for (k = first; k < last; k++)
                  mod(tmp->array.SZ[k-first], D->array.SZ[k][k]);
               valuation(tmp, D, &val);
               pos = mpz_get_ui(&val);

               /* paranoia test */
               if (ker_list[pos] != 1 || (list[pos] != 0 && abs(list[pos]) != anz[0])){
                  fprintf(stderr, "ERROR in orbit_ker!\n");
                  exit(7);
               }

               if (list[pos] == 0){
                  /* new element in the orbit */
                  list[pos] = -anz[0];
                  counter++;
                  length[0][anz[0]]++;
               }
               free_mat(tmp);
	       /*
               if (list[pos] == -anz[0]){
                  m = min(pos, m);
               }
	       */
            }
            list[i] = anz[0];
	    /*
	    if (m == i){
	       m = coho_size;
	    }
	    else{
               i = m - 1;
	    }
	    */
	    i = -1;
         }
      }
      i = search_next(list, ker_list, coho_size);
      anz[0]++;
   }

   /* save representatives */
   orbit_rep = (matrix_TYP **)calloc(anz[0], sizeof(matrix_TYP *));
   for (i = 0; i < anz[0]; i++){
      orbit_rep[i] = copy_mat(ker_elements[smallest[i]]);
      if (GRAPH_DEBUG)
         ttt += length[0][i];
   }

   if (GRAPH_DEBUG){
      if (ttt != ker_order){
         fprintf(stderr, "NEIN!!!!!!!!!!!!\n");
         exit(5);
      }
   }

   /* clean */
   mpz_clear(&val);
   free(smallest);
   free(list);

   return(orbit_rep);
}




/* ----------------------------------------------------------------------------- */
/* orbitalgorithm on a ksi + ker, where ker is a subgroup of H^1                 */
/* ker_elements: list of elements in the subgroup                                */
/* H^1 is given by D, the second return value of cohomology                      */
/* coho_size: size of H^1                                                        */
/* N: operating group (IS CHANGED !!!)                                           */
/* no: number of elements in N                                                   */
/* ker_list: list, which element of the cohomology group is in the subgroup      */
/* ker_order: order of the subgroup                                              */
/* anz: save the number of orbits here                                           */
/* length: save the length of an orbit here                                      */
/* ----------------------------------------------------------------------------- */
/* in orbit 0 is only the 0-element 						 */
/* ----------------------------------------------------------------------------- */
matrix_TYP **orbit_ksi_plus_ker(matrix_TYP *ksi,
                                matrix_TYP **ker_elements,
                                int ker_order,
                                matrix_TYP *D,
                                matrix_TYP **N,
                                int no,
                                int coho_size,
                                int *ker_list,
                                int *anz,
                                int **length)
{
   int i, j, k, /* m, */ rows, first, last, diff, pos,
       ttt = 0,
       *list, *smallest, counter;

   matrix_TYP *tmp, **orbit_rep;

   rational eins, minuseins;

   MP_INT val;



   rows = ker_elements[0]->rows - 1;

   if (ker_order > 1){
      eins.z = eins.n = minuseins.n = 1;
      minuseins.z = -1;
      for (first = 0; first < D->cols && D->array.SZ[first][first] == 1; first++);
      for (last = first; last < D->cols && D->array.SZ[last][last] != 0; last++);
      diff = last - first;

      for (i = 0; i < no; i++){
         tmp = mat_mul(N[i], ksi);
         mat_addeq(tmp, ksi, eins, minuseins);
         for (j = 0; j < rows; j++){
            tmp->array.SZ[j][0] %= D->array.SZ[j + first][j + first];
            if (tmp->array.SZ[j][0] < 0){
               tmp->array.SZ[j][0] += D->array.SZ[j + first][j + first];
            }
         }
         real_mat(N[i], rows + 1, rows);
         real_mat(N[i], rows + 1, rows + 1);
         N[i]->array.SZ[rows][rows] = 1;
         for (j = 0; j < rows; j++){
            N[i]->array.SZ[j][rows] = tmp->array.SZ[j][0];
         }
         free_mat(tmp);
      }

      mpz_init(&val);
      list = (int *)calloc(coho_size, sizeof(int));
      smallest = (int *)calloc(coho_size, sizeof(int));
      length[0] = (int *)calloc(coho_size, sizeof(int));
      counter = anz[0] = 1;
      i = 0;

      /* orbit algorithm */
      while (i != -1){
         list[i] = -anz[0];
         smallest[anz[0] - 1] = i;
         length[0][anz[0] - 1] = 1;
         /* m = coho_size; */
         for (; i < coho_size && counter < ker_order; i++){
            if (list[i] == -anz[0]){
               for (j = 0; j < no; j++){
                  tmp = mat_mul(N[j], ker_elements[i]);
                  for (k = first; k < last; k++)
                     mod(tmp->array.SZ[k-first], D->array.SZ[k][k]);
                  valuation(tmp, D, &val);
                  pos = mpz_get_ui(&val);

                  /* paranoia test */
                  if (ker_list[pos] != 1 || (list[pos] != 0 && abs(list[pos]) != anz[0])){
                     fprintf(stderr, "ERROR in orbit_ksi_plus_ker!\n");
                     exit(7);
                  }

                  if (list[pos] == 0){
                     /* new element in the orbit */
                     list[pos] = -anz[0];
                     counter++;
                     length[0][anz[0] - 1]++;
                  }
                  free_mat(tmp);
		  /*
                  if (list[pos] == -anz[0]){
                     m = min(pos, m);
                  }
		  */
               }
               list[i] = anz[0];
	       /*
	       if (m == i){
	          m = coho_size;
	       }
	       else{
                  i = m - 1;
               }
	       */
	       i = -1;
            }
         }
         i = search_next(list, ker_list, coho_size);
         anz[0]++;
      }

      /* save representatives */
      anz[0]--;
      orbit_rep = (matrix_TYP **)calloc(anz[0], sizeof(matrix_TYP *));
      for (i = 0; i < anz[0]; i++){
         orbit_rep[i] = copy_mat(ker_elements[smallest[i]]);
         real_mat(orbit_rep[i], rows, 1);
         if (GRAPH_DEBUG){
            ttt += length[0][i];
         }
      }

      if (GRAPH_DEBUG){
         if (ttt != ker_order){
            fprintf(stderr, "NEIN!!!!!!!!!!!!\n");
            exit(5);
         }
      }

      /* clean */
      mpz_clear(&val);
      free(smallest);
      free(list);
   }
   else{
      /* trivial part */
      length[0] = (int *)calloc(1, sizeof(int));
      length[0][0] = anz[0] = 1;
      orbit_rep = (matrix_TYP **)calloc(1, sizeof(matrix_TYP *));
      orbit_rep[0] = init_mat(rows, 1, "");
   }

   return(orbit_rep);
}



/* -------------------------------------------------------------------- */
static int search_pos(matrix_TYP *mat, matrix_TYP **menge, int anz)
{
   int i;

   for (i = 0; i < anz; i++){
      if (cmp_mat(mat, menge[i]) == 0)
         return(i);
   }
   return(-1);
}



/* -------------------------------------------------------------------- */
/* orbit_on_lattices	                                                */
/* Let R be a (affine) space group and P its pointgroup.                */
/* Calculation of the orbits of the P-invariant max. sublattices or min.*/
/* superlattices (saved in gitter) under the operation of the linear    */
/* part of the affine normalizer of R.                                  */
/* return number of orbits                                              */
/* -------------------------------------------------------------------- */
/* gitter: lattices in normal form		                    	*/
/* gitter_no: number of lattices (>= 1 !!!)                             */
/* N: affine_normalizer of R (linear part)                              */
/* N_no: number of generators of the normalizer                      	*/
/* list: list, which lattice is in which orbit                          */
/* length: list with the lengthes of the orbits                         */
/* smallest: for each orbit, lowest number of the elements in it        */
/* conj: if conj, save conjugating matrices here                        */
/* -------------------------------------------------------------------- */
int orbit_on_lattices(matrix_TYP **gitter,
                      int gitter_no,
                      matrix_TYP **N,
                      int N_no,
                      int *list,
                      int *length,
                      int *smallest,
                      matrix_TYP **conj)
{
   int i, j, m, pos, counter, anz, dim;

   matrix_TYP *tmp;


   dim = gitter[0]->rows;
   anz = i = counter = 0;

   /* orbit algorithm */
   while (i != -1){
      anz++;
      counter++;
      list[i] = -anz;
      if (conj)
         conj[i] = init_mat(dim, dim, "1");
      length[anz] = 1;
      smallest[anz] = i;
      for (; i < gitter_no && counter < gitter_no; i++){
         if (list[i] == -anz){
            m = gitter_no;
            for (j = 0; j < N_no; j++){
               tmp = mat_mul(N[j], gitter[i]);
               long_col_hnf(tmp);
               pos = search_pos(tmp, gitter, gitter_no);

               /* paranoia test */
               if (pos == -1 || (list[pos] != 0 && abs(list[pos]) != anz)){
                  fprintf(stderr, "ERROR in operiere_gitter!\n");
                  exit(7);
               }

               if (list[pos] == 0){
                  /* new element in the orbit */
                  list[pos] = -anz;
                  counter++;
                  length[anz]++;
                  if (conj){
                     conj[pos] = mat_mul(N[j], conj[i]);
                  }
               }
               free_mat(tmp);
               if (list[pos] == -anz){
                  m = min(pos, m);
               }
            }
            list[i] = anz;
            i = m - 1;
         }
      }
      i = search_next(list, NULL, gitter_no);
   }

   for (i = 0; i < gitter_no; i++){
      if (list[i] < 0) list[i] *= (-1);
   }

   return(anz);
}



/* -------------------------------------------------------------------- */
static matrix_TYP *add_mod_lattice(matrix_TYP *A,
                                   matrix_TYP *B,
                                   matrix_TYP *D,
                                   matrix_TYP *Ti,
                                   matrix_TYP *T,
                                   int dim,
                                   int erz_no)
{
   int i, j;

   rational eins;

   matrix_TYP *C, *tmp, *tmp2;


   eins.z = eins.n = 1;

   C = mat_add(A, B, eins, eins);
   tmp = mat_mul(Ti, C);
   for (i = 0; i < erz_no; i++){
      for (j = 0; j < dim; j++){
         tmp->array.SZ[i * dim + j][0] %= D->array.SZ[j][j];
         if (tmp->array.SZ[i * dim + j][0] < 0)
            tmp->array.SZ[i * dim + j][0] += D->array.SZ[j][j];
      }
   }

   tmp2 = mat_mul(T, tmp);
   free_mat(tmp);
   free_mat(C);
   return(tmp2);
}



/* -------------------------------------------------------------------- */
/* calculate information about Ker phi | B^1 				*/
/* -------------------------------------------------------------------- */
matrix_TYP **kernel_factor_fct(matrix_TYP **translationen,
                               int translanz,
                               int erz_no,
                               matrix_TYP *lattice,
                               int *anz)
{
   int i, j, dim, pos, size = 1024;

   matrix_TYP *tmp, **elem, *D, *R, *Li, *Li_d, *L_d, *L;


   anz[0] = 1;
   dim = lattice->rows;
   elem = (matrix_TYP **)calloc(size, sizeof(matrix_TYP *));
   elem[0] = init_mat(erz_no * dim, 1, "");
   L = init_mat(dim, dim, "1");
   R = init_mat(dim, dim, "1");
   D = long_elt_mat(L, lattice, R);
   Li = mat_inv(L);
   L_d = matrix_on_diagonal(L, erz_no);
   Li_d = matrix_on_diagonal(Li, erz_no);

   for (i = 0; i < anz[0]; i++){
      for (j = 0; j < translanz; j++){
         tmp = add_mod_lattice(elem[i], translationen[j], D, L_d, Li_d, dim, erz_no);
         pos = search_pos(tmp, elem, anz[0]);
         if (pos == -1){
            if (anz[0] >= size){
               size += 1024;
               elem = realloc(elem, size * sizeof(matrix_TYP *));
            }
            elem[anz[0]] = tmp;
            tmp = NULL;
            anz[0]++;
         }
         else{
            free_mat(tmp);
         }
      }
   }

   /* clean */
   free_mat(D);
   free_mat(R);
   free_mat(L);
   free_mat(Li);
   free_mat(L_d);
   free_mat(Li_d);

   return(elem);
}