#include "weyl.h"
#include "queue.h"

#include <stdio.h>
#include <memory.h>
#include <stdlib.h>

#define BIT(n) ((uint64_t)1 << (n))

typedef struct {
  weylid_t id;
  int position;
} weylid_lookup_t;

static int search(const void *key, const void *base, size_t nmem, size_t size, int (*compar) (const void *, const void *, void *), void *arg);
static int compare_root_vectors(int rank, const int *x, const int *y);
static int compare_root_vectors_qsort(const void *x, const void *y, void *arg);
static int compare_weylid_lookup(const void *x, const void *y);
static int lookup_id(weylid_t id, weylid_lookup_t *list, int len);
static weylid_t multiply_generator(int s, weylid_t w, const int *simple, const int *mapping, int rank, int positive);
static void reflect_root_vector(const int *cartan, int rank, int i, int *old, int *new);

/***************** simple helper functions **********************************/

// glibc search function, but with user pointer and returning index (or -1 if not found)
static int search (const void *key, const void *base, size_t nmemb, size_t size, int (*compar) (const void *, const void *, void *), void *arg)
{
  size_t l, u, idx;
  const void *p;
  int comparison;

  l = 0;
  u = nmemb;
  while (l < u) {
    idx = (l + u) / 2;
    p = (void *) (((const char *) base) + (idx * size));
    comparison = (*compar) (key, p, arg);
    if (comparison < 0)
      u = idx;
    else if (comparison > 0)
      l = idx + 1;
    else
      return idx;
  }

  return -1;
}

// maybe we want a different ordering here?
static int compare_root_vectors(int rank, const int *x, const int *y)
{
  for(int i = 0; i < rank; i++)
    if(x[i] != y[i])
      return x[i] - y[i];

  return 0;
}

static int compare_root_vectors_qsort(const void *x, const void *y, void *arg)
{
  return compare_root_vectors(*((int*)arg), x, y);
}

static int compare_weylid(const void *x, const void *y)
{
  weylid_t u = *((weylid_t*)x);
  weylid_t v = *((weylid_t*)y);

  return u > v ? 1 : u < v ? -1 : 0;
}

static int compare_weylid_lookup(const void *x, const void *y)
{
  weylid_t u = ((weylid_lookup_t*)x)->id;
  weylid_t v = ((weylid_lookup_t*)y)->id;

  return u > v ? 1 : u < v ? -1 : 0;
}

static int lookup_id(weylid_t id, weylid_lookup_t *list, int len)
{
  weylid_lookup_t key;
  key.id = id;
  weylid_lookup_t *p = (weylid_lookup_t*)bsearch(&key, list, len, sizeof(weylid_lookup_t), compare_weylid_lookup);
  return p->position;
}

static weylid_t multiply_generator(int s, weylid_t w, const int* simple, const int* mapping, int rank, int positive)
{
  weylid_t sw = 0;

  for(int i = 0; i < positive; i++) {
    if(w & BIT(i))
      if(mapping[i*rank+s] != -1)
	sw |= BIT(mapping[i*rank+s]);
  }

  if(w & BIT(simple[s]))
    return sw;
  else
    return sw | BIT(simple[s]);
}

static void reflect_root_vector(const int *cartan, int rank, int i, int *old, int *new)
{
  memcpy(new, old, rank*sizeof(int));
  for(int j = 0; j < rank; j++)
    new[i] -= cartan[i*rank + j]*old[j];
}

/************* Weyl group infos ************************/

int weyl_rank(semisimple_type_t type)
{
  int rank = 0;
  for(int i = 0; i < type.n; i++)
    rank += type.factors[i].rank;
  return rank;
}

int weyl_order(semisimple_type_t type)
{
  int order = 1;
  for(int i = 0; i < type.n; i++) {
    switch(type.factors[i].series) {

    case 'A':
      for(int j = 1; j <= type.factors[i].rank + 1; j++)
	order *= j;
      break;

    case 'B': case 'C':
      for(int j = 1; j <= type.factors[i].rank; j++)
	order *= 2*j;
      break;

    case 'D':
      for(int j = 2; j <= type.factors[i].rank; j++)
	order *= 2*j;
      break;

    case '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 Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
      break;

    case 'F':
      ERROR(type.factors[i].rank != 4, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
      order *= 1152;
      break;

    case 'G':
      ERROR(type.factors[i].rank != 2, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
      order *= 12;
      break;

    default:
      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
    }
  }

  return order;
}

int weyl_hyperplanes(semisimple_type_t type)
{
  int hyperplanes = 0;

  for(int i = 0; i < type.n; i++) {
    switch(type.factors[i].series) {
    case 'A':
      hyperplanes += (type.factors[i].rank * (type.factors[i].rank + 1)) / 2;
      break;

    case 'B': case 'C':
      hyperplanes += type.factors[i].rank * type.factors[i].rank;
      break;

    case 'D':
      hyperplanes += type.factors[i].rank * (type.factors[i].rank - 1);
      break;

    case '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 Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
      break;

    case 'F':
      ERROR(type.factors[i].rank != 4, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
      hyperplanes += 24;
      break;

    case 'G':
      ERROR(type.factors[i].rank != 2, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
      hyperplanes += 6;
      break;

    default:
      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
    }
  }

  return hyperplanes;
}

int weyl_opposition(semisimple_type_t type, int simple_root)
{
  int offset = 0;
  int factor = 0;
  int r, iota_r;

  for(factor = 0; factor < type.n; factor++)
    if(simple_root < offset + type.factors[factor].rank)
      break;
    else
      offset += type.factors[factor].rank;
  r = simple_root - offset;

  switch(type.factors[factor].series) {
  case 'A':
    iota_r = type.factors[factor].rank - 1 - r;
    break;

  case 'B': case 'C':
    iota_r = r;
    break;

  case 'D':
    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
    break;

  case 'E':
    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
    break;

  case 'F':
    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
    break;

  case 'G':
    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
    break;

  default:
    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
  }

  return iota_r + offset;
}

void weyl_cartan_matrix(semisimple_type_t type, int *m)
{
  int offset = 0;
  int rank = weyl_rank(type);

  int **A = (int**)malloc(rank*sizeof(int*));

  memset(m, 0, rank*rank*sizeof(int));
  for(int i = 0; i < rank; i++)
    m[i*rank+i] = 2;

  for(int k = 0; k < type.n; k++) {
    for(int i = 0; i < type.factors[k].rank; i++)  // A is the submatrix corresponding to the current simple factor
      A[i] = &m[(i+offset)*rank + offset];

    switch(type.factors[k].series) {
    case 'A':
      for(int i = 1; i < type.factors[k].rank; i++) {
	A[i][i-1] = -1;
	A[i-1][i] = -1;
      }
      break;

    case 'B':        // not sure at all about the order of B and C
      if(type.factors[k].rank >= 2) {
	A[0][1] = -2;
	A[1][0] = -1;
      }
      for(int i = 2; i < type.factors[k].rank; i++) {
	A[i][i-1] = -1;
	A[i-1][i] = -1;
      }
      break;

    case 'C':
      if(type.factors[k].rank >= 2) {
	A[0][1] = -1;
	A[1][0] = -2;
      }
      for(int i = 2; i < type.factors[k].rank; i++) {
	A[i][i-1] = -1;
	A[i-1][i] = -1;
      }
      break;

    case 'D':
      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
      break;

    case 'E':
      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
      break;

    case 'F':
      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
      break;

    case 'G':
      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
      break;

    default:
      ERROR(1, "A Weyl 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;
  }

  free(A);
}

/************ memory allocation ********************/

weylgroup_element_t *weyl_alloc(semisimple_type_t type)
{
  int rank = weyl_rank(type);
  int order = weyl_order(type);

  int *left = (int*)malloc(rank*order*sizeof(int));
  int *right = (int*)malloc(rank*order*sizeof(int));
  weylgroup_element_t *group = (weylgroup_element_t*)malloc(order*sizeof(weylgroup_element_t));

  for(int i = 0; i < order; i++) {
    group[i].left = &left[i*rank];
    group[i].right = &right[i*rank];
  }

  return group;
}

void weyl_free(weylgroup_element_t *x)
{
  free(x[0].left);
  free(x[0].right);
  free(x);
}

void weyl_generate(semisimple_type_t type, weylgroup_element_t *group)
{
  int rank, order, positive;
  queue_t queue;
  int current;
  int roots_known, elements, length_elements, nextids_count;
  int *cartan_matrix;
  int *root_vectors;
  int *vector;
  int *simple_roots;
  int *root_mapping;
  weylid_t *ids, *edges, *nextids;
  weylid_lookup_t *lookup;

  rank = weyl_rank(type);
  order = weyl_order(type);
  positive = weyl_hyperplanes(type);

  ERROR(positive > 64, "We can't handle root systems with more than 64 positive roots!\n");

  cartan_matrix =      (int*)malloc(rank*rank      *sizeof(int));
  root_vectors =       (int*)malloc(2*positive*rank*sizeof(int));
  vector =             (int*)malloc(rank           *sizeof(int));
  root_mapping =       (int*)malloc(positive*rank  *sizeof(int));
  simple_roots =       (int*)malloc(rank           *sizeof(int));
  ids =           (weylid_t*)malloc(order          *sizeof(weylid_t));
  edges =         (weylid_t*)malloc(rank*order     *sizeof(weylid_t));
  nextids =       (weylid_t*)malloc(rank*order     *sizeof(weylid_t));
  lookup = (weylid_lookup_t*)malloc(order          *sizeof(weylid_lookup_t));

  weyl_cartan_matrix(type, cartan_matrix);

  // enumerate roots
  memset(root_vectors, 0, 2*positive*rank*sizeof(int));

  // first the simple roots
  queue_init(&queue);
  for(int i = 0; i < rank; i++) {
    root_vectors[rank*i + i] = 1;
    queue_put(&queue, i);
  }

  // and then we get all others by reflecting
  roots_known = rank;
  while((current = queue_get(&queue)) != -1) {
    for(int i = 0; i < rank; i++) {
      reflect_root_vector(cartan_matrix, rank, i, &root_vectors[rank*current], vector);
      int j;
      for(j = 0; j < roots_known; j++)
	if(compare_root_vectors(rank, &root_vectors[rank*j], vector) == 0)
	  break;
      if(j == roots_known) {
	memcpy(&root_vectors[rank*roots_known], vector, rank*sizeof(int));
	queue_put(&queue, roots_known);
	roots_known++;
      }
    }
  }

  ERROR(roots_known != 2*positive, "Number of roots does not match!\n")

  // sort roots and restrict to positives
  qsort_r(root_vectors, 2*positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
  memcpy(root_vectors, &root_vectors[positive*rank], positive*rank*sizeof(int));

  for(int i = 0; i < positive; i++) {
    for(int j = 0; j < rank; j++) {
      reflect_root_vector(cartan_matrix, rank, j, &root_vectors[rank*i], vector);
      root_mapping[i*rank+j] =
	search(vector, root_vectors, positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
    }
  }

  // where in the list are the simple roots?
  for(int i = 0; i < rank; i++) {
    memset(vector, 0, rank*sizeof(int));
    vector[i] = 1;
    simple_roots[i] = search(vector, root_vectors, positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
  }

  // enumerate weyl group elements using difference sets
  nextids[0] = 0;
  nextids_count = 1;
  elements = 0;
  for(int len = 0; len <= positive; len++) {
    length_elements = 0;

    // find unique ids in edges added in the last iteration
    qsort(nextids, nextids_count, sizeof(weylid_t), compare_weylid);
    for(int i = 0; i < nextids_count; i++)
      if(i == 0 || nextids[i] != nextids[i-1])
	ids[elements + length_elements++] = nextids[i];

    // add new edges
    nextids_count = 0;
    for(int i = elements; i < elements + length_elements; i++)
      for(int j = 0; j < rank; j++) {
	edges[i*rank+j] = multiply_generator(j, ids[i], simple_roots, root_mapping, rank, positive);
	if(!(ids[i] & BIT(simple_roots[j]))) // the new element is longer then the old one
	  nextids[nextids_count++] = edges[i*rank+j];
      }

    elements += length_elements;
  }

  // translate the ids to list positions (i.e. local continuous ids)
  for(int i = 0; i < order; i++) {
    lookup[i].id = ids[i];
    lookup[i].position = i;
  }
  qsort(lookup, order, sizeof(weylid_lookup_t), compare_weylid_lookup);

  for(int i = 0; i < order; i++) {
    group[i].id = ids[i];
    for(int j = 0; j < rank; j++)
      group[i].left[j] = lookup_id(edges[i*rank+j], lookup, order);
  }

  free(cartan_matrix);
  free(root_vectors);
  free(vector);
  free(root_mapping);
  free(simple_roots);
  free(ids);
  free(edges);
  free(nextids);
  free(lookup);
}