florian
/
triangle_group_limit_set
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

limit set cairo

This commit is contained in:
Florian Stecker 2018-08-08 16:24:03 +02:00
commit a373a8ccf7
9 changed files with 1040 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
Makefile Normal file
View File

@ -0,0 +1,30 @@
HEADERS=triangle.h linalg.h queue.h initcairo.h
#SPECIAL_OPTIONS=-O0 -g -D_DEBUG
#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -Wall -Wno-unused-function -Werror=implicit-function-declaration
#SPECIAL_OPTIONS=
CAIRO_OPTIONS=$(shell pkg-config --cflags cairo)
GENERAL_OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE
OPTIONS=$(GENERAL_OPTIONS) $(CAIRO_OPTIONS) $(SPECIAL_OPTIONS)
all: limit_set
limit_set: limit_set.o linalg.o triangle.o initcairo.o
gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o -lm -lgsl -lcblas -lcairo -lX11
linalg.o: linalg.c $(HEADERS)
gcc $(OPTIONS) -c linalg.c
triangle.o: triangle.c $(HEADERS)
gcc $(OPTIONS) -c triangle.c
limit_set.o: limit_set.c $(HEADERS)
gcc $(OPTIONS) -c limit_set.c
initcairo.o: initcairo.c $(HEADERS)
gcc $(OPTIONS) -c initcairo.c
clean:
rm -f limit_set linalg.o triangle.o limit_set.o

203
initcairo.c Normal file
View File

@ -0,0 +1,203 @@
#include <sys/stat.h>
#include <sys/time.h>
#include <math.h>
#include <cairo-xlib.h>
#include "initcairo.h"
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); return 0;}
static Bool alwaysTruePredicate(Display *display, XEvent *event, XPointer arg)
{
return true;
}
int initCairo(int screen, int mask, int width, int height, const char *name, GraphicsInfo *info)
{
mask |= StructureNotifyMask | ExposureMask | KeyPressMask | ButtonPressMask | PointerMotionMask;
info->display = XOpenDisplay(NULL);
ERROR(!info->display, "Failed to open X display\n");
Visual *visual = XDefaultVisual(info->display, screen);
int depth = XDefaultDepth(info->display, screen);
Window root = XRootWindow(info->display, screen);
XSetWindowAttributes swa;
info->cmap = XCreateColormap(info->display, root, visual, AllocNone);
swa.colormap = info->cmap;
swa.background_pixmap = None;
swa.border_pixel = 0;
swa.event_mask = mask;
info->win = XCreateWindow(info->display, root, 0, 0, width, height, 0, depth, InputOutput, visual, CWBorderPixel|CWColormap|CWEventMask, &swa);
ERROR(!info->win, "Failed to create window.\n");
XStoreName(info->display, info->win, name);
XMapWindow(info->display, info->win);
info->sfc = cairo_xlib_surface_create(info->display, info->win, visual, width, height);
cairo_xlib_surface_set_size(info->sfc, width, height);
info->ctx = cairo_create(info->sfc);
cairo_matrix_init_identity(&info->matrix);
info->scalefactor = 1.0;
info->wm_delete_window = XInternAtom(info->display, "WM_DELETE_WINDOW", 0);
info->wm_protocols = XInternAtom(info->display, "WM_PROTOCOLS", 0);
XSetWMProtocols(info->display, info->win, &info->wm_delete_window, 1);
return 1;
}
void destroyCairo(GraphicsInfo *info)
{
XDestroyWindow(info->display, info->win);
XFreeColormap(info->display, info->cmap);
XCloseDisplay(info->display);
}
void startTimer(GraphicsInfo *info)
{
gettimeofday(&info->start_time, NULL);
info->frames = 0;
}
void waitUpdateTimer(GraphicsInfo *info)
{
struct timeval current_time;
double new_elapsed;
gettimeofday(&current_time, NULL);
new_elapsed = current_time.tv_sec - info->start_time.tv_sec;
new_elapsed += (current_time.tv_usec - info->start_time.tv_usec) * 1e-6;
info->frametime = new_elapsed - info->elapsed;
if(info->frametime < 0.01) { // frames < 10ms are considered too short; sleep a while and then measure again
usleep(10000 - info->frametime*1e6);
gettimeofday(&current_time, NULL);
new_elapsed = current_time.tv_sec - info->start_time.tv_sec;
new_elapsed += (current_time.tv_usec - info->start_time.tv_usec) * 1e-6;
info->frametime = new_elapsed - info->elapsed;
}
info->elapsed = new_elapsed;
info->frames++;
}
int processStandardEvent(GraphicsInfo *info, XEvent *ev, void (*draw)(void*), void *data)
{
int state;
static int last_x, last_y;
switch(ev->type) {
case Expose:
draw(data);
break;
case ConfigureNotify:
printf("ConfigureNotify Event, new dimensions: %d %d %d %d\n", ev->xconfigure.x, ev->xconfigure.y, ev->xconfigure.width, ev->xconfigure.height);
info->width = ev->xconfigure.width;
info->height = ev->xconfigure.height;
cairo_xlib_surface_set_size(info->sfc, info->width, info->height);
break;
case KeyPress:
state = ev->xkey.state & (ShiftMask | LockMask | ControlMask);
if(state == 0 && ev->xkey.keycode == 24)
return true;
break;
case ButtonPress:
if(ev->xbutton.button == 4) {
cairo_matrix_t transform;
cairo_matrix_init_identity(&transform);
cairo_matrix_translate(&transform, ev->xbutton.x, ev->xbutton.y);
cairo_matrix_scale(&transform, 5.0/4.0, 5.0/4.0);
cairo_matrix_translate(&transform, -ev->xbutton.x, -ev->xbutton.y);
cairo_matrix_multiply(&info->matrix, &info->matrix, &transform);
info->scalefactor *= 5.0/4.0;
} else if(ev->xbutton.button == 5) {
cairo_matrix_t transform;
cairo_matrix_init_identity(&transform);
cairo_matrix_translate(&transform, ev->xbutton.x, ev->xbutton.y);
cairo_matrix_scale(&transform, 4.0/5.0, 4.0/5.0);
cairo_matrix_translate(&transform, -ev->xbutton.x, -ev->xbutton.y);
cairo_matrix_multiply(&info->matrix, &info->matrix, &transform);
info->scalefactor *= 4.0/5.0;
}
last_x = ev->xbutton.x;
last_y = ev->xbutton.y;
draw(data);
break;
case MotionNotify:
if(ev->xmotion.state & Button1Mask && !(ev->xmotion.state & ShiftMask)) {
int dx = ev->xmotion.x - last_x;
int dy = ev->xmotion.y - last_y;
cairo_matrix_t transform;
cairo_matrix_init_identity(&transform);
cairo_matrix_translate(&transform, dx, dy);
cairo_matrix_multiply(&info->matrix, &info->matrix, &transform);
printf("Button1Motion Event, dx: %d, dy: %d\n", dx, dy);
draw(data);
} else if(ev->xmotion.state & Button1Mask && ev->xmotion.state & ShiftMask) {
double angle =
atan2((double)ev->xmotion.y - (double)info->height/2, (double)ev->xmotion.x - (double)info->width/2)-
atan2((double)last_y - (double)info->height/2, (double)last_x - (double)info->width/2);
cairo_matrix_t transform;
cairo_matrix_init_identity(&transform);
cairo_matrix_translate(&transform, (double)info->width/2, (double)info->height/2);
cairo_matrix_rotate(&transform, angle);
cairo_matrix_translate(&transform, -(double)info->width/2, -(double)info->height/2);
cairo_matrix_multiply(&info->matrix, &info->matrix, &transform);
printf("Button1Motion Event, angle: %f\n", angle);
draw(data);
}
last_x = ev->xmotion.x;
last_y = ev->xmotion.y;
break;
case ClientMessage:
if((Atom)ev->xclient.message_type == info->wm_protocols && (Atom)ev->xclient.data.l[0] == info->wm_delete_window) {
// printf("Window closed\n");
return true;
}
break;
}
return 0;
}
int checkEvents(GraphicsInfo *info, int (*process)(GraphicsInfo*, XEvent*, void*), void (*draw)(void*), void *data) // get any events from the queue and the server, process them if neccessary, quit if wanted
{
XEvent ev;
// while(XCheckIfEvent(info->display, &ev, alwaysTruePredicate, NULL)) { // we essentially want XCheckWindowEvent, but we want to avoid that events for other windows fill up the queue
XNextEvent(info->display, &ev);
if(ev.xany.window != info->win)
return 0;
/* if(ev.type == KeyRelease) { // deal with autorepeat
XEvent nev;
if(XCheckIfEvent(info->display, &nev, alwaysTruePredicate, NULL)) { // is there another event?
if (nev.type == KeyPress && nev.xkey.time == ev.xkey.time && nev.xkey.keycode == ev.xkey.keycode) // which is equal, but KeyPress? Then it's just auto-repeat
return 0; // so we ignore both
XPutBackEvent(info->display, &nev); // otherwise put the event back, we will consider it in the next round
}
} */
if(processStandardEvent(info, &ev, draw, data))
return 1;
if(process(info, &ev, data))
return 1; // quit event queue and application
return 0;
}

39
initcairo.h Normal file
View File

@ -0,0 +1,39 @@
#ifndef INITCAIRO_H
#define INITCAIRO_H
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <unistd.h>
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <cairo.h>
#include <cairo-xlib.h>
typedef struct {
Display *display;
Window win;
Colormap cmap;
Atom wm_protocols;
Atom wm_delete_window;
cairo_surface_t *sfc;
cairo_t *ctx;
struct timeval start_time;
unsigned long frames;
double elapsed, frametime;
unsigned int width;
unsigned int height;
cairo_matrix_t matrix;
double scalefactor;
} GraphicsInfo;
int initCairo(int screen, int mask, int width, int height, const char *name, GraphicsInfo *info);
void destroyCairo(GraphicsInfo *info);
void startTimer(GraphicsInfo *info);
void waitUpdateTimer(GraphicsInfo *info);
int checkEvents(GraphicsInfo *info, int (*process)(GraphicsInfo*, XEvent*, void*), void (*draw)(void*), void *data);
#endif

353
limit_set.c Normal file
View File

@ -0,0 +1,353 @@
#include <math.h>
#include <stdio.h>
#include <sys/time.h>
#include <X11/XKBlib.h>
#include "initcairo.h"
#include "triangle.h"
#include "linalg.h"
#define LOOP(i) for(int i = 0; i < 3; i++)
typedef struct {
double x;
double y;
} point_t;
int processEvent(GraphicsInfo *info, XEvent *ev, void *data);
void setup();
void destroy();
void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s);
void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan);
int compareAngle(const void *x, const void *y);
int computeLimitCurve();
void draw(void *data);
point_t intersect(point_t a, point_t b, point_t c, point_t d);
GraphicsInfo G;
double parameter = 2.8;
int n_group_elements;
double *limit_curve; // x, y, angle triples
int limit_curve_valid = 0;
groupelement_t *group;
workspace_t *ws;
gsl_matrix *gen[3];
gsl_matrix **matrices;
gsl_matrix *cartan, *cob, *coxeter, *coxeter_fixedpoints, *fixedpoints;
int processEvent(GraphicsInfo *info, XEvent *ev, void *data)
{
int state;
unsigned long key;
switch(ev->type) {
case KeyPress:
state = ev->xkey.state & (ShiftMask | LockMask | ControlMask);
key = XkbKeycodeToKeysym(ev->xkey.display, ev->xkey.keycode, 0, !!(state & ShiftMask));
printf("Key pressed: %ld\n", key);
switch(key) {
case '<':
parameter /= exp(0.002);
limit_curve_valid = computeLimitCurve();
break;
case '>':
parameter *= exp(0.002);
limit_curve_valid = computeLimitCurve();
break;
case ' ':
parameter = 2.890053638;
limit_curve_valid = computeLimitCurve();
break;
case XK_Return:
parameter = 2.76375163;
limit_curve_valid = computeLimitCurve();
break;
}
draw(0);
}
return 0;
}
void setup()
{
n_group_elements = 10000;
limit_curve = malloc(3*n_group_elements*sizeof(double));
group = malloc(n_group_elements*sizeof(groupelement_t));
ws = workspace_alloc(3);
LOOP(i) gen[i] = gsl_matrix_alloc(3, 3);
matrices = malloc(n_group_elements*sizeof(gsl_matrix*));
for(int i = 0; i < n_group_elements; i++)
matrices[i] = gsl_matrix_alloc(3, 3);
cartan = gsl_matrix_alloc(3, 3);
cob = gsl_matrix_alloc(3, 3);
coxeter = gsl_matrix_alloc(3, 3);
coxeter_fixedpoints = gsl_matrix_alloc(3, 3);
fixedpoints = gsl_matrix_alloc(3, 3);
}
void destroy()
{
free(limit_curve);
free(group);
workspace_free(ws);
LOOP(i) gsl_matrix_free(gen[i]);
for(int i = 0; i < n_group_elements; i++)
gsl_matrix_free(matrices[i]);
free(matrices);
gsl_matrix_free(cartan);
gsl_matrix_free(cob);
gsl_matrix_free(coxeter);
gsl_matrix_free(coxeter_fixedpoints);
gsl_matrix_free(fixedpoints);
}
void cartanMatrix(gsl_matrix *cartan, double a1, double a2, double a3, double s)
{
gsl_matrix_set(cartan, 0, 0, -2);
gsl_matrix_set(cartan, 0, 1, 2*s*cos(a3));
gsl_matrix_set(cartan, 0, 2, 2/s*cos(a2));
gsl_matrix_set(cartan, 1, 0, 2/s*cos(a3));
gsl_matrix_set(cartan, 1, 1, -2);
gsl_matrix_set(cartan, 1, 2, 2*s*cos(a1));
gsl_matrix_set(cartan, 2, 0, 2*s*cos(a2));
gsl_matrix_set(cartan, 2, 1, 2/s*cos(a1));
gsl_matrix_set(cartan, 2, 2, -2);
}
void initializeTriangleGenerators(gsl_matrix **gen, gsl_matrix *cartan)
{
LOOP(i) gsl_matrix_set_identity(gen[i]);
LOOP(i) LOOP(j) *gsl_matrix_ptr(gen[i], i, j) += gsl_matrix_get(cartan, i, j);
}
int compareAngle(const void *x, const void *y)
{
return ((double*)x)[2] > ((double*)y)[2] ? 1 : -1;
}
// intersect the lines a-b and c-d
point_t intersect(point_t a, point_t b, point_t c, point_t d)
{
point_t res;
double t = ((c.y-d.y)*(c.x-a.x) - (c.x-d.x)*(c.y-a.y)) / ((b.x-a.x)*(c.y-d.y) - (b.y-a.y)*(c.x-d.x));
res.x = (1-t)*a.x + t*b.x;
res.y = (1-t)*a.y + t*b.y;
return res;
}
void draw_segment(point_t a, point_t b)
{
cairo_move_to(G.ctx, a.x, a.y);
cairo_line_to(G.ctx, b.x, b.y);
cairo_stroke(G.ctx);
}
void draw_polygon(int sides, point_t a, point_t b, point_t c, point_t d) // only quadrilaterals so far
{
draw_segment(a, b);
draw_segment(b, c);
draw_segment(c, d);
draw_segment(d, a);
}
void draw_box(const char *word, char base)
{
// lower stack[i] reserved for linalg.c is used by multiply_*()
gsl_matrix *tmp1 = ws->stack[10];
gsl_matrix *conj = ws->stack[12];
gsl_matrix *conjinv = ws->stack[13];
gsl_matrix *coxeter[3] = { ws->stack[14], ws->stack[15], ws->stack[16] }; // fixed points of coxeter elements
gsl_matrix_set_identity(conj);
gsl_matrix_set_identity(conjinv);
for(int i = 0; i < strlen(word); i++) {
if(word[i] == ' ')
continue;
multiply_right(conj, gen[word[i]-'a'], ws);
multiply_left(gen[word[i]-'a'], conjinv, ws);
}
for(int i = 0; i < 3; i++) {
multiply_many(ws, tmp1, 5, conj, gen[(base-'A'+i)%3], gen[(base-'A'+i+1)%3], gen[(base-'A'+i+2)%3], conjinv);
eigenvectors(tmp1, coxeter[i], ws);
multiply_left(cob, coxeter[i], ws);
}
point_t p_a[3],p_m[3],p_r[3],p_i[4];
LOOP(i) p_a[i].x = gsl_matrix_get(coxeter[i], 0, 0) / gsl_matrix_get(coxeter[i], 2, 0);
LOOP(i) p_a[i].y = gsl_matrix_get(coxeter[i], 1, 0) / gsl_matrix_get(coxeter[i], 2, 0);
LOOP(i) p_m[i].x = gsl_matrix_get(coxeter[i], 0, 1) / gsl_matrix_get(coxeter[i], 2, 1);
LOOP(i) p_m[i].y = gsl_matrix_get(coxeter[i], 1, 1) / gsl_matrix_get(coxeter[i], 2, 1);
LOOP(i) p_r[i].x = gsl_matrix_get(coxeter[i], 0, 2) / gsl_matrix_get(coxeter[i], 2, 2);
LOOP(i) p_r[i].y = gsl_matrix_get(coxeter[i], 1, 2) / gsl_matrix_get(coxeter[i], 2, 2);
p_i[0] = intersect(p_a[1], p_r[1], p_m[0], p_r[0]);
p_i[1] = intersect(p_a[0], p_r[0], p_a[2], p_r[2]);
p_i[2] = intersect(p_a[1], p_m[1], p_a[2], p_r[2]);
p_i[3] = intersect(p_a[0], p_r[0], p_a[1], p_m[1]);
draw_polygon(4, p_a[1], p_i[0], p_r[0], p_i[3]);
}
int computeLimitCurve()
{
int p = 5, q = 5, r = 5;
int k = 2, l = 2, m = 2;
generate_triangle_group(group, n_group_elements, p, q, r);
// do first in the Fuchsian case to get the angles
cartanMatrix(cartan, M_PI/p, M_PI/q, M_PI/r, 1.0);
initializeTriangleGenerators(gen, cartan);
gsl_matrix_set_identity(matrices[0]);
for(int i = 1; i < n_group_elements; i++)
multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
diagonalize_symmetric_form(cartan, cob, ws);
multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]);
if(!eigenvectors(coxeter, coxeter_fixedpoints, ws))
return 0;
for(int i = 0; i < n_group_elements; i++) {
multiply_many(ws, fixedpoints, 3, cob, matrices[i], coxeter_fixedpoints);
limit_curve[3*i+2] = atan2(
gsl_matrix_get(fixedpoints, 0, 0)/gsl_matrix_get(fixedpoints, 2, 0),
gsl_matrix_get(fixedpoints, 1, 0)/gsl_matrix_get(fixedpoints, 2, 0));
}
// now to it again to calculate x and y coordinates
cartanMatrix(cartan, M_PI*k/p, M_PI*l/q, M_PI*m/r, parameter);
initializeTriangleGenerators(gen, cartan);
gsl_matrix_set_identity(matrices[0]);
for(int i = 1; i < n_group_elements; i++)
multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
gsl_matrix_memcpy(cob, cartan); // is this a good choice of basis
// gsl_matrix_set_identity(cob);
multiply_many(ws, coxeter, 3, gen[0], gen[1], gen[2]);
if(!eigenvectors(coxeter, coxeter_fixedpoints, ws))
return 0;
for(int i = 0; i < n_group_elements; i++) {
multiply_many(ws, fixedpoints, 3, cob, matrices[i], coxeter_fixedpoints);
limit_curve[3*i] = gsl_matrix_get(fixedpoints, 0, 0)/gsl_matrix_get(fixedpoints, 2, 0);
limit_curve[3*i+1] = gsl_matrix_get(fixedpoints, 1, 0)/gsl_matrix_get(fixedpoints, 2, 0);
}
qsort(limit_curve, n_group_elements, 3*sizeof(double), compareAngle);
return 1;
}
void draw(void *data)
{
struct timeval current_time;
gettimeofday(&current_time, 0);
double start_time = current_time.tv_sec + current_time.tv_usec*1e-6;
cairo_push_group(G.ctx);
cairo_set_source_rgb(G.ctx, 1, 1, 1);
cairo_paint(G.ctx);
cairo_set_matrix(G.ctx, &G.matrix);
if(limit_curve_valid) {
int last_inside = 0;
for(int i = 0; i < n_group_elements; i++) {
double x = limit_curve[3*i];
double y = limit_curve[3*i+1];
cairo_user_to_device(G.ctx, &x, &y);
if(-x < G.width && x < 3*G.width && -y < G.height && y < 3*G.height) {
if(!last_inside) {
cairo_move_to(G.ctx, limit_curve[3*i], limit_curve[3*i+1]);
} else {
cairo_line_to(G.ctx, limit_curve[3*i], limit_curve[3*i+1]);
}
last_inside = 1;
} else {
last_inside = 0;
}
}
cairo_set_line_width(G.ctx, 1.0/G.scalefactor);
cairo_set_line_join(G.ctx, CAIRO_LINE_JOIN_BEVEL);
cairo_set_source_rgb(G.ctx, 0, 0, 0);
cairo_stroke(G.ctx);
cairo_set_source_rgb(G.ctx, 1, 0, 0);
// draw_box("c", 'C');
draw_box("", 'B');
draw_box("a", 'A');
// draw_box("", 'C');
cairo_set_source_rgb(G.ctx, 0, 0, 1);
draw_box("ca", 'A');
draw_box("cac", 'C');
draw_box("caca", 'A');
draw_box("cacac", 'C');
draw_box("acac", 'A');
draw_box("aca", 'C');
cairo_set_source_rgb(G.ctx, 0, 1, 0);
draw_box("ac", 'C');
draw_box("aca", 'A');
draw_box("acac", 'C');
draw_box("acaca", 'A');
draw_box("caca", 'C');
draw_box("cac", 'A');
}
cairo_identity_matrix(G.ctx);
cairo_move_to(G.ctx, 15, 30);
cairo_set_source_rgb(G.ctx, 0, 0, 0);
cairo_set_font_size(G.ctx, 16);
char buf[100];
sprintf(buf, "t = %.8f", parameter);
cairo_show_text(G.ctx, buf);
cairo_pop_group_to_source(G.ctx);
cairo_paint(G.ctx);
cairo_surface_flush(G.sfc);
gettimeofday(&current_time, 0);
double end_time = current_time.tv_sec + current_time.tv_usec*1e-6;
printf("draw() finished in %g seconds\n", end_time - start_time);
}
int main()
{
setup();
limit_curve_valid = computeLimitCurve();
if(!initCairo(0, KeyPressMask, 200, 200, "Triangle group", &G))
return 1;
cairo_matrix_init_identity(&G.matrix);
G.matrix.xx = G.matrix.yy = G.scalefactor = 1100.0;
G.matrix.x0 = 1150.0;
G.matrix.y0 = 900.0;
startTimer(&G);
while(!checkEvents(&G, processEvent, draw, 0)) {
waitUpdateTimer(&G);
}
destroyCairo(&G);
destroy();
return 0;
}

228
linalg.c Normal file
View File

@ -0,0 +1,228 @@
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <memory.h>
#include "linalg.h"
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
/*********************************************** temporary storage ********************************************************/
workspace_t *workspace_alloc(int n)
{
workspace_t *result = (workspace_t*)malloc(sizeof(workspace_t));
result->n = n;
result->work_nonsymmv = gsl_eigen_nonsymmv_alloc(n);
result->work_symmv = gsl_eigen_symmv_alloc(n);
result->work_sv = gsl_vector_alloc(n);
result->eval_complex = gsl_vector_complex_alloc(n);
result->evec_complex = gsl_matrix_complex_alloc(n, n);
result->eval_real = gsl_vector_alloc(n);
result->evec_real = gsl_matrix_alloc(n, n);
result->tmp = gsl_matrix_alloc(n, n);
result->permutation = gsl_permutation_alloc(n);
for(int i = 0; i < 20; i++)
result->stack[i] = gsl_matrix_alloc(n, n);
return result;
}
void workspace_free(workspace_t *workspace)
{
gsl_eigen_nonsymmv_free(workspace->work_nonsymmv);
gsl_eigen_symmv_free(workspace->work_symmv);
gsl_vector_free(workspace->work_sv);
gsl_vector_complex_free(workspace->eval_complex);
gsl_matrix_complex_free(workspace->evec_complex);
gsl_vector_free(workspace->eval_real);
gsl_matrix_free(workspace->evec_real);
gsl_matrix_free(workspace->tmp);
gsl_permutation_free(workspace->permutation);
for(int i = 0; i < 20; i++)
gsl_matrix_free(workspace->stack[i]);
}
/************************************************** basic operations ********************************************************/
void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws)
{
int s;
gsl_matrix_memcpy(ws->tmp, in);
gsl_linalg_LU_decomp(ws->tmp, ws->permutation, &s);
gsl_linalg_LU_invert(ws->tmp, ws->permutation, out);
}
void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws)
{
invert(conjugator, out, ws); // use out to temporarily store inverse conjugator
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, in, out, 0.0, ws->tmp); // in * conjugator^{-1}
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, conjugator, ws->tmp, 0.0, out);
}
void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out)
{
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, out);
}
void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws)
{
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, ws->stack[0]);
gsl_matrix_memcpy(a, ws->stack[0]);
}
void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws)
{
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, ws->stack[0]);
gsl_matrix_memcpy(b, ws->stack[0]);
}
void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...)
{
va_list args;
va_start(args, n);
gsl_matrix_set_identity(out);
for(int i = 0; i < n; i++) {
gsl_matrix *cur = va_arg(args, gsl_matrix *);
multiply_right(out, cur, ws);
}
va_end(args);
}
void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
{
gsl_matrix_memcpy(ws->tmp, g);
gsl_linalg_SV_decomp(ws->tmp, ws->evec_real, ws->eval_real, ws->work_sv);
for(int i = 0; i < ws->n - 1; i++)
mu[i] = log(gsl_vector_get(ws->eval_real, i) / gsl_vector_get(ws->eval_real, i+1));
}
void initialize(gsl_matrix *g, double *data, int x, int y)
{
gsl_matrix_view view = gsl_matrix_view_array(data, x, y);
gsl_matrix_memcpy(g, &view.matrix);
}
void rotation_matrix(gsl_matrix *g, double *vector)
{
double normalized[3];
double norm = sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]);
for(int i = 0; i < 3; i++)
normalized[i] = vector[i] / norm;
gsl_matrix_set_identity(g);
gsl_matrix_set(g, 0, 0, cos(norm));
gsl_matrix_set(g, 0, 1, -sin(norm) * normalized[2]);
gsl_matrix_set(g, 0, 2, +sin(norm) * normalized[1]);
gsl_matrix_set(g, 1, 0, +sin(norm) * normalized[2]);
gsl_matrix_set(g, 1, 1, cos(norm));
gsl_matrix_set(g, 1, 2, -sin(norm) * normalized[0]);
gsl_matrix_set(g, 2, 0, -sin(norm) * normalized[1]);
gsl_matrix_set(g, 2, 1, +sin(norm) * normalized[0]);
gsl_matrix_set(g, 2, 2, cos(norm));
for(int i = 0; i < 3; i++)
for(int j = 0; j < 3; j++)
g->data[i * g->tda + j] += (1 - cos(norm)) * normalized[i] * normalized[j];
}
double trace(gsl_matrix *g)
{
return gsl_matrix_get(g, 0, 0) + gsl_matrix_get(g, 1, 1) + gsl_matrix_get(g, 2, 2);
}
double determinant(gsl_matrix *g, workspace_t *ws)
{
int s;
gsl_matrix_memcpy(ws->tmp, g);
gsl_linalg_LU_decomp(ws->tmp, ws->permutation, &s);
return gsl_linalg_LU_det(ws->tmp, s);
}
int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
{
gsl_eigen_nonsymmv_params(1, ws->work_nonsymmv);
gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
int real = 1;
for(int i = 0; i < 4; i++)
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
real = 0;
if(!real)
return 0; // non-real eigenvalues!
for(int i = 0; i < ws->n - 1; i++) {
mu[i] = log(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i+1)));
}
return 1;
}
int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
{
gsl_matrix_memcpy(ws->stack[0], g);
gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
int r = gsl_eigen_nonsymmv(ws->stack[0], 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 = 1;
for(int i = 0; i < ws->n; i++)
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
real = 0;
if(!real)
return 0; // non-real eigenvalues!
for(int i = 0; i < ws->n; i++)
for(int j = 0; j < ws->n; j++)
gsl_matrix_set(evec_real, i, j, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
return 1;
}
void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws)
{
gsl_matrix_memcpy(ws->stack[0], g);
int r = gsl_eigen_symmv (ws->stack[0], eval, evec, ws->work_symmv);
ERROR(r, "gsl_eigen_symmv failed!\n");
gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
}
// returns number of positive directions and matrix transforming TO diagonal basis
int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
{
gsl_matrix_memcpy(ws->stack[0], A);
int r = gsl_eigen_symmv (ws->stack[0], ws->eval_real, cob, ws->work_symmv);
ERROR(r, "gsl_eigen_symmv failed!\n");
gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC);
gsl_matrix_transpose(cob);
int positive = 0;
for(int i = 0; i < ws->n; i++) {
if(gsl_vector_get(ws->eval_real, i) > 0)
positive++;
for(int j = 0; j < ws->n; j++)
*gsl_matrix_ptr(cob, i, j) *= sqrt(fabs(gsl_vector_get(ws->eval_real, i)));
}
return positive;
// printf("Eigenvalues: %.10f, %.10f, %.10f\n", gsl_vector_get(ws->eval_real, 0), gsl_vector_get(ws->eval_real, 1), gsl_vector_get(ws->eval_real, 2));
// return 0;
}

48
linalg.h Normal file
View File

@ -0,0 +1,48 @@
#ifndef LINALG_H
#define LINALG_H
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <memory.h>
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
typedef struct _workspace {
int n;
gsl_eigen_nonsymmv_workspace *work_nonsymmv;
gsl_eigen_symmv_workspace *work_symmv;
gsl_vector *work_sv;
gsl_vector_complex *eval_complex;
gsl_matrix_complex *evec_complex;
gsl_vector *eval_real;
gsl_matrix *evec_real;
gsl_matrix *tmp;
gsl_permutation *permutation;
gsl_matrix *stack[20];
} workspace_t;
workspace_t *workspace_alloc(int n);
void workspace_free(workspace_t *workspace);
void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws);
void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws);
void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out);
void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws);
void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws);
void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...);
void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
void initialize(gsl_matrix *g, double *data, int x, int y);
void rotation_matrix(gsl_matrix *g, double *vector);
int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
double trace(gsl_matrix *g);
double determinant(gsl_matrix *g, workspace_t *ws);
int eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws);
int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws);
#endif

46
queue.h Normal file
View File

@ -0,0 +1,46 @@
#ifndef QUEUE_H
#define QUEUE_H
#include <stdio.h>
#include <stdlib.h>
#define QUEUE_SIZE 100000
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
#ifdef _DEBUG
#define LOG(msg, ...) fprintf(stderr, msg, ##__VA_ARGS__)
#else
#define LOG(msg, ...)
#endif
typedef struct {
unsigned int start;
unsigned int end;
int data[QUEUE_SIZE];
} queue_t;
static void queue_init(queue_t *q)
{
q->start = q->end = 0;
}
static void queue_put(queue_t *q, int x)
{
q->data[q->end++] = x;
q->end %= QUEUE_SIZE;
ERROR(q->start == q->end, "The queue is full! Increase QUEUE_SIZE\n");
}
static int queue_get(queue_t *q)
{
if(q->start == q->end)
return -1;
int result = q->data[q->start++];
q->start %= QUEUE_SIZE;
return result;
}
#endif

72
triangle.c Normal file
View File

@ -0,0 +1,72 @@
#include <stdio.h>
#include <string.h>
#include <memory.h>
#include "triangle.h"
#include "queue.h"
static groupelement_t *find_adj(groupelement_t *cur, int gen, int other, int order)
{
groupelement_t *next = cur->adj[other];
for(int i = 0; i < order - 1; i++) {
if(!next)
break;
next = next->adj[gen];
if(!next)
break;
next = next->adj[other];
}
return next;
}
int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3)
{
queue_t q;
int id, n;
groupelement_t *cur, *inv;
int k[3] = {k1, k2, k3};
memset(group, 0, size*sizeof(groupelement_t));
group[0].letter = -1;
n = 1;
queue_init(&q);
queue_put(&q, 0);
while((id = queue_get(&q)) != -1) {
cur = &group[id];
for(int gen = 0; gen < 3; gen++) {
if(!cur->adj[gen])
cur->adj[gen] = find_adj(cur, gen, (gen+2)%3, k[(gen+1)%3]);
if(!cur->adj[gen])
cur->adj[gen] = find_adj(cur, gen, (gen+1)%3, k[(gen+2)%3]);
if(!cur->adj[gen] && n < size) {
cur->adj[gen] = &group[n];
cur->adj[gen]->id = n;
cur->adj[gen]->parent = cur;
cur->adj[gen]->letter = gen;
cur->adj[gen]->length = cur->length + 1;
queue_put(&q, n);
n++;
}
if(cur->adj[gen])
cur->adj[gen]->adj[gen] = cur;
}
}
for(int i = 0; i < size; i++) {
cur = &group[i];
inv = &group[0];
while(cur->parent && inv) {
inv = inv->adj[cur->letter];
cur = cur->parent;
}
group[i].inverse = inv;
}
return n;
}

21
triangle.h Normal file
View File

@ -0,0 +1,21 @@
#ifndef TRIANGLE_H
#define TRIANGLE_H
#include <stdio.h>
#include <string.h>
#include <memory.h>
#include "queue.h"
typedef struct _groupelement {
int id;
int length;
struct _groupelement *adj[3];
struct _groupelement *parent;
struct _groupelement *inverse;
int letter;
} groupelement_t;
int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3);
#endif