From a373a8ccf7c047c62c4cf0beed1f476ea735e74c Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 8 Aug 2018 16:24:03 +0200 Subject: [PATCH] limit set cairo --- Makefile | 30 +++++ initcairo.c | 203 ++++++++++++++++++++++++++++++ initcairo.h | 39 ++++++ limit_set.c | 353 ++++++++++++++++++++++++++++++++++++++++++++++++++++ linalg.c | 228 +++++++++++++++++++++++++++++++++ linalg.h | 48 +++++++ queue.h | 46 +++++++ triangle.c | 72 +++++++++++ triangle.h | 21 ++++ 9 files changed, 1040 insertions(+) create mode 100644 Makefile create mode 100644 initcairo.c create mode 100644 initcairo.h create mode 100644 limit_set.c create mode 100644 linalg.c create mode 100644 linalg.h create mode 100644 queue.h create mode 100644 triangle.c create mode 100644 triangle.h diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..033c2c7 --- /dev/null +++ b/Makefile @@ -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 diff --git a/initcairo.c b/initcairo.c new file mode 100644 index 0000000..296efa1 --- /dev/null +++ b/initcairo.c @@ -0,0 +1,203 @@ +#include +#include +#include +#include + +#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(¤t_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(¤t_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; +} diff --git a/initcairo.h b/initcairo.h new file mode 100644 index 0000000..409a7c0 --- /dev/null +++ b/initcairo.h @@ -0,0 +1,39 @@ +#ifndef INITCAIRO_H +#define INITCAIRO_H + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +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 diff --git a/limit_set.c b/limit_set.c new file mode 100644 index 0000000..e527b6b --- /dev/null +++ b/limit_set.c @@ -0,0 +1,353 @@ +#include +#include +#include +#include + +#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(¤t_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(¤t_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; +} diff --git a/linalg.c b/linalg.c new file mode 100644 index 0000000..472e8ff --- /dev/null +++ b/linalg.c @@ -0,0 +1,228 @@ +#include +#include +#include +#include +#include +#include +#include + +#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; +} diff --git a/linalg.h b/linalg.h new file mode 100644 index 0000000..880d160 --- /dev/null +++ b/linalg.h @@ -0,0 +1,48 @@ +#ifndef LINALG_H +#define LINALG_H + +#include +#include +#include +#include +#include +#include +#include +#include + +#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 diff --git a/queue.h b/queue.h new file mode 100644 index 0000000..88ba8d0 --- /dev/null +++ b/queue.h @@ -0,0 +1,46 @@ +#ifndef QUEUE_H +#define QUEUE_H + +#include +#include + +#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 diff --git a/triangle.c b/triangle.c new file mode 100644 index 0000000..770c7fe --- /dev/null +++ b/triangle.c @@ -0,0 +1,72 @@ +#include +#include +#include + +#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; +} diff --git a/triangle.h b/triangle.h new file mode 100644 index 0000000..076f6d4 --- /dev/null +++ b/triangle.h @@ -0,0 +1,21 @@ +#ifndef TRIANGLE_H +#define TRIANGLE_H + +#include +#include +#include + +#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