delete 'old' folder (still available in old_experiments branch together with other things)
This commit is contained in:
		
							
								
								
									
										312
									
								
								old/coxeter.c
									
									
									
									
									
								
							
							
						
						
									
										312
									
								
								old/coxeter.c
									
									
									
									
									
								
							@@ -1,312 +0,0 @@
 | 
			
		||||
#include "coxeter.h"
 | 
			
		||||
#include "queue.h"
 | 
			
		||||
 | 
			
		||||
#include <memory.h>
 | 
			
		||||
#include <gsl/gsl_linalg.h>
 | 
			
		||||
#include <gsl/gsl_blas.h>
 | 
			
		||||
 | 
			
		||||
static void applyGenerator(int i, double *x, double *y, int rank, double **schlaefliMatrix)
 | 
			
		||||
{
 | 
			
		||||
  memcpy(y, x, rank*sizeof(double));
 | 
			
		||||
  for(int j = 0; j < rank; j++)
 | 
			
		||||
    y[i] -= schlaefliMatrix[i][j]*x[j];
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
static int check_equal_sector(gsl_matrix *x1, gsl_matrix *x2, gsl_matrix *yinv, gsl_matrix *schlaefli, gsl_vector *starting, gsl_vector *i1, gsl_vector *i2)
 | 
			
		||||
{
 | 
			
		||||
  double scalarProduct;
 | 
			
		||||
  int rank = x1->size1;
 | 
			
		||||
 | 
			
		||||
  // calculate <y^{-1} * x1 * x2 * s, e_i> for all i, or identically the components of schlaefli*y^{-1}*x*s
 | 
			
		||||
  gsl_blas_dgemv(CblasNoTrans, 1.0, x2, starting, 0.0, i1);
 | 
			
		||||
  gsl_blas_dgemv(CblasNoTrans, 1.0, x1, i1, 0.0, i2);
 | 
			
		||||
  gsl_blas_dgemv(CblasNoTrans, 1.0, yinv, i2, 0.0, i1);
 | 
			
		||||
  gsl_blas_dgemv(CblasNoTrans, 1.0, schlaefli, i1, 0.0, i2);
 | 
			
		||||
  for(int i = 0; i < rank; i++)
 | 
			
		||||
    if(gsl_vector_get(i2, i) < 0)
 | 
			
		||||
      return 0;
 | 
			
		||||
  return 1;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
static void invert(gsl_matrix *x, gsl_matrix *xinv)
 | 
			
		||||
{
 | 
			
		||||
  int rank = x->size1;
 | 
			
		||||
  gsl_matrix *xLU = gsl_matrix_alloc(rank, rank);
 | 
			
		||||
  gsl_permutation *p = gsl_permutation_alloc(rank);
 | 
			
		||||
  int s;
 | 
			
		||||
  gsl_matrix_memcpy(xLU, x);
 | 
			
		||||
  gsl_linalg_LU_decomp(xLU, p, &s);
 | 
			
		||||
  gsl_linalg_LU_invert(xLU, p, xinv);
 | 
			
		||||
  gsl_permutation_free(p);
 | 
			
		||||
  gsl_matrix_free(xLU);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
static void generate_schlaefli_matrix(semisimple_type_t type, gsl_matrix *mat)
 | 
			
		||||
{
 | 
			
		||||
  gsl_matrix_set_zero(mat);
 | 
			
		||||
  int offset = 0;
 | 
			
		||||
 | 
			
		||||
  for(int k = 0; k < type.n; k++) {
 | 
			
		||||
    if(type.factors[k].series == 'A') {
 | 
			
		||||
      for(int i = 0; i < type.factors[k].rank; i++)
 | 
			
		||||
	for(int j = 0; j < type.factors[k].rank; j++)
 | 
			
		||||
	  if(i == j)
 | 
			
		||||
	    gsl_matrix_set(mat, offset + i, offset + j, 2.0);
 | 
			
		||||
	  else if(i - j == 1 || i - j == -1)
 | 
			
		||||
	    gsl_matrix_set(mat, offset + i, offset + j, -1.0);
 | 
			
		||||
    } else if(type.factors[k].series == 'B' || type.factors[k].series == 'C') {
 | 
			
		||||
      for(int i = 0; i < type.factors[k].rank; i++)
 | 
			
		||||
	for(int j = 0; j < type.factors[k].rank; j++)
 | 
			
		||||
	  if(i == j)
 | 
			
		||||
	    gsl_matrix_set(mat, offset + i, offset + j, 2.0);
 | 
			
		||||
	  else if(i == 0 && j == 1 || i == 1 && j == 0)
 | 
			
		||||
	    gsl_matrix_set(mat, offset + i, offset + j, -sqrt(2.0));
 | 
			
		||||
	  else if(i - j == 1 || i - j == -1)
 | 
			
		||||
	    gsl_matrix_set(mat, offset + i, offset + j, -1.0);
 | 
			
		||||
    } else if(type.factors[k].series == 'D') {
 | 
			
		||||
      for(int i = 0; i < type.factors[k].rank; i++)
 | 
			
		||||
	for(int j = 0; j < type.factors[k].rank; j++)
 | 
			
		||||
	  if(i == j)
 | 
			
		||||
	    gsl_matrix_set(mat, offset + i, offset + j, 2.0);
 | 
			
		||||
	  else if(i > 0 && j > 0 && (i - j == 1 || i - j == -1) || i == 0 && j == 2 || i == 2 && j == 0)
 | 
			
		||||
	    gsl_matrix_set(mat, offset + i, offset + j, -1.0);
 | 
			
		||||
    } else if(type.factors[k].series == 'F') {
 | 
			
		||||
      ERROR(type.factors[k].rank != 4, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
 | 
			
		||||
      for(int i = 0; i < 4; i++)
 | 
			
		||||
	gsl_matrix_set(mat, offset + i, offset + i, 2.0);
 | 
			
		||||
      for(int i = 0; i < 3; i++) {
 | 
			
		||||
	gsl_matrix_set(mat, offset + i, offset + i + 1, i == 1 ? -sqrt(2.0) : -1.0);
 | 
			
		||||
	gsl_matrix_set(mat, offset + i + 1, offset + i, i == 1 ? -sqrt(2.0) : -1.0);
 | 
			
		||||
      }
 | 
			
		||||
    } else if(type.factors[k].series == 'G') {
 | 
			
		||||
      ERROR(type.factors[k].rank != 2, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
 | 
			
		||||
      gsl_matrix_set(mat, offset + 0, offset + 0, 2.0);
 | 
			
		||||
      gsl_matrix_set(mat, offset + 0, offset + 1, -sqrt(3.0));
 | 
			
		||||
      gsl_matrix_set(mat, offset + 1, offset + 0, -sqrt(3.0));
 | 
			
		||||
      gsl_matrix_set(mat, offset + 1, offset + 1, 2.0);
 | 
			
		||||
    } else {
 | 
			
		||||
      ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    offset += type.factors[k].rank;
 | 
			
		||||
  }
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
int coxeter_rank(semisimple_type_t type)
 | 
			
		||||
{
 | 
			
		||||
  int rank = 0;
 | 
			
		||||
  for(int i = 0; i < type.n; i++)
 | 
			
		||||
    rank += type.factors[i].rank;
 | 
			
		||||
  return rank;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
int coxeter_order(semisimple_type_t type)
 | 
			
		||||
{
 | 
			
		||||
  int order = 1;
 | 
			
		||||
  for(int i = 0; i < type.n; i++) {
 | 
			
		||||
    if(type.factors[i].series == 'A') { // (rank+1)!
 | 
			
		||||
      for(int j = 1; j <= type.factors[i].rank + 1; j++)
 | 
			
		||||
	order *= j;
 | 
			
		||||
    } else if(type.factors[i].series == 'B' || type.factors[i].series == 'C') { // 2^rank * rank!
 | 
			
		||||
      for(int j = 1; j <= type.factors[i].rank; j++)
 | 
			
		||||
	order *= 2*j;
 | 
			
		||||
    } else if(type.factors[i].series == 'D') { // 2^(rank-1) * rank!
 | 
			
		||||
      for(int j = 2; j <= type.factors[i].rank; j++)
 | 
			
		||||
	order *= 2*j;
 | 
			
		||||
    } else if(type.factors[i].series == 'E') {
 | 
			
		||||
      if(type.factors[i].rank == 6)
 | 
			
		||||
	order *= 51840;
 | 
			
		||||
      else if(type.factors[i].rank == 7)
 | 
			
		||||
	order *= 2903040;
 | 
			
		||||
      else if(type.factors[i].rank == 8)
 | 
			
		||||
	order *= 696729600;
 | 
			
		||||
      else
 | 
			
		||||
	ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
    } else if(type.factors[i].series == 'F') {
 | 
			
		||||
      ERROR(type.factors[i].rank != 4, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
      order *= 1152;
 | 
			
		||||
    } else if(type.factors[i].series == 'G') {
 | 
			
		||||
      ERROR(type.factors[i].rank != 2, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
      order *= 12;
 | 
			
		||||
    } else if(type.factors[i].series == 'H') {
 | 
			
		||||
      if(type.factors[i].rank == 2)
 | 
			
		||||
	order *= 10;
 | 
			
		||||
      else if(type.factors[i].rank == 3)
 | 
			
		||||
	order *= 120;
 | 
			
		||||
      else if(type.factors[i].rank == 4)
 | 
			
		||||
	order *= 14400;
 | 
			
		||||
      else
 | 
			
		||||
	ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
    } else {
 | 
			
		||||
      ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
    }
 | 
			
		||||
  }
 | 
			
		||||
  return order;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
int coxeter_hyperplanes(semisimple_type_t type)
 | 
			
		||||
{
 | 
			
		||||
  int hyperplanes = 0;
 | 
			
		||||
  for(int i = 0; i < type.n; i++) {
 | 
			
		||||
    if(type.factors[i].series == 'A') // rank*(rank+1)/2
 | 
			
		||||
      hyperplanes += (type.factors[i].rank * (type.factors[i].rank + 1)) / 2;
 | 
			
		||||
    else if(type.factors[i].series == 'B' || type.factors[i].series == 'C') // rank * rank
 | 
			
		||||
      hyperplanes += type.factors[i].rank * type.factors[i].rank;
 | 
			
		||||
    else if(type.factors[i].series == 'D') // rank * (rank - 1)
 | 
			
		||||
      hyperplanes += type.factors[i].rank * (type.factors[i].rank - 1);
 | 
			
		||||
    else if(type.factors[i].series == 'E') {
 | 
			
		||||
      if(type.factors[i].rank == 6)
 | 
			
		||||
	hyperplanes += 36;
 | 
			
		||||
      else if(type.factors[i].rank == 7)
 | 
			
		||||
	hyperplanes += 63;
 | 
			
		||||
      else if(type.factors[i].rank == 8)
 | 
			
		||||
	hyperplanes += 120;
 | 
			
		||||
      else
 | 
			
		||||
	ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
    } else if(type.factors[i].series == 'F') {
 | 
			
		||||
      ERROR(type.factors[i].rank != 4, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
      hyperplanes += 24;
 | 
			
		||||
    } else if(type.factors[i].series == 'G') {
 | 
			
		||||
      ERROR(type.factors[i].rank != 2, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
      hyperplanes += 6;
 | 
			
		||||
    } else {
 | 
			
		||||
      ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
    }
 | 
			
		||||
  }
 | 
			
		||||
  return hyperplanes;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
unsigned long opposition_involution(semisimple_type_t type, unsigned long theta)
 | 
			
		||||
{
 | 
			
		||||
  int offset = 0;
 | 
			
		||||
  unsigned long result = 0;
 | 
			
		||||
  for(int i = 0; i < type.n; i++) {
 | 
			
		||||
    unsigned long current = (theta >> offset) & ((1 << type.factors[i].rank) - 1);
 | 
			
		||||
    unsigned long iota_current;
 | 
			
		||||
 | 
			
		||||
    if(type.factors[i].series == 'B' || type.factors[i].series == 'C' || type.factors[i].series == 'F' || type.factors[i].series == 'G') {
 | 
			
		||||
      iota_current = current;
 | 
			
		||||
    } else if(type.factors[i].series == 'A') {
 | 
			
		||||
      iota_current = 0;
 | 
			
		||||
      for(int j = 0; j < type.factors[i].rank; j++)
 | 
			
		||||
	iota_current += ((current >> j) & 1) << (type.factors[i].rank - 1 - j);
 | 
			
		||||
    } else if(type.factors[i].series == 'D') {
 | 
			
		||||
      if(type.factors[i].rank % 2 == 0) {
 | 
			
		||||
	iota_current = current;
 | 
			
		||||
      } else {
 | 
			
		||||
	ERROR(1, "The opposition involution for type %c%d is not yet implemented!\n", type.factors[i].series, type.factors[i].rank);
 | 
			
		||||
      }
 | 
			
		||||
    } else if(type.factors[i].series == 'E') {
 | 
			
		||||
      ERROR(1, "The opposition involution for En is not yet implemented!\n");
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    result += iota_current << offset;
 | 
			
		||||
    offset += type.factors[i].rank;
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  return result;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
static void generate_starting_vector(int rank, gsl_matrix *schlaefli, gsl_vector *result)
 | 
			
		||||
{
 | 
			
		||||
  gsl_matrix *schlaefliLU = gsl_matrix_alloc(rank, rank);
 | 
			
		||||
  gsl_vector *diagonal = gsl_vector_alloc(rank);
 | 
			
		||||
  gsl_permutation *p = gsl_permutation_alloc(rank);
 | 
			
		||||
  int s;
 | 
			
		||||
  for(int i = 0; i < rank; i++)
 | 
			
		||||
    gsl_vector_set(diagonal, i, 1.0);
 | 
			
		||||
  gsl_matrix_memcpy(schlaefliLU, schlaefli);
 | 
			
		||||
  gsl_linalg_LU_decomp(schlaefliLU, p, &s);
 | 
			
		||||
  gsl_linalg_LU_solve(schlaefliLU, p, diagonal, result);
 | 
			
		||||
  gsl_permutation_free(p);
 | 
			
		||||
  gsl_vector_free(diagonal);
 | 
			
		||||
  gsl_matrix_free(schlaefliLU);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
void generate_coxeter_graph(semisimple_type_t type, int *result)
 | 
			
		||||
{
 | 
			
		||||
  int rank = coxeter_rank(type);
 | 
			
		||||
  int order = coxeter_order(type);
 | 
			
		||||
  int element_count;
 | 
			
		||||
  gsl_matrix *schlaefliMatrix = gsl_matrix_alloc(rank, rank);
 | 
			
		||||
  gsl_vector *startingVector = gsl_vector_alloc(rank);
 | 
			
		||||
  gsl_matrix *generators = (gsl_matrix*)malloc(rank*sizeof(gsl_matrix));
 | 
			
		||||
  double *generators_data = (double*)malloc(rank*rank*rank*sizeof(double));
 | 
			
		||||
  gsl_matrix *elements = (gsl_matrix*)malloc(order*sizeof(gsl_matrix));
 | 
			
		||||
  double *elements_data = (double*)malloc(rank*rank*order*sizeof(double));
 | 
			
		||||
  gsl_matrix *inverses = (gsl_matrix*)malloc(order*sizeof(gsl_matrix));
 | 
			
		||||
  double *inverses_data = (double*)malloc(rank*rank*order*sizeof(double));
 | 
			
		||||
  gsl_vector *i1 = gsl_vector_alloc(rank);
 | 
			
		||||
  gsl_vector *i2 = gsl_vector_alloc(rank);
 | 
			
		||||
  gsl_matrix *current_element = gsl_matrix_alloc(rank, rank);
 | 
			
		||||
  queue_t queue;
 | 
			
		||||
  int current;
 | 
			
		||||
 | 
			
		||||
  for(int i = 0; i < rank; i++)
 | 
			
		||||
    generators[i] = gsl_matrix_view_array(generators_data + i*rank*rank, rank, rank).matrix;
 | 
			
		||||
  for(int i = 0; i < order; i++)
 | 
			
		||||
    elements[i] = gsl_matrix_view_array(elements_data + i*rank*rank, rank, rank).matrix;
 | 
			
		||||
  for(int i = 0; i < order; i++)
 | 
			
		||||
    inverses[i] = gsl_matrix_view_array(inverses_data + i*rank*rank, rank, rank).matrix;
 | 
			
		||||
 | 
			
		||||
  generate_schlaefli_matrix(type, schlaefliMatrix);
 | 
			
		||||
  generate_starting_vector(rank, schlaefliMatrix, startingVector);
 | 
			
		||||
 | 
			
		||||
  for(int i = 0; i < rank; i++) {
 | 
			
		||||
    gsl_matrix_set_identity(&generators[i]);
 | 
			
		||||
    for(int j = 0; j < rank; j++)
 | 
			
		||||
      if(i == j)
 | 
			
		||||
	gsl_matrix_set(&generators[i], i, j, -1.0);
 | 
			
		||||
      else
 | 
			
		||||
	gsl_matrix_set(&generators[i], i, j, -gsl_matrix_get(schlaefliMatrix, i, j));
 | 
			
		||||
    //    gsl_matrix_fprintf(stdout, &generators[i], "%f");
 | 
			
		||||
    //    printf("\n");
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  queue_init(&queue);
 | 
			
		||||
  queue_put(&queue, 0);
 | 
			
		||||
  element_count = 1;
 | 
			
		||||
  gsl_matrix_set_identity(&elements[0]);
 | 
			
		||||
  gsl_matrix_set_identity(&inverses[0]);
 | 
			
		||||
  while((current = queue_get(&queue)) != -1) {
 | 
			
		||||
    for(int i = 0; i < rank; i++) {
 | 
			
		||||
      int j;
 | 
			
		||||
      for(j = 0; j < element_count; j++) {
 | 
			
		||||
	if(check_equal_sector(&generators[i], &elements[current], &inverses[j], schlaefliMatrix, startingVector, i1, i2)) { // generators[i] * elements[current] = elements[j]
 | 
			
		||||
	  result[rank*current + i] = j;
 | 
			
		||||
	  break;
 | 
			
		||||
	}
 | 
			
		||||
      }
 | 
			
		||||
      // if no existing element equals generators[i] * elements[current], create one
 | 
			
		||||
      if(j == element_count) {
 | 
			
		||||
	ERROR(element_count >= order, "Got more elements than the order of the group should be!\n");
 | 
			
		||||
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &generators[i], &elements[current], 0.0, &elements[element_count]);
 | 
			
		||||
	invert(&elements[element_count], &inverses[element_count]);
 | 
			
		||||
	result[rank*current + i] = element_count;
 | 
			
		||||
	queue_put(&queue, element_count);
 | 
			
		||||
	element_count++;
 | 
			
		||||
      }
 | 
			
		||||
    }
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  ERROR(element_count != order, "Something went wrong building the Coxeter group. Found %d elements, %d expected\n", element_count, order);
 | 
			
		||||
 | 
			
		||||
  /*
 | 
			
		||||
  for(int i = 0; i < order; i++) {
 | 
			
		||||
    printf("%d: ", i);
 | 
			
		||||
    for(int j = 0; j < rank; j++)
 | 
			
		||||
      printf("%d ", result[rank*i + j]);
 | 
			
		||||
    printf("\n");
 | 
			
		||||
  }
 | 
			
		||||
  */
 | 
			
		||||
 | 
			
		||||
  gsl_vector_free(startingVector);
 | 
			
		||||
  gsl_vector_free(i1);
 | 
			
		||||
  gsl_vector_free(i2);
 | 
			
		||||
  gsl_matrix_free(schlaefliMatrix);
 | 
			
		||||
  gsl_matrix_free(current_element);
 | 
			
		||||
  free(generators);
 | 
			
		||||
  free(elements);
 | 
			
		||||
  free(generators_data);
 | 
			
		||||
  free(elements_data);
 | 
			
		||||
}
 | 
			
		||||
@@ -1,20 +0,0 @@
 | 
			
		||||
#ifndef COXETER_H
 | 
			
		||||
#define COXETER_H
 | 
			
		||||
 | 
			
		||||
typedef struct {
 | 
			
		||||
  char series;
 | 
			
		||||
  int rank;
 | 
			
		||||
} simple_type_t;
 | 
			
		||||
 | 
			
		||||
typedef struct {
 | 
			
		||||
  int n;
 | 
			
		||||
  simple_type_t *factors;
 | 
			
		||||
} semisimple_type_t;
 | 
			
		||||
 | 
			
		||||
void generate_coxeter_graph(semisimple_type_t type, int *result);
 | 
			
		||||
int coxeter_order(semisimple_type_t type);
 | 
			
		||||
int coxeter_hyperplanes(semisimple_type_t type);
 | 
			
		||||
int coxeter_rank(semisimple_type_t type);
 | 
			
		||||
unsigned long opposition_involution(semisimple_type_t type, unsigned long theta);
 | 
			
		||||
 | 
			
		||||
#endif
 | 
			
		||||
@@ -1,193 +0,0 @@
 | 
			
		||||
#include <stdio.h>
 | 
			
		||||
#include <string.h>
 | 
			
		||||
#include <sys/stat.h>
 | 
			
		||||
 | 
			
		||||
#include "thickenings.h"
 | 
			
		||||
#include "coxeter.h"
 | 
			
		||||
#include "queue.h"
 | 
			
		||||
 | 
			
		||||
#define SWAP(t, a, b) do {t tmp = a; a = b; b = tmp;} while(0)
 | 
			
		||||
 | 
			
		||||
char *stringify_SLn1_permutation(int *word, int wordlength, int rank, char *str)
 | 
			
		||||
{
 | 
			
		||||
  for(int i = 0; i <= rank; i++)
 | 
			
		||||
    str[i] = '1' + i;
 | 
			
		||||
  str[rank+1] = 0;
 | 
			
		||||
 | 
			
		||||
  for(int i = 0; i < wordlength; i++) {
 | 
			
		||||
    char tmp = str[word[i]];
 | 
			
		||||
    str[word[i]] = str[word[i]+1];
 | 
			
		||||
    str[word[i]+1] = tmp;
 | 
			
		||||
  }
 | 
			
		||||
  return str;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
char *stringify_Onn1_permutation(int *word, int wordlength, int rank, char *str)
 | 
			
		||||
{
 | 
			
		||||
  for(int i = 0; i <= rank*2; i++)
 | 
			
		||||
    str[i] = '1' + i;
 | 
			
		||||
  str[2*rank+1] = 0;
 | 
			
		||||
 | 
			
		||||
  for(int i = 0; i < wordlength; i++) {
 | 
			
		||||
    if(word[i] == 0)
 | 
			
		||||
      SWAP(char, str[rank-1], str[rank+1]);
 | 
			
		||||
    else {
 | 
			
		||||
      SWAP(char, str[rank-word[i]], str[rank-word[i]-1]);
 | 
			
		||||
      SWAP(char, str[rank+word[i]], str[rank+word[i]+1]);
 | 
			
		||||
    }
 | 
			
		||||
  }
 | 
			
		||||
  return str;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
int main(int argc, const char *argv[])
 | 
			
		||||
{
 | 
			
		||||
  FILE *infile;
 | 
			
		||||
  struct stat st;
 | 
			
		||||
  int rank, order, hyperplanes;
 | 
			
		||||
  semisimple_type_t type;
 | 
			
		||||
  int n;
 | 
			
		||||
  signed char *level;
 | 
			
		||||
  node_t *graph;
 | 
			
		||||
  int *left, *right;
 | 
			
		||||
  int left_invariant, right_invariant;
 | 
			
		||||
  int left_invariant_wanted = 0, right_invariant_wanted = 0;
 | 
			
		||||
 | 
			
		||||
  unsigned long left_invariance, right_invariance;
 | 
			
		||||
  edgelist_t *edgelists;
 | 
			
		||||
  int *words;
 | 
			
		||||
 | 
			
		||||
  queue_t queue;
 | 
			
		||||
  int current;
 | 
			
		||||
  int *seen;
 | 
			
		||||
  int *generators;
 | 
			
		||||
  int ngens;
 | 
			
		||||
 | 
			
		||||
  char string_buffer1[1000];
 | 
			
		||||
  const char *alphabet = "abcdefghijklmnopqrstuvwxyz";
 | 
			
		||||
 | 
			
		||||
  // parse arguments
 | 
			
		||||
 | 
			
		||||
  if(argc < 2)
 | 
			
		||||
    infile = stdin;
 | 
			
		||||
  else {
 | 
			
		||||
    if(strcmp(argv[1], "-") == 0)
 | 
			
		||||
      infile = stdin;
 | 
			
		||||
    else
 | 
			
		||||
      infile = fopen(argv[1], "rb");
 | 
			
		||||
 | 
			
		||||
    if(argc >= 4) {
 | 
			
		||||
      if(strcmp(argv[2], "-") != 0)
 | 
			
		||||
	for(int i = 0; i < strlen(argv[2]); i++)
 | 
			
		||||
	  left_invariant_wanted |= (1 << (argv[2][i] - 'a'));
 | 
			
		||||
      if(strcmp(argv[3], "-") != 0)
 | 
			
		||||
	for(int i = 0; i < strlen(argv[3]); i++)
 | 
			
		||||
	  right_invariant_wanted |= (1 << (argv[3][i] - 'a'));
 | 
			
		||||
    }
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  fread(&type.n, sizeof(int), 1, infile);  // we completely trust the input data
 | 
			
		||||
  type.factors = malloc(type.n * sizeof(simple_type_t));
 | 
			
		||||
  fread(type.factors, sizeof(simple_type_t), type.n, infile);
 | 
			
		||||
  fread(&left_invariance, sizeof(simple_type_t), type.n, infile);
 | 
			
		||||
  fread(&right_invariance, sizeof(simple_type_t), type.n, infile);
 | 
			
		||||
 | 
			
		||||
  // get graph
 | 
			
		||||
 | 
			
		||||
  rank = coxeter_rank(type);
 | 
			
		||||
  order = coxeter_order(type);
 | 
			
		||||
  hyperplanes = coxeter_hyperplanes(type);
 | 
			
		||||
  ERROR(strlen(alphabet) < rank, "The alphabet has too few letters\n");
 | 
			
		||||
  seen = (int*)malloc(order*sizeof(int));
 | 
			
		||||
  generators = (int*)malloc(order*sizeof(int));
 | 
			
		||||
  level = (signed char*)malloc(order*sizeof(int));
 | 
			
		||||
 | 
			
		||||
  graph = graph_alloc(type);
 | 
			
		||||
  prepare_graph(type, graph);
 | 
			
		||||
 | 
			
		||||
  // finally do stuff
 | 
			
		||||
 | 
			
		||||
  int counter = 0;
 | 
			
		||||
 | 
			
		||||
  while(fread(level, sizeof(signed char), order, infile) == order) {
 | 
			
		||||
 | 
			
		||||
    /*
 | 
			
		||||
    if((counter++) % 100000 == 0)
 | 
			
		||||
      print_thickening(rank, order, level, 0, 0, 0, alphabet, stdout);
 | 
			
		||||
    continue;
 | 
			
		||||
    */
 | 
			
		||||
 | 
			
		||||
    left_invariant = right_invariant = -1; // all 1s
 | 
			
		||||
    for(int j = 0; j < order; j++) {
 | 
			
		||||
      for(int k = 0; k < rank; k++) {
 | 
			
		||||
	if(level[j] > 0 && level[graph[j].left[k]] < 0 || level[j] < 0 && level[graph[j].left[k]] > 0) {
 | 
			
		||||
	  left_invariant &= ~(1 << k);
 | 
			
		||||
	}
 | 
			
		||||
	if(level[j] > 0 && level[graph[j].right[k]] < 0 || level[j] < 0 && level[graph[j].right[k]] > 0) {
 | 
			
		||||
	  right_invariant &= ~(1 << k);
 | 
			
		||||
	}
 | 
			
		||||
      }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    if((~left_invariant & left_invariant_wanted) == 0 && (~right_invariant & right_invariant_wanted) == 0) {
 | 
			
		||||
      ngens = 0;
 | 
			
		||||
      memset(generators, 0, order*sizeof(int));
 | 
			
		||||
      for(int j = 0; j < order; j++) {
 | 
			
		||||
	if(level[j] == HEAD_MARKER && generators[j] == 0) { // ignore the generator, if it is equivalent to one already seen
 | 
			
		||||
	  ngens++;
 | 
			
		||||
	  queue_init(&queue);
 | 
			
		||||
	  queue_put(&queue, j);
 | 
			
		||||
	  while((current = queue_get(&queue)) != -1) {
 | 
			
		||||
	    if(generators[current] == 0) { // visit everyone only once
 | 
			
		||||
	      generators[current] = ngens;
 | 
			
		||||
	      for(int k = 0; k < rank; k++) {
 | 
			
		||||
		if(left_invariant & (1 << k))
 | 
			
		||||
		queue_put(&queue, graph[current].left[k]);
 | 
			
		||||
		if(right_invariant & (1 << k))
 | 
			
		||||
		  queue_put(&queue, graph[current].right[k]);
 | 
			
		||||
	      }
 | 
			
		||||
	    }
 | 
			
		||||
	  }
 | 
			
		||||
	}
 | 
			
		||||
      }
 | 
			
		||||
 | 
			
		||||
      printf("left: ");
 | 
			
		||||
      for(int j = 0; j < rank; j++)
 | 
			
		||||
	printf("%c", left_invariant & (1 << j) ? alphabet[j] : ' ');
 | 
			
		||||
      printf(" right: ");
 | 
			
		||||
      for(int j = 0; j < rank; j++)
 | 
			
		||||
	printf("%c", right_invariant & (1 << j) ? alphabet[j] : ' ');
 | 
			
		||||
      printf(" generators: ");
 | 
			
		||||
 | 
			
		||||
      memset(seen, 0, order*sizeof(int));
 | 
			
		||||
      for(int i = 0; i < order; i++) {
 | 
			
		||||
	if(generators[i] != 0 && seen[generators[i]-1] == 0) {
 | 
			
		||||
	  seen[generators[i]-1] = 1;
 | 
			
		||||
	  //	if(type.n == 1 && type.factors[0].series == 'A')
 | 
			
		||||
	  //	  printf("%s ", stringify_SLn1_permutation(graph[i].word, graph[i].wordlength, rank, string_buffer1));
 | 
			
		||||
	  //	else if(type.n == 1 && (type.factors[0].series == 'B' || type.factors[0].series == 'C'))
 | 
			
		||||
	  //	  printf("%s ", stringify_Onn1_permutation(graph[i].word, graph[i].wordlength, rank, string_buffer1));
 | 
			
		||||
	  //	else
 | 
			
		||||
	  if(i == 0)
 | 
			
		||||
	    printf("1 ");
 | 
			
		||||
	  else
 | 
			
		||||
	    printf("%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1));
 | 
			
		||||
	}
 | 
			
		||||
      }
 | 
			
		||||
 | 
			
		||||
      printf("\n");
 | 
			
		||||
    }
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  if(infile != stdin)
 | 
			
		||||
    fclose(infile);
 | 
			
		||||
 | 
			
		||||
  // cleanup
 | 
			
		||||
 | 
			
		||||
  graph_free(type, graph);
 | 
			
		||||
  free(seen);
 | 
			
		||||
  free(generators);
 | 
			
		||||
  free(type.factors);
 | 
			
		||||
 | 
			
		||||
  return 0;
 | 
			
		||||
}
 | 
			
		||||
							
								
								
									
										189
									
								
								old/process.c
									
									
									
									
									
								
							
							
						
						
									
										189
									
								
								old/process.c
									
									
									
									
									
								
							@@ -1,189 +0,0 @@
 | 
			
		||||
#include <stdio.h>
 | 
			
		||||
#include <string.h>
 | 
			
		||||
#include <sys/stat.h>
 | 
			
		||||
 | 
			
		||||
#include "thickenings.h"
 | 
			
		||||
#include "weyl.h"
 | 
			
		||||
#include "queue.h"
 | 
			
		||||
 | 
			
		||||
int main(int argc, const char *argv[])
 | 
			
		||||
{
 | 
			
		||||
  FILE *infile;
 | 
			
		||||
  semisimple_type_t type;
 | 
			
		||||
  unsigned long left_invariance, right_invariance; // these are the invariances we have already modded out
 | 
			
		||||
  unsigned long left_invariant, right_invariant; // these are the invariances of the thickening under consideration
 | 
			
		||||
  int rank, cosets;
 | 
			
		||||
  node_t *graph;
 | 
			
		||||
  signed char *thickening;
 | 
			
		||||
  int *seen, *generators;
 | 
			
		||||
  queue_t queue;
 | 
			
		||||
  int ngenerators;
 | 
			
		||||
  int current;
 | 
			
		||||
 | 
			
		||||
  char string_buffer1[1000];
 | 
			
		||||
  const char *alphabet = "abcdefghijklmnopqrstuvwxyz";
 | 
			
		||||
 | 
			
		||||
  if(argc < 2)
 | 
			
		||||
    infile = stdin;
 | 
			
		||||
  else
 | 
			
		||||
    infile = fopen(argv[1], "rb");
 | 
			
		||||
 | 
			
		||||
  // we completely trust the input data
 | 
			
		||||
  ERROR(fread(&type.n, sizeof(int), 1, infile) == 0, "The input file seems to be empty!\n");
 | 
			
		||||
  type.factors = malloc(type.n * sizeof(simple_type_t));
 | 
			
		||||
  fread(type.factors, sizeof(simple_type_t), type.n, infile);
 | 
			
		||||
  fread(&left_invariance, sizeof(simple_type_t), type.n, infile);
 | 
			
		||||
  fread(&right_invariance, sizeof(simple_type_t), type.n, infile);
 | 
			
		||||
 | 
			
		||||
  rank = weyl_rank(type);
 | 
			
		||||
  graph = graph_alloc(type);
 | 
			
		||||
  cosets = prepare_simplified_graph(type, left_invariance, right_invariance, graph);
 | 
			
		||||
 | 
			
		||||
  thickening = (signed char*)malloc(cosets*sizeof(signed char));
 | 
			
		||||
  generators = (int*)malloc(cosets*sizeof(int));
 | 
			
		||||
  seen = (int*)malloc(cosets*sizeof(int));
 | 
			
		||||
 | 
			
		||||
  while(fread(thickening, sizeof(signed char), cosets, infile) == cosets) {
 | 
			
		||||
 | 
			
		||||
    // determine invariances of this thickening
 | 
			
		||||
    left_invariant = right_invariant = -1; // set all bits to 1
 | 
			
		||||
    for(int j = 0; j < cosets; j++) {
 | 
			
		||||
      for(int k = 0; k < rank; k++) {
 | 
			
		||||
	if(thickening[j] > 0 && thickening[graph[j].left[k]] < 0 ||
 | 
			
		||||
	   thickening[j] < 0 && thickening[graph[j].left[k]] > 0)
 | 
			
		||||
	  left_invariant &= ~(1 << k);
 | 
			
		||||
	if(thickening[j] > 0 && thickening[graph[j].right[k]] < 0 ||
 | 
			
		||||
	   thickening[j] < 0 && thickening[graph[j].right[k]] > 0)
 | 
			
		||||
	  right_invariant &= ~(1 << k);
 | 
			
		||||
      }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    // print this stuff
 | 
			
		||||
    printf("left: ");
 | 
			
		||||
    for(int j = 0; j < rank; j++)
 | 
			
		||||
      printf("%c", left_invariant & (1 << j) ? alphabet[j] : ' ');
 | 
			
		||||
    printf(" right: ");
 | 
			
		||||
    for(int j = 0; j < rank; j++)
 | 
			
		||||
      printf("%c", right_invariant & (1 << j) ? alphabet[j] : ' ');
 | 
			
		||||
    printf(" generators: ");
 | 
			
		||||
 | 
			
		||||
    // find a minimal set of weyl group elements such that the union of the ideals generated by their cosets wrt the invariances determined above gives the thickening
 | 
			
		||||
    // in the first step, mark everything which is equivalent to a "head" by a generator id
 | 
			
		||||
    ngenerators = 0;
 | 
			
		||||
    memset(generators, 0, cosets*sizeof(int));
 | 
			
		||||
    for(int j = 0; j < cosets; j++) {
 | 
			
		||||
      if(thickening[j] == HEAD_MARKER && generators[j] == 0) { // ignore the generator, if it is equivalent to one already seen
 | 
			
		||||
	ngenerators++;
 | 
			
		||||
	queue_init(&queue);
 | 
			
		||||
	queue_put(&queue, j);
 | 
			
		||||
	while((current = queue_get(&queue)) != -1) {
 | 
			
		||||
	  if(generators[current] == 0) { // visit everyone only once
 | 
			
		||||
	    generators[current] = ngenerators;
 | 
			
		||||
	    for(int k = 0; k < rank; k++) {
 | 
			
		||||
	      if(left_invariant & (1 << k))
 | 
			
		||||
		queue_put(&queue, graph[current].left[k]);
 | 
			
		||||
	      if(right_invariant & (1 << k))
 | 
			
		||||
		queue_put(&queue, graph[current].right[k]);
 | 
			
		||||
	    }
 | 
			
		||||
	  }
 | 
			
		||||
	}
 | 
			
		||||
      }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    // in the second step, go through the list in ascending word length order and print the first appearance of each generator id
 | 
			
		||||
    memset(seen, 0, cosets*sizeof(int));
 | 
			
		||||
    for(int i = 0; i < cosets; i++) {
 | 
			
		||||
      if(generators[i] != 0 && seen[generators[i]-1] == 0) {
 | 
			
		||||
	seen[generators[i]-1] = 1;
 | 
			
		||||
	printf("%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1));
 | 
			
		||||
      }
 | 
			
		||||
    }
 | 
			
		||||
 | 
			
		||||
    printf("\n");
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  if(infile != stdin)
 | 
			
		||||
    fclose(infile);
 | 
			
		||||
 | 
			
		||||
  graph_free(type, graph);
 | 
			
		||||
  free(type.factors);
 | 
			
		||||
  free(thickening);
 | 
			
		||||
}
 | 
			
		||||
  /*******************************************************************************************
 | 
			
		||||
 | 
			
		||||
  seen = (int*)malloc(order*sizeof(int));
 | 
			
		||||
  generators = (int*)malloc(order*sizeof(int));
 | 
			
		||||
  level = (signed char*)malloc(order*sizeof(int));
 | 
			
		||||
 | 
			
		||||
  graph = graph_alloc(type);
 | 
			
		||||
  prepare_graph(type, graph);
 | 
			
		||||
 | 
			
		||||
  // finally do stuff
 | 
			
		||||
 | 
			
		||||
  int counter = 0;
 | 
			
		||||
 | 
			
		||||
  while(fread(level, sizeof(signed char), order, infile) == order) {
 | 
			
		||||
 | 
			
		||||
    if((~left_invariant & left_invariant_wanted) == 0 && (~right_invariant & right_invariant_wanted) == 0) {
 | 
			
		||||
      ngens = 0;
 | 
			
		||||
      memset(generators, 0, order*sizeof(int));
 | 
			
		||||
      for(int j = 0; j < order; j++) {
 | 
			
		||||
	if(level[j] == HEAD_MARKER && generators[j] == 0) { // ignore the generator, if it is equivalent to one already seen
 | 
			
		||||
	  ngens++;
 | 
			
		||||
	  queue_init(&queue);
 | 
			
		||||
	  queue_put(&queue, j);
 | 
			
		||||
	  while((current = queue_get(&queue)) != -1) {
 | 
			
		||||
	    if(generators[current] == 0) { // visit everyone only once
 | 
			
		||||
	      generators[current] = ngens;
 | 
			
		||||
	      for(int k = 0; k < rank; k++) {
 | 
			
		||||
		if(left_invariant & (1 << k))
 | 
			
		||||
		queue_put(&queue, graph[current].left[k]);
 | 
			
		||||
		if(right_invariant & (1 << k))
 | 
			
		||||
		  queue_put(&queue, graph[current].right[k]);
 | 
			
		||||
	      }
 | 
			
		||||
	    }
 | 
			
		||||
	  }
 | 
			
		||||
	}
 | 
			
		||||
      }
 | 
			
		||||
 | 
			
		||||
      printf("left: ");
 | 
			
		||||
      for(int j = 0; j < rank; j++)
 | 
			
		||||
	printf("%c", left_invariant & (1 << j) ? alphabet[j] : ' ');
 | 
			
		||||
      printf(" right: ");
 | 
			
		||||
      for(int j = 0; j < rank; j++)
 | 
			
		||||
	printf("%c", right_invariant & (1 << j) ? alphabet[j] : ' ');
 | 
			
		||||
      printf(" generators: ");
 | 
			
		||||
 | 
			
		||||
      memset(seen, 0, order*sizeof(int));
 | 
			
		||||
      for(int i = 0; i < order; i++) {
 | 
			
		||||
	if(generators[i] != 0 && seen[generators[i]-1] == 0) {
 | 
			
		||||
	  seen[generators[i]-1] = 1;
 | 
			
		||||
	  //	if(type.n == 1 && type.factors[0].series == 'A')
 | 
			
		||||
	  //	  printf("%s ", stringify_SLn1_permutation(graph[i].word, graph[i].wordlength, rank, string_buffer1));
 | 
			
		||||
	  //	else if(type.n == 1 && (type.factors[0].series == 'B' || type.factors[0].series == 'C'))
 | 
			
		||||
	  //	  printf("%s ", stringify_Onn1_permutation(graph[i].word, graph[i].wordlength, rank, string_buffer1));
 | 
			
		||||
	  //	else
 | 
			
		||||
	  if(i == 0)
 | 
			
		||||
	    printf("1 ");
 | 
			
		||||
	  else
 | 
			
		||||
	    printf("%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1));
 | 
			
		||||
	}
 | 
			
		||||
      }
 | 
			
		||||
 | 
			
		||||
      printf("\n");
 | 
			
		||||
    }
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  if(infile != stdin)
 | 
			
		||||
    fclose(infile);
 | 
			
		||||
 | 
			
		||||
  // cleanup
 | 
			
		||||
 | 
			
		||||
  graph_free(type, graph);
 | 
			
		||||
  free(seen);
 | 
			
		||||
  free(generators);
 | 
			
		||||
  free(type.factors);
 | 
			
		||||
 | 
			
		||||
  return 0;
 | 
			
		||||
}
 | 
			
		||||
*/
 | 
			
		||||
							
								
								
									
										41
									
								
								old/test.c
									
									
									
									
									
								
							
							
						
						
									
										41
									
								
								old/test.c
									
									
									
									
									
								
							@@ -1,41 +0,0 @@
 | 
			
		||||
#include <stdio.h>
 | 
			
		||||
 | 
			
		||||
#include "weyl.h"
 | 
			
		||||
 | 
			
		||||
static void print_element(weylgroup_element_t *e)
 | 
			
		||||
{
 | 
			
		||||
  if(e->wordlength == 0)
 | 
			
		||||
    printf("1");
 | 
			
		||||
  else
 | 
			
		||||
    for(int j = 0; j < e->wordlength; j++)
 | 
			
		||||
      printf("%c", e->word[j] + 'a');
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
int main()
 | 
			
		||||
{
 | 
			
		||||
  semisimple_type_t type;
 | 
			
		||||
  simple_type_t t;
 | 
			
		||||
  type.n = 1;
 | 
			
		||||
  type.factors = &t;
 | 
			
		||||
  t.series = 'A';
 | 
			
		||||
  t.rank = 3;
 | 
			
		||||
 | 
			
		||||
  int order = weyl_order(type);
 | 
			
		||||
  doublequotient_t *dq = weyl_generate_bruhat(type, 0x02, 0x03);
 | 
			
		||||
 | 
			
		||||
  for(int i = 0; i < dq->count; i++) {
 | 
			
		||||
    print_element(dq->cosets[i].min);
 | 
			
		||||
    printf(" -> ");
 | 
			
		||||
    for(doublecoset_list_t *current = dq->cosets[i].bruhat_lower; current; current = current->next) {
 | 
			
		||||
      print_element(current->to->min);
 | 
			
		||||
      printf(" ");
 | 
			
		||||
    }
 | 
			
		||||
    printf("| ");
 | 
			
		||||
    print_element(dq->cosets[i].opposite->min);
 | 
			
		||||
    printf("\n");
 | 
			
		||||
  }
 | 
			
		||||
 | 
			
		||||
  weyl_destroy_bruhat(dq);
 | 
			
		||||
 | 
			
		||||
  return 0;
 | 
			
		||||
}
 | 
			
		||||
		Reference in New Issue
	
	Block a user