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

Compare commits

base: florian:429f0890d6b02a7aed8d42c763048e5de47441f0
Branches Tags
florian:master
..
compare: florian:3309c37955053930947b0d087e80d3552c146096
Branches Tags
florian:master

No commits in common. "429f0890d6b02a7aed8d42c763048e5de47441f0" and "3309c37955053930947b0d087e80d3552c146096" have entirely different histories.
429f0890d6 ... 3309c37955

20 changed files with 65 additions and 933 deletions
Show Stats Expand all files Collapse all files

20
.gitignore vendored
View File

@ -1,7 +1,17 @@
*.o
.\#*
\#*\#
gmon.out
core
triangle_group/singular_values
.#*
singular_values
singular_values_mpi
singular_values_barbot
output/
complex_anosov
special_element
max_slope_picture/generate
convert
billiard_words
*.pnm
*.png
*.hi
gmon.out
restart
core

54
Makefile
View File

@ -1,27 +1,50 @@
HEADERS=mat.h coxeter.h enumerate.h generators.h qext.h
HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h parallel.h qext.h
#SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT
SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT
#SPECIAL_OPTIONS=-O3 -pg -g -funroll-loops -fno-inline
SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT
#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT
#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL
#SPECIAL_OPTIONS=
OPTIONS=-I../mps/include -L../mps/lib -pthread -m64 -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
CC=gcc
all: complex_anosov
all: singular_values special_element convert billiard_words
complex_anosov: complex_anosov.o mat.o coxeter.o enumerate.o generators.o qext.o
$(CC) $(OPTIONS) -o complex_anosov -lm complex_anosov.o mat.o coxeter.o enumerate.o generators.o qext.o -lgmp -lmps
convert: convert.hs
ghc --make -dynamic convert.hs
complex_anosov.o: complex_anosov.c $(HEADERS)
gcc $(OPTIONS) -c complex_anosov.c
billiard_words: billiard_words.hs
ghc --make -dynamic billiard_words.hs
enumerate.o: enumerate.c $(HEADERS)
gcc $(OPTIONS) -c enumerate.c
singular_values: singular_values.o coxeter.o mat.o enumerate_triangle_group.o parallel.o
mpicc $(OPTIONS) -o singular_values coxeter.o singular_values.o mat.o enumerate_triangle_group.o parallel.o -lm -lgmp -lmps
generators.o: generators.c $(HEADERS)
gcc $(OPTIONS) -c generators.c
singular_values_barbot: singular_values_barbot.o coxeter.o mat.o enumerate_triangle_group.o parallel.o qext.o
mpicc $(OPTIONS) -o singular_values_barbot coxeter.o singular_values_barbot.o mat.o enumerate_triangle_group.o parallel.o qext.o -lm -lgmp -lmps
#singular_values_mpi: singular_values_mpi.o coxeter.o mat.o
# mpicc $(OPTIONS) -o singular_values_mpi coxeter.o singular_values_mpi.o mat.o -lm -lgmp -lmps
special_element: special_element.o coxeter.o linalg.o mat.o enumerate_triangle_group.o
gcc $(OPTIONS) -o special_element coxeter.o linalg.o special_element.o mat.o enumerate_triangle_group.o -lm -lgmp -lmps -lgsl -lcblas
singular_values.o: singular_values.c $(HEADERS)
gcc $(OPTIONS) -c singular_values.c
singular_values_barbot.o: singular_values_barbot.c $(HEADERS)
gcc $(OPTIONS) -c singular_values_barbot.c
#singular_values_mpi.o: singular_values_mpi.c $(HEADERS)
# mpicc $(OPTIONS) -c singular_values_mpi.c
special_element.o: special_element.c $(HEADERS)
gcc $(OPTIONS) -c special_element.c
enumerate_triangle_group.o: enumerate_triangle_group.c $(HEADERS)
gcc $(OPTIONS) -c enumerate_triangle_group.c
linalg.o: linalg.c $(HEADERS)
gcc $(OPTIONS) -c linalg.c
coxeter.o: coxeter.c $(HEADERS)
gcc $(OPTIONS) -c coxeter.c
@ -29,8 +52,11 @@ coxeter.o: coxeter.c $(HEADERS)
mat.o: mat.c $(HEADERS)
gcc $(OPTIONS) -c mat.c
parallel.o: parallel.c $(HEADERS)
gcc $(OPTIONS) -c parallel.c
qext.o: qext.c $(HEADERS)
gcc $(OPTIONS) -c qext.c
clean:
rm -f complex_anosov complex_anosov.o mat.o coxeter.o enumerate.o generators.o qext.o
rm -f singular_values singular_values_barbot special_element singular_values_mpi coxeter.o linalg.o singular_values.o singular_values_barbot.o singular_values_mpi.o mat.o special_element.o convert.hi convert.o convert billiard_words.hi billiard_words.o billiard_words enumerate_triangle_group.o parallel.o qext.o

304
complex_anosov.c
View File

@ -1,304 +0,0 @@
#include "coxeter.h"
#include "enumerate.h"
#include "generators.h"
#include "mat.h"
#include "qext.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#define LOOP(i,n) for(int i = 0; i < (n); i++)
#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
//#define INFO(msg, ...) fprintf(stderr, "[%10.3f] " msg, runtime(), ##__VA_ARGS__)
#define INFO(msg, ...)
/*
Number of elements
up to length 0: 1
up to length 1: 4
up to length 2: 10
up to length 3: 22
up to length 4: 46
up to length 5: 91
up to length 6: 175
up to length 7: 334
up to length 8: 634
up to length 9: 1198
up to length 10: 2260
up to length 11: 4261
up to length 12: 8029
up to length 13: 15124
up to length 14: 28486
up to length 15: 53650
up to length 16: 101038
up to length 17: 190279
up to length 18: 358339
up to length 19: 674830
up to length 20: 1270846
up to length 21: 2393266
up to length 22: 4507012
up to length 23: 8487625
*/
static double gaussian_sqrt5_real(NUMBER x)
{
double result = 0.0;
mpq_t tmp;
mpq_init(tmp);
// a_0 + sqrt(5)a_1 + 4a_2 + 2sqrt(5)a_3
mpq_set_si(tmp, 4, 1);
mpq_mul(tmp, tmp, x->a[2]);
mpq_add(tmp, tmp, x->a[0]);
result = mpq_get_d(tmp);
mpq_set_si(tmp, 2, 1);
mpq_mul(tmp, tmp, x->a[3]);
mpq_add(tmp, tmp, x->a[1]);
result += mpq_get_d(tmp)*sqrt(5);
mpq_clear(tmp);
return result;
}
static double gaussian_sqrt5_imag(NUMBER x)
{
double result = 0.0;
mpq_t tmp;
mpq_init(tmp);
// a_1 + 2sqrt(5)a_2 + 14a_3
mpq_set_si(tmp, 14, 1);
mpq_mul(tmp, tmp, x->a[3]);
mpq_add(tmp, tmp, x->a[1]);
result = mpq_get_d(tmp);
mpq_set_si(tmp, 2, 1);
mpq_mul(tmp, tmp, x->a[2]);
result += mpq_get_d(tmp)*sqrt(5);
mpq_clear(tmp);
return result;
}
static int read_idlist(const char *filename, int *list)
{
FILE *f = fopen(filename, "r");
if(f == NULL) {
fprintf(stderr, "Could not open %s\n", filename);
exit(1);
}
char *line = NULL;
size_t len = 0;
ssize_t read;
int n = 0;
while ((read = getline(&line, &len, f)) != -1)
list[n++] = atoi(line);
if (line)
free(line);
fclose(f);
return n;
}
static int compare_int(const void *a, const void *b)
{
const int *aa = a;
const int *bb = b;
return (*aa > *bb) - (*aa < *bb);
}
static double runtime()
{
static struct timespec starttime;
static int started = 0;
if(!started) {
clock_gettime(CLOCK_MONOTONIC, &starttime);
started = 1;
}
struct timespec curtime;
double diff;
clock_gettime(CLOCK_MONOTONIC, &curtime);
return (curtime.tv_sec - starttime.tv_sec) + (curtime.tv_nsec - starttime.tv_nsec) / 1e9;
}
enum mode {
MODE_HELP,
MODE_EIGENVALUES,
MODE_SUMMARY,
MODE_TRACE_IDS
};
int main(int argc, char *argv[])
{
char buf[100];
int mode;
runtime(); // start timer
// parse arguments
if(argc < 2 || strcmp(argv[1], "help") == 0)
mode = MODE_HELP;
else if(strcmp(argv[1], "evs") == 0)
mode = MODE_EIGENVALUES;
else if(strcmp(argv[1], "summary") == 0)
mode = MODE_SUMMARY;
else if(strcmp(argv[1], "trace_ids") == 0)
mode = MODE_TRACE_IDS;
else
mode = MODE_HELP;
if(mode == MODE_HELP) {
fprintf(stderr, "Usage: %s <help|evs|summary|trace_ids> [arguments]\n", argv[0]);
fprintf(stderr, "%s help display this page\n", argv[0]);
fprintf(stderr, "%s evs <n> <q1> <q2> <q3> <treal> <timag> enumerate group and output unique (log) eigenvalue triples\n", argv[0]);
fprintf(stderr, "%s summary <n> <q1> <q2> <q3> <treal> <timag> only output max slope etc.\n", argv[0]);
fprintf(stderr, "%s trace_ids <n> <q1> <q2> <q3> <treal> <timag> list of ids of unique traces\n", argv[0]);
return 0;
}
if(argc < 8) {
fprintf(stderr, "Not enough arguments!\n");
return 0;
}
int n = atoi(argv[2]);
int nlist, nuniq;
int q[3];
mpq_t treal, timag;
q[0] = atoi(argv[3]);
q[1] = atoi(argv[4]);
q[2] = atoi(argv[5]);
mpq_inits(treal, timag, NULL);
mpq_set_str(treal, argv[6], 10);
mpq_set_str(timag, argv[7], 10);
mpq_canonicalize(treal);
mpq_canonicalize(timag);
// enumerate group
INFO("generate group\n");
group_t *group = coxeter_init_triangle(5, 5, 5, n);
// read list of elements we need to compute, or just fill it with all elements
INFO("prepare list\n");
int *idlist = malloc(n*sizeof(int));
char *id_file = getenv("IDLIST");
if(id_file != NULL) {
nlist = read_idlist(id_file, idlist);
// sort and symmetrize the list
qsort(idlist, nlist, sizeof(int), compare_int);
int ninverses = 0;
LOOP(i, nlist) {
int id = idlist[i];
int invid = group->elements[id].inverse->id;
if(!bsearch(&invid, idlist, nlist, sizeof(int), compare_int)) {
idlist[nlist+ninverses] = invid;
ninverses++;
}
}
nlist += ninverses;
qsort(idlist, nlist, sizeof(int), compare_int);
} else {
// just list all elements which have inverses
nlist = 0;
LOOP(i, n) {
if(group->elements[i].inverse)
idlist[nlist++] = i;
}
}
// get generator matrices
INFO("make generators\n");
mat gen[3];
LOOP(i, 3) mat_init(gen[i], 3, QT_GAUSS_SQRT5);
generators_triangle_reflection_group_555_complex(gen, q[0], q[1], q[2], treal, timag);
// compute the traces of all elements in idlist
INFO("enumerate traces\n");
struct tracedata *traces;
nuniq = enumerate_coxeter_group_traces(group, gen, &traces, idlist, nlist, 1);
// compute eigenvalues out of traces
INFO("compute eigenvalues\n");
mps_context *solver = mps_context_new();
mps_monomial_poly *poly = mps_monomial_poly_new(solver, 3);
mps_context_set_output_prec(solver, 20); // relative precision
mps_context_set_output_goal(solver, MPS_OUTPUT_GOAL_APPROXIMATE);
double ev_real[3], ev_imag[3], ev_abs2[3];
double max_slope = 0;
int max_slope_id = 0;
LOOP(i, nuniq) {
solve_characteristic_polynomial_d(solver, poly,
gaussian_sqrt5_real(traces[i].tr),
gaussian_sqrt5_imag(traces[i].tr),
gaussian_sqrt5_real(traces[i].trinv),
gaussian_sqrt5_imag(traces[i].trinv),
ev_real, ev_imag);
LOOP(j, 3) ev_abs2[j] = ev_real[j]*ev_real[j] + ev_imag[j]*ev_imag[j];
if(fabs(ev_abs2[0]) < fabs(ev_abs2[1]))
SWAP(double, ev_abs2[0], ev_abs2[1]);
if(fabs(ev_abs2[1]) < fabs(ev_abs2[2]))
SWAP(double, ev_abs2[1], ev_abs2[2]);
if(fabs(ev_abs2[0]) < fabs(ev_abs2[1]))
SWAP(double, ev_abs2[0], ev_abs2[1]);
if(mode == MODE_TRACE_IDS) {
// we only want to record unordered pairs here
if(CMP(traces[i].tr, traces[i].trinv) >= 0)
printf("%d\n", traces[i].id);
continue;
}
if(log(ev_abs2[0]) < 1e-6) // we regard this as a finite order element
continue;
double slope = - log(ev_abs2[0]) / log(ev_abs2[2]);
if(slope > max_slope) {
max_slope = slope;
max_slope_id = traces[i].id;
}
if(mode == MODE_EIGENVALUES) {
printf("%d %f %f %f\n",
traces[i].id, log(ev_abs2[0])/2, log(ev_abs2[1])/2, log(ev_abs2[2])/2);
}
}
// output summary
coxeter_snprint(buf, sizeof(buf), &group->elements[max_slope_id]);
if(mode == MODE_SUMMARY) {
gmp_fprintf(stdout, "%f %f %d %d %d %f %s\n", mpq_get_d(treal), mpq_get_d(timag), n, nlist, nuniq, max_slope, buf);
} else if(mode == MODE_EIGENVALUES) {
gmp_fprintf(stderr, "q = %.2f + i*%.2f\tn = %d\tnlist = %d\tnuniq = %d\tMaximal slope: %f at %s\n", mpq_get_d(treal), mpq_get_d(timag), n, nlist, nuniq, max_slope, buf);
} else if(mode == MODE_TRACE_IDS) {
gmp_fprintf(stderr, "q = %.2f + i*%.2f\tElements: %d\tTraces: %d\n", mpq_get_d(treal), mpq_get_d(timag), n, nuniq);
}
// clean up
INFO("clean up\n");
// mps_monomial_poly_free(solver, MPS_POLYNOMIAL(poly));
mps_context_free(solver);
enumerate_tracedata_clear(traces, nuniq);
free(idlist);
LOOP(i, 3) mat_clear(gen[i]);
coxeter_clear(group);
mpq_clears(treal, timag, NULL);
return 0;
}

9
compute_picture.sh
View File

@ -1,9 +0,0 @@
#!/bin/bash
for i in $(seq -50 50); do
for j in $(seq 0 50); do
if [ $i -ne 0 ] || [ $j -ne 0 ]; then
IDLIST=output/idlist_13 ./complex_anosov summary 15124 1 1 1 $i/50 $j/50
fi
done
done

16
coxeter.c
View File

@ -5,8 +5,6 @@
#include "coxeter.h"
#define LOOP(i,n) for(int i = 0; i < (n); i++)
#define MIN(x,y) ((x) > (y) ? (y) : (x))
#define MAX(x,y) ((x) > (y) ? (x) : (y))
group_t *coxeter_init_triangle(int p, int q, int r, int nmax)
{
@ -187,17 +185,3 @@ void coxeter_clear(group_t *g)
free(g->lists);
free(g);
}
int coxeter_snprint(char *str, int buflen, groupelement_t *g)
{
int n = MIN(g->length, buflen-1); // number of characters without null byte
str[n] = 0;
int i = n-1;
while(g->parent) {
str[i--] = 'a' + g->letter;
g = g->parent;
}
return n;
}

2
coxeter.h
View File

@ -13,7 +13,6 @@ typedef struct _groupelement {
struct _groupelement *inverse;
struct _groupelement **left;
struct _groupelement **right;
struct _groupelement *conjugacy_class;
} groupelement_t;
typedef struct _group {
@ -27,6 +26,5 @@ typedef struct _group {
group_t *coxeter_init(int rank, int *coxeter_matrix, int nmax);
group_t *coxeter_init_triangle(int p, int q, int r, int nmax);
void coxeter_clear(group_t *g);
int coxeter_snprint(char *str, int buflen, groupelement_t *g);
#endif

231
enumerate.c
View File

@ -1,231 +0,0 @@
#include "enumerate.h"
#include <stdlib.h>
#define LOOP(i,n) for(int i = 0; i < (n); i++)
// #define INFO(msg, ...) fprintf(stderr, "[%10.3f] " msg, runtime(), ##__VA_ARGS__)
#define INFO(msg, ...)
static int compare_int(const void *a, const void *b)
{
const int *aa = a;
const int *bb = b;
return (*aa > *bb) - (*aa < *bb);
}
static int index_in_list(const int *list, int n, int key)
{
int *ptr = bsearch(&key, list, n, sizeof(int), compare_int);
if(ptr)
return ptr - list;
else
return -1;
}
// idlist is assumed to be sorted and symmetric
void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices, const int *idlist, int nlist)
{
TYPE type = GETTYPE(gen[0]->x[0]);
mat_workspace *ws;
// mark elements which we need to compute
int ncompute = 0;
LOOP(i, group->size) group->elements[i].need_to_compute = 0;
LOOP(i, nlist) {
groupelement_t *cur = &group->elements[idlist[i]];
while(cur && cur->need_to_compute == 0) {
cur->need_to_compute = 1;
ncompute++;
cur = cur->parent;
}
}
INFO("need to compute %d elements for %d elements in list\n", ncompute, nlist);
// compute them depth first
int maxlen = group->elements[group->size-1].length;
groupelement_t **stack = malloc((maxlen+1)*sizeof(groupelement_t));
int *letter_stack = malloc((maxlen+1)*sizeof(int));
mat *matrix_stack = malloc((maxlen+1)*sizeof(mat));
int level = 0;
LOOP(i, maxlen+1) mat_init(matrix_stack[i], 3, type);
ws = mat_workspace_init(3, type);
stack[0] = &group->elements[0];
mat_identity(matrix_stack[0]);
letter_stack[0] = 0;
while(level >= 0) {
int letter = letter_stack[level];
groupelement_t *next = stack[level]->left[letter];
if(next && next->length > level && next->need_to_compute) {
mat_multiply(ws, matrix_stack[level+1], gen[letter], matrix_stack[level]);
int id = next->id;
int listid = index_in_list(idlist, nlist, id);
if(listid != -1)
mat_copy(matrices[listid], matrix_stack[level+1]);
next->need_to_compute = 0;
stack[level+1] = next;
letter_stack[level+1] = 0;
level++;
} else {
letter = ++letter_stack[level];
while(letter >= group->rank) { // done with this level, go back down
level--;
if(level < 0)
break;
letter = ++letter_stack[level];
}
}
}
LOOP(i, maxlen+1) mat_clear(matrix_stack[i]);
mat_workspace_clear(ws);
}
static int compare_tracedata(const void *a_, const void *b_)
{
int d = 0;
struct tracedata **a = (struct tracedata **)a_;
struct tracedata **b = (struct tracedata **)b_;
d = CMP((*a)->tr,(*b)->tr);
if(d == 0) {
d = CMP((*a)->trinv, (*b)->trinv);
}
return d;
}
static int compare_tracedata_id(const void *a, const void *b)
{
int ida = (*(struct tracedata **)a)->id;
int idb = (*(struct tracedata **)b)->id;
return ida > idb ? 1 : ida < idb ? -1 : 0;
}
int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out, const int *idlist, int nlist, int unique)
{
TYPE t = GETTYPE(gen[0]->x[0]);
int nuniq;
mat *matrices = malloc(nlist*sizeof(mat));
struct tracedata *traces = malloc(nlist*sizeof(struct tracedata));
struct tracedata **distinct_traces = malloc(nlist*sizeof(struct tracedata*));
LOOP(i, nlist) mat_init(matrices[i], 3, t);
enumerate_coxeter_group(group, gen, matrices, idlist, nlist);
LOOP(i, nlist) {
INIT(traces[i].tr, t);
INIT(traces[i].trinv, t);
int inv_id_in_list = index_in_list(idlist, nlist, group->elements[idlist[i]].inverse->id);
mat_trace(traces[i].tr, matrices[i]);
mat_trace(traces[i].trinv, matrices[inv_id_in_list]);
traces[i].id = idlist[i];
distinct_traces[i] = &traces[i];
}
if(unique) {
qsort(distinct_traces, nlist, sizeof(struct tracedata*), compare_tracedata);
nuniq = 0;
LOOP(i, nlist) {
if(i == 0 || compare_tracedata(&distinct_traces[i], &distinct_traces[nuniq-1]) != 0) {
distinct_traces[nuniq] = distinct_traces[i];
nuniq++;
} else {
int oldlength = group->elements[distinct_traces[nuniq-1]->id].length;
int newlength = group->elements[distinct_traces[i]->id].length;
if(newlength < oldlength)
distinct_traces[nuniq-1]->id = distinct_traces[i]->id;
}
}
} else {
nuniq = nlist;
}
qsort(distinct_traces, nuniq, sizeof(struct tracedata*), compare_tracedata_id);
struct tracedata *unique_traces = malloc(nuniq*sizeof(struct tracedata));
LOOP(i, nuniq) {
INIT(unique_traces[i].tr, t);
INIT(unique_traces[i].trinv, t);
SET(unique_traces[i].tr, distinct_traces[i]->tr);
SET(unique_traces[i].trinv, distinct_traces[i]->trinv);
unique_traces[i].id = distinct_traces[i]->id;
}
LOOP(i, nlist) mat_clear(matrices[i]);
LOOP(i, nlist) {
CLEAR(traces[i].tr);
CLEAR(traces[i].trinv);
}
free(matrices);
free(traces);
free(distinct_traces);
*traces_out = unique_traces;
return nuniq;
}
void enumerate_tracedata_clear(struct tracedata *traces, int n)
{
LOOP(i, n) {
CLEAR(traces[i].tr);
CLEAR(traces[i].trinv);
}
free(traces);
}
// returns squares of absolute values
int solve_characteristic_polynomial_d(mps_context *solv, mps_monomial_poly *poly, double tr_real, double tr_imag, double trinv_real, double trinv_imag, double *eigenvalues_real, double *eigenvalues_imag)
{
// mpq_t coeff1, coeff2, zero;
cplx_t *roots;
double radii[3];
double *radii_p[3];
mps_boolean errors;
int result = 0;
/*
mpq_inits(coeff1, coeff2, zero, NULL);
mpq_set(coeff1, trinv);
mpq_sub(coeff2, zero, tr);
*/
mps_monomial_poly_set_coefficient_d(solv, poly, 0, -1, 0);
mps_monomial_poly_set_coefficient_d(solv, poly, 1, trinv_real, trinv_imag);
mps_monomial_poly_set_coefficient_d(solv, poly, 2, -tr_real, -tr_imag);
// mps_monomial_poly_set_coefficient_q(solv, poly, 1, coeff1, zero);
// mps_monomial_poly_set_coefficient_q(solv, poly, 2, coeff2, zero);
mps_monomial_poly_set_coefficient_d(solv, poly, 3, 1, 0);
mps_context_set_input_poly(solv, (mps_polynomial*)poly);
mps_mpsolve(solv);
roots = cplx_valloc(3);
for(int i = 0; i < 3; i++)
radii_p[i] = &(radii[i]);
mps_context_get_roots_d(solv, &roots, radii_p);
errors = mps_context_has_errors(solv);
if(errors) {
result = 1;
} else {
for(int i = 0; i < 3; i++) {
eigenvalues_real[i] = cplx_Re(roots[i]);
eigenvalues_imag[i] = cplx_Im(roots[i]);
// if(fabs(cplx_Im(roots[i])) > radii[i]) // non-real root
// result = 2;
}
}
cplx_vfree(roots);
// mpq_clears(coeff1, coeff2, zero, NULL);
return result;
}

22
enumerate.h
View File

@ -1,22 +0,0 @@
#ifndef ENUMERATE_H
#define ENUMERATE_H
#include "mat.h"
#include "coxeter.h"
#include <mps/mps.h>
struct tracedata {
int id;
NUMBER tr;
NUMBER trinv;
};
void enumerate_coxeter_group(group_t *group, mat *gen, mat *matrices, const int *idlist, int nlist);
int enumerate_coxeter_group_traces(group_t *group, mat *gen, struct tracedata **traces_out, const int *idlist, int nlist, int unique);
void enumerate_tracedata_clear(struct tracedata *traces, int n);
int solve_characteristic_polynomial_d(mps_context *solv, mps_monomial_poly *poly, double tr_real, double tr_imag, double trinv_real, double trinv_imag, double *eigenvalues_real, double *eigenvalues_imag);
#endif

0
old/enumerate_triangle_group.c → enumerate_triangle_group.c
View File

0
old/enumerate_triangle_group.h → enumerate_triangle_group.h
View File

148
generators.c
View File

@ -1,148 +0,0 @@
#include "generators.h"
#include "mat.h"
#define LOOP(i,n) for(int i = 0; i < (n); i++)
void generators_triangle_reflection_generic(mat *gen, NUMBER sqrho1, NUMBER sqrho2, NUMBER sqrho3, NUMBER t)
{
NUMBER tmp;
NUMBER tinv;
TYPE type = GETTYPE(sqrho1);
INIT(tmp, type);
INIT(tinv, type);
INV(tinv, t);
LOOP(i, 3) {
mat_init(gen[i], 3, type);
mat_zero(gen[i]);
LOOP(j, 3) {
SET_INT(*mat_ref(gen[i], j, j), i == j ? 1 : -1);
}
}
MUL(tmp, t, sqrho1);
NEG(*mat_ref(gen[2], 1, 2), tmp);
NEG(*mat_ref(gen[0], 2, 0), sqrho2);
NEG(*mat_ref(gen[1], 0, 1), sqrho3);
MUL(tmp, tinv, sqrho1);
NEG(*mat_ref(gen[1], 2, 1), tmp);
NEG(*mat_ref(gen[2], 0, 2), sqrho2);
NEG(*mat_ref(gen[0], 1, 0), sqrho3);
CLEAR(tmp);
CLEAR(tinv);
}
// warning: this is not compatible anymore!
int generators_triangle_reflection_group(mat *gen, int p1, int p2, int p3, int q1, int q2, int q3, mpq_t t)
{
int p[3] = {p1, p2, p3};
int q[3] = {q1, q2, q3};
TYPE type;
LOOP(i, 3)
if(p[i] < 2 || p[i] > 6)
return 0;
#ifdef QEXT
type = QT_TRIVIAL;
LOOP(i, 3)
if(p[i] == 5)
type = QT_SQRT5;
#else
LOOP(i, 3)
if(p[i] == 5)
return 0;
#endif
NUMBER rho[3];
LOOP(i, 3) INIT(rho[i], type);
// compute rho1, rho2, rho3
LOOP(i, 3) {
if(p[i] == 2 && q[i] == 1) {
SET_INT(rho[i], 0);
} else if(p[i] == 3 && q[i] == 1) {
SET_INT(rho[i], 1);
} else if(p[i] == 4 && q[i] == 1) {
SET_INT(rho[i], 2);
} else if(p[i] == 4 && q[i] == 2) {
SET_INT(rho[i], 0);
} else if(p[i] == 6 && q[i] == 1) {
SET_INT(rho[i], 3);
} else if(p[i] == 6 && q[i] == 2) {
SET_INT(rho[i], 1);
} else if(p[i] == 6 && q[i] == 3) {
SET_INT(rho[i], 0);
#ifdef QEXT
} else if(p[i] == 5 && q[i] == 1) {
mpq_set_si(rho[i]->a[0], 3, 2);
mpq_set_si(rho[i]->a[1], 1, 2);
} else if(p[i] == 5 && q[i] == 2) {
mpq_set_si(rho[i]->a[0], 3, 2);
mpq_set_si(rho[i]->a[1], -1, 2);
#endif
} else {
return 0;
}
}
NUMBER param;
INIT(param, type);
SET_Q(param, t);
generators_triangle_reflection_generic(gen, rho[0], rho[1], rho[2], param);
LOOP(i, 3) CLEAR(rho[i]);
CLEAR(param);
return 1;
}
int generators_triangle_reflection_group_555_complex(mat *gen, int q1, int q2, int q3, mpq_t treal, mpq_t timag)
{
int q[3] = {q1, q2, q3};
NUMBER sqrho[3];
LOOP(i, 3) INIT(sqrho[i], QT_GAUSS_SQRT5);
// compute sqrho1, sqrho2, sqrho3
LOOP(i, 3) {
if(q[i] == 1) {
// 2cos(pi/5) = 1/2 + 1/2*sqrt(5) = 1/2 + 7/12 a - 1/24 a^3
mpq_set_si(sqrho[i]->a[0], 1, 2);
mpq_set_si(sqrho[i]->a[1], 7, 12);
mpq_set_si(sqrho[i]->a[2], 0, 1);
mpq_set_si(sqrho[i]->a[3], -1, 24);
} else if(q[i] == 2) {
// 2cos(pi/5) = -1/2 + 1/2*sqrt(5)
mpq_set_si(sqrho[i]->a[0], -1, 2);
mpq_set_si(sqrho[i]->a[1], 7, 12);
mpq_set_si(sqrho[i]->a[2], 0, 1);
mpq_set_si(sqrho[i]->a[3], -1, 24);
} else {
return 0;
}
}
NUMBER param;
INIT(param, QT_GAUSS_SQRT5);
mpq_set(param->a[0], treal);
mpq_set_si(param->a[1], -1, 6);
mpq_mul(param->a[1], param->a[1], timag);
mpq_set_si(param->a[2], 0, 1);
mpq_set_si(param->a[3], 1, 12);
mpq_mul(param->a[3], param->a[3], timag);
generators_triangle_reflection_generic(gen, sqrho[0], sqrho[1], sqrho[2], param);
LOOP(i, 3) CLEAR(sqrho[i]);
CLEAR(param);
return 1;
}

13
generators.h
View File

@ -1,13 +0,0 @@
#ifndef GENERATORS_H
#define GENERATORS_H
#include "mat.h"
void generators_triangle_reflection_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER q);
int generators_triangle_reflection_group(mat *gen, int p1, int p2, int p3, int q1, int q2, int q3, mpq_t q);
#ifdef QEXT
int generators_triangle_reflection_group_555_complex(mat *gen, int q1, int q2, int q3, mpq_t real, mpq_t imag);
#endif
#endif

67
make_picture.py
View File

@ -1,67 +0,0 @@
#!/usr/bin/python
import sys
import math
def sharpen(x):
alpha = 1
y = abs(2*x-1)**alpha # between 0 and 1
if x > 0.499 and x < 0.501:
return 0.5
elif x > 0.5:
return (y+1)/2
else:
return (1-y)/2
f = open(sys.argv[1])
lines = [x.split() for x in f]
f.close()
res1 = int(sys.argv[2]) # 400 pixel per unit
res2 = int(sys.argv[3]) # 50 pixel in picture (acutally one quadrant)
data = {(round(float(l[0])*res1), round(float(l[1])*res1)) : float(l[5]) for l in lines}
data[(0,0)] = 2.0
print("P3")
print("1000 1000")
print("255")
for i in range (-500,500):
for j in range(-500,500):
x = j/500*res2
y = i/500*res2
if y < 0:
y = -y
x0 = math.floor(x)
y0 = math.floor(y)
x1 = math.ceil(x)
y1 = math.ceil(y)
tx = sharpen(x-x0)
ty = sharpen(y-y0)
if not (x0,y0) in data or not (x0,y1) in data or not (x1,y0) in data or not (x1,y1) in data:
value = 0
else:
value = 0.0
value += data[(x0,y0)]*(1-tx)*(1-ty)
value += data[(x0,y1)]*(1-tx)*ty
value += data[(x1,y0)]*tx*(1-ty)
value += data[(x1,y1)]*tx*ty
value -= 1
value = 0 if value < 0 else 1 if value > 1 else value
r = round(255*math.sqrt(value))
g = round(255*(value)**3)
b = round(255*math.sin(2*math.pi*(value)))
r = 0 if r < 0 else 255 if r > 255 else r
g = 0 if g < 0 else 255 if g > 255 else g
b = 0 if b < 0 else 255 if b > 255 else b
print("%d %d %d" % (r,g,b))
#print(data)

4
mat.h
View File

@ -24,8 +24,7 @@
#define SUB qext_sub
#define NEG qext_neg
#define MUL qext_mul
#define DIV qext_div
#define INV qext_inv
// #define DIV qext_div
#define CMP qext_cmp
#define PRINT qext_print
#define SNPRINT qext_snprint
@ -47,7 +46,6 @@
#define NEG mpq_neg
#define MUL mpq_mul
#define DIV mpq_div
#define INV mpq_inv
#define PRINT(x) gmp_printf("%Qd", x)
#define SNPRINT(out, size, x) gmp_snprintf(out, size, "%Qd", x)
#define TYPE int

0
parallelization/parallel.c → parallel.c
View File

0
parallelization/parallel.h → parallel.h
View File

102
qext.c
View File

@ -1,7 +1,5 @@
#include "qext.h"
#include <assert.h>
int qext_rank;
mpq_t *qext_coefficient;
@ -11,20 +9,13 @@ mpq_t *qext_coefficient;
#define LOOP(i,n) for(int i = 0; i < (n); i++)
#define RANK(x) ((x)->type->rank)
struct qext_type *QT_TRIVIAL = &(struct qext_type) {
.rank = 1,
.integer_coeffs = (const int[]) {-1, 1}
};
const int qext_type_trivial_coeffs[] = {-1, 1};
struct qext_type qext_type_trivial_v = {1, qext_type_trivial_coeffs, NULL};
struct qext_type *QT_TRIVIAL = &qext_type_trivial_v;
struct qext_type *QT_SQRT5 = &(struct qext_type) {
.rank = 2,
.integer_coeffs = (const int[]) {-5, 0, 1}
};
struct qext_type *QT_GAUSS_SQRT5 = &(struct qext_type) {
.rank = 4,
.integer_coeffs = (const int[]) {36, 0, -8, 0, 1}
};
const int qext_type_sqrt5_coeffs[] = {-5, 0, 1};
struct qext_type qext_type_sqrt5_v = {2, qext_type_sqrt5_coeffs, NULL};
struct qext_type *QT_SQRT5 = &qext_type_sqrt5_v;
static void qext_init_type(struct qext_type *type)
{
@ -170,86 +161,7 @@ void qext_mul(qext_number out, qext_number x, qext_number y)
mpq_clear(tmp);
}
void qext_inv(qext_number out, qext_number in)
{
qext_number one;
qext_init(one, in->type);
qext_set_int(one, 1);
qext_div(out, one, in);
qext_clear(one);
}
void qext_div(qext_number out, qext_number x, qext_number y)
{
int rk = RANK(x);
struct qext_type *type = x->type;
mpq_t tmp;
mpq_t result[rk];
mpq_t matrix[rk*rk];
mpq_init(tmp);
LOOP(i, rk) mpq_init(result[i]);
LOOP(i, rk*rk) mpq_init(matrix[i]);
// initialize a matrix which represents multiplication by y
// that means the i-th column is y*a^i
LOOP(i, rk)
mpq_set(matrix[i*rk], y->a[i]);
LOOP(j, rk-1) { // columns
mpq_mul(matrix[j+1], matrix[(rk-1)*rk+j], type->coeffs[0]);
LOOP(i, rk-1) { // rows
mpq_mul(matrix[(i+1)*rk+(j+1)], matrix[(rk-1)*rk+j], type->coeffs[i+1]);
mpq_add(matrix[(i+1)*rk+(j+1)], matrix[(i+1)*rk+(j+1)], matrix[i*rk+j]);
}
}
// initialize result as x
LOOP(i, rk) mpq_set(result[i], x->a[i]);
// use Gaussian elimination to solve the system matrix * out = result
LOOP(j, rk)
{
// find nonzero entry
int row = j;
while(row < rk && mpq_sgn(matrix[row*rk+j]) == 0) row++;
// if the entire column is zero, the matrix was not invertible
// this could happen if the polynomial is not irreducible / we're not in a field
assert(row != rk);
// permute rows
if(row != j) {
mpq_swap(result[row], result[j]);
LOOP(k, rk) {
mpq_swap(matrix[row*rk+k], matrix[j*rk+k]);
}
}
// normalize
LOOP(k, rk)
if(k != j)
mpq_div(matrix[j*rk+k], matrix[j*rk+k], matrix[j*rk+j]);
mpq_div(result[j], result[j], matrix[j*rk+j]);
mpq_set_ui(matrix[j*rk+j], 1, 1);
// subtract
LOOP(i, rk) {
if(i == j)
continue;
mpq_mul(tmp, matrix[i*rk+j], result[j]);
mpq_sub(result[i], result[i], tmp);
for(int k = j+1; k < rk; k++) {
mpq_mul(tmp, matrix[i*rk+j], matrix[j*rk+k]);
mpq_sub(matrix[i*rk+k], matrix[i*rk+k], tmp);
}
mpq_set_ui(matrix[i*rk+j], 0, 1); // it isn't strictly necessary to do this as we won't use this entry again, but helps debugging
}
}
LOOP(i, rk) mpq_set(out->a[i], result[i]);
mpq_clear(tmp);
LOOP(i, rk) mpq_clear(result[i]);
LOOP(i, rk*rk) mpq_clear(matrix[i]);
// todo: implement this by solving a linear system
}

6
qext.h
View File

@ -5,8 +5,8 @@
struct qext_type {
int rank;
const int *integer_coeffs; // actually the rank+1 coefficients of the polynomial
mpq_t *coeffs; // components of a^rank
const int *integer_coeffs;
mpq_t *coeffs;
};
struct qext_number_internal {
@ -18,7 +18,6 @@ typedef struct qext_number_internal qext_number[1];
extern struct qext_type *QT_TRIVIAL;
extern struct qext_type *QT_SQRT5;
extern struct qext_type *QT_GAUSS_SQRT5;
struct qext_type *qext_newtype(int rank, const int *coeffs);
void qext_init(qext_number x, struct qext_type *type);
@ -30,7 +29,6 @@ void qext_add(qext_number result, qext_number x, qext_number y);
void qext_sub(qext_number result, qext_number x, qext_number y);
void qext_neg(qext_number result, qext_number x);
void qext_mul(qext_number out, qext_number x, qext_number y);
void qext_inv(qext_number out, qext_number in);
void qext_div(qext_number out, qext_number x, qext_number y);
int qext_cmp(qext_number x, qext_number y);
void qext_print(qext_number x);

0
old/singular_values.c → singular_values.c
View File

0
old/singular_values_barbot.c → singular_values_barbot.c
View File