triangle_reflection_complex/singular_values.c

#include "triangle.h"
#include "linalg.h"

#define MAX_ELEMENTS 14000
//#define DRAW_PICTURE 1

void initialize_triangle_generators(workspace_t *ws, gsl_matrix **gen, double a1, double a2, double a3, double s, double t)
{
	gsl_matrix **tmp;
	tmp = getTempMatrices(ws, 6);

	double rho1sqrt = sqrt(s*(s+2*cos(2*M_PI*a1))+1);
	double rho2sqrt = sqrt(s*(s+2*cos(2*M_PI*a2))+1);
	double rho3sqrt = sqrt(s*(s+2*cos(2*M_PI*a3))+1);

	for(int i = 0; i < 6; i++)
		gsl_matrix_set_zero(tmp[i]);

	gsl_matrix_set(tmp[0], 0, 0, 1);
	gsl_matrix_set(tmp[0], 1, 0, -rho3sqrt/t);
	gsl_matrix_set(tmp[0], 2, 0, -rho2sqrt*t);
	gsl_matrix_set(tmp[0], 1, 1, -1);
	gsl_matrix_set(tmp[0], 2, 2, -1);

	gsl_matrix_set(tmp[1], 0, 0, -1);
	gsl_matrix_set(tmp[1], 0, 1, -rho3sqrt*t);
	gsl_matrix_set(tmp[1], 1, 1, 1);
	gsl_matrix_set(tmp[1], 2, 1, -rho1sqrt/t);
	gsl_matrix_set(tmp[1], 2, 2, -1);

	gsl_matrix_set(tmp[2], 0, 0, -1);
	gsl_matrix_set(tmp[2], 1, 1, -1);
	gsl_matrix_set(tmp[2], 0, 2, -rho2sqrt/t);
	gsl_matrix_set(tmp[2], 1, 2, -rho1sqrt*t);
	gsl_matrix_set(tmp[2], 2, 2, 1);

	gsl_matrix_set(tmp[3], 0, 0, 1);
	gsl_matrix_set(tmp[3], 1, 1, s);
	gsl_matrix_set(tmp[3], 2, 2, 1/s);

	gsl_matrix_set(tmp[4], 0, 0, 1/s);
	gsl_matrix_set(tmp[4], 1, 1, 1);
	gsl_matrix_set(tmp[4], 2, 2, s);

	gsl_matrix_set(tmp[5], 0, 0, s);
	gsl_matrix_set(tmp[5], 1, 1, 1/s);
	gsl_matrix_set(tmp[5], 2, 2, 1);

	multiply_many(ws, gen[0], 3, tmp[2], tmp[3], tmp[1]);
	multiply_many(ws, gen[1], 3, tmp[0], tmp[4], tmp[2]);
	multiply_many(ws, gen[2], 3, tmp[1], tmp[5], tmp[0]);

	invert(gen[0], gen[3], ws);
	invert(gen[1], gen[4], ws);
	invert(gen[2], gen[5], ws);

	releaseTempMatrices(ws, 6);
}

char *print_word(groupelement_t *g, char *str)
{
  int i = g->length - 1;

  str[g->length] = 0;
  while(g->parent) {
    str[i--] = 'a' + g->letter;
    g = g->parent;
  }

  return str;
}

int main(int argc, char *argv[])
{
  groupelement_t *group;
  workspace_t *ws;
  gsl_matrix **matrices;
  gsl_matrix *gen[6];
  gsl_matrix *tmp, *tmp2, *coxeter;
  double mu[6], tr, trinv;
  char buf[100];

  /*
  if(argc < 2) {
    fprintf(stderr, "Too few arguments!\n");
    return 1;
  }
  */

  // allocate stuff
  group = malloc(MAX_ELEMENTS*sizeof(groupelement_t));
  ws = workspace_alloc(3);
  matrices = malloc(MAX_ELEMENTS*sizeof(gsl_matrix*));
  for(int i = 0; i < MAX_ELEMENTS; i++)
    matrices[i] = gsl_matrix_alloc(3, 3);
  for(int i = 0; i < 6; i++)
    gen[i] = gsl_matrix_alloc(3, 3);
  tmp = gsl_matrix_alloc(3, 3);
  tmp2 = gsl_matrix_alloc(3, 3);
  coxeter = gsl_matrix_alloc(3, 3);

#ifdef DRAW_PICTURE
  int width = 800;
  int height = 800;
  printf("P3\n%d %d\n255\n", width, height);

  for(int y = 0; y < height; y++) {
	  for(int x = 0; x < width; x++) {
		  double s = exp(((double)x/width-0.5)*4);
		  double t = exp(((double)y/height-0.5)*4);
#else
		  double s = atof(argv[1]);
		  double t = atof(argv[2]);
#endif

		  // the real action
		  generate_triangle_group(group, MAX_ELEMENTS, 5, 6, 7);
		  initialize_triangle_generators(ws, gen, 1.0/5.0, 1.0/6.0, 1.0/7.0, s, t);

		  gsl_matrix_set_identity(matrices[0]);
		  for(int i = 1; i < MAX_ELEMENTS; i++) {
			  if(group[i].length % 2 != 0)
				  continue;

			  int parent = group[i].parent->id;
			  int grandparent = group[i].parent->parent->id;
			  int letter;

			  if(group[parent].letter == 0 && group[i].letter == 1)
				  letter = 3; // p = ab
			  else if(group[parent].letter == 1 && group[i].letter == 2)
				  letter = 4; // q = bc
			  else if(group[parent].letter == 2 && group[i].letter == 0)
				  letter = 5; // r = ca
			  else if(group[parent].letter == 1 && group[i].letter == 0)
				  letter = 0; // p^{-1} = ba
			  else if(group[parent].letter == 2 && group[i].letter == 1)
				  letter = 1; // q^{-1} = cb
			  else if(group[parent].letter == 0 && group[i].letter == 2)
				  letter = 2; // r^{-1} = ac

			  multiply(matrices[grandparent], gen[letter], matrices[i]);
		  }

		  double excentricity;
		  double max_excentricity = 0;
		  int max_excentricity_index = 0;

		  for(int i = 0; i < MAX_ELEMENTS; i++) {
			  if(group[i].length % 2 != 0)
				  continue;

			  // jordan projection
			  gsl_matrix_memcpy(tmp, matrices[i]);
			  gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
			  int r = gsl_eigen_nonsymmv(tmp, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
			  ERROR(r, "gsl_eigen_nonsymmv failed!\n");
			  gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);

			  int real_evs = 0;

			  for(int j = 0; j < 3; j++)
				  if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, j)), 0) == 0)
					  real_evs++;

			  if(real_evs != 3)
				  continue;

			  mu[0] = fabs(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 0)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 1)));
			  mu[1] = fabs(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 1)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 2)));

			  excentricity = log(mu[1])/log(mu[0]);
			  if(excentricity > max_excentricity + 1e-3) {
				  max_excentricity = excentricity;
				  max_excentricity_index = i;
			  }

#ifndef DRAW_PICTURE
			  printf("%d %f %f %f %f 0x0000ff %d %s\n", group[i].length, mu[0], mu[1], tr, trinv, i, print_word(&group[i], buf));
#endif
		  }

#ifdef DRAW_PICTURE
		  int color = (int)((max_excentricity-1)/(max_excentricity+4)*255);
		  if(color < 0)
			  color = 0;
		  if(color > 255)
			  color = 255;
		  printf("%d %d %d\n", color, color, color);
	  }
	  fprintf(stderr, "y = %d\n", y);
  }
#else
  fprintf(stderr, "max = %f %s\n", max_excentricity, print_word(&group[max_excentricity_index], buf));
#endif

  // free stuff
  for(int i = 0; i < MAX_ELEMENTS; i++)
    gsl_matrix_free(matrices[i]);
  for(int i = 0; i < 6; i++)
    gsl_matrix_free(gen[i]);
  gsl_matrix_free(tmp);
  gsl_matrix_free(tmp2);
  gsl_matrix_free(coxeter);
  free(matrices);
  free(group);
  workspace_free(ws);

  return 0;
}