Compare commits
No commits in common. "429f0890d6b02a7aed8d42c763048e5de47441f0" and "3309c37955053930947b0d087e80d3552c146096" have entirely different histories.
429f0890d6
...
3309c37955
20
.gitignore
vendored
20
.gitignore
vendored
@ -1,7 +1,17 @@
|
|||||||
*.o
|
*.o
|
||||||
.\#*
|
triangle_group/singular_values
|
||||||
\#*\#
|
.#*
|
||||||
gmon.out
|
singular_values
|
||||||
core
|
singular_values_mpi
|
||||||
|
singular_values_barbot
|
||||||
output/
|
output/
|
||||||
complex_anosov
|
special_element
|
||||||
|
max_slope_picture/generate
|
||||||
|
convert
|
||||||
|
billiard_words
|
||||||
|
*.pnm
|
||||||
|
*.png
|
||||||
|
*.hi
|
||||||
|
gmon.out
|
||||||
|
restart
|
||||||
|
core
|
||||||
|
|||||||
54
Makefile
54
Makefile
@ -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 -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=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL
|
||||||
#SPECIAL_OPTIONS=
|
#SPECIAL_OPTIONS=
|
||||||
|
|
||||||
OPTIONS=-I../mps/include -L../mps/lib -pthread -m64 -std=gnu99 -D_GNU_SOURCE $(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
|
convert: convert.hs
|
||||||
$(CC) $(OPTIONS) -o complex_anosov -lm complex_anosov.o mat.o coxeter.o enumerate.o generators.o qext.o -lgmp -lmps
|
ghc --make -dynamic convert.hs
|
||||||
|
|
||||||
complex_anosov.o: complex_anosov.c $(HEADERS)
|
billiard_words: billiard_words.hs
|
||||||
gcc $(OPTIONS) -c complex_anosov.c
|
ghc --make -dynamic billiard_words.hs
|
||||||
|
|
||||||
enumerate.o: enumerate.c $(HEADERS)
|
singular_values: singular_values.o coxeter.o mat.o enumerate_triangle_group.o parallel.o
|
||||||
gcc $(OPTIONS) -c enumerate.c
|
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)
|
singular_values_barbot: singular_values_barbot.o coxeter.o mat.o enumerate_triangle_group.o parallel.o qext.o
|
||||||
gcc $(OPTIONS) -c generators.c
|
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)
|
coxeter.o: coxeter.c $(HEADERS)
|
||||||
gcc $(OPTIONS) -c coxeter.c
|
gcc $(OPTIONS) -c coxeter.c
|
||||||
@ -29,8 +52,11 @@ coxeter.o: coxeter.c $(HEADERS)
|
|||||||
mat.o: mat.c $(HEADERS)
|
mat.o: mat.c $(HEADERS)
|
||||||
gcc $(OPTIONS) -c mat.c
|
gcc $(OPTIONS) -c mat.c
|
||||||
|
|
||||||
|
parallel.o: parallel.c $(HEADERS)
|
||||||
|
gcc $(OPTIONS) -c parallel.c
|
||||||
|
|
||||||
qext.o: qext.c $(HEADERS)
|
qext.o: qext.c $(HEADERS)
|
||||||
gcc $(OPTIONS) -c qext.c
|
gcc $(OPTIONS) -c qext.c
|
||||||
|
|
||||||
clean:
|
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
304
complex_anosov.c
@ -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;
|
|
||||||
}
|
|
||||||
@ -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
16
coxeter.c
@ -5,8 +5,6 @@
|
|||||||
#include "coxeter.h"
|
#include "coxeter.h"
|
||||||
|
|
||||||
#define LOOP(i,n) for(int i = 0; i < (n); i++)
|
#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)
|
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->lists);
|
||||||
free(g);
|
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;
|
|
||||||
}
|
|
||||||
|
|||||||
@ -13,7 +13,6 @@ typedef struct _groupelement {
|
|||||||
struct _groupelement *inverse;
|
struct _groupelement *inverse;
|
||||||
struct _groupelement **left;
|
struct _groupelement **left;
|
||||||
struct _groupelement **right;
|
struct _groupelement **right;
|
||||||
struct _groupelement *conjugacy_class;
|
|
||||||
} groupelement_t;
|
} groupelement_t;
|
||||||
|
|
||||||
typedef struct _group {
|
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(int rank, int *coxeter_matrix, int nmax);
|
||||||
group_t *coxeter_init_triangle(int p, int q, int r, int nmax);
|
group_t *coxeter_init_triangle(int p, int q, int r, int nmax);
|
||||||
void coxeter_clear(group_t *g);
|
void coxeter_clear(group_t *g);
|
||||||
int coxeter_snprint(char *str, int buflen, groupelement_t *g);
|
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
|||||||
231
enumerate.c
231
enumerate.c
@ -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
22
enumerate.h
@ -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
|
|
||||||
148
generators.c
148
generators.c
@ -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
13
generators.h
@ -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
|
|
||||||
@ -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
4
mat.h
@ -24,8 +24,7 @@
|
|||||||
#define SUB qext_sub
|
#define SUB qext_sub
|
||||||
#define NEG qext_neg
|
#define NEG qext_neg
|
||||||
#define MUL qext_mul
|
#define MUL qext_mul
|
||||||
#define DIV qext_div
|
// #define DIV qext_div
|
||||||
#define INV qext_inv
|
|
||||||
#define CMP qext_cmp
|
#define CMP qext_cmp
|
||||||
#define PRINT qext_print
|
#define PRINT qext_print
|
||||||
#define SNPRINT qext_snprint
|
#define SNPRINT qext_snprint
|
||||||
@ -47,7 +46,6 @@
|
|||||||
#define NEG mpq_neg
|
#define NEG mpq_neg
|
||||||
#define MUL mpq_mul
|
#define MUL mpq_mul
|
||||||
#define DIV mpq_div
|
#define DIV mpq_div
|
||||||
#define INV mpq_inv
|
|
||||||
#define PRINT(x) gmp_printf("%Qd", x)
|
#define PRINT(x) gmp_printf("%Qd", x)
|
||||||
#define SNPRINT(out, size, x) gmp_snprintf(out, size, "%Qd", x)
|
#define SNPRINT(out, size, x) gmp_snprintf(out, size, "%Qd", x)
|
||||||
#define TYPE int
|
#define TYPE int
|
||||||
|
|||||||
102
qext.c
102
qext.c
@ -1,7 +1,5 @@
|
|||||||
#include "qext.h"
|
#include "qext.h"
|
||||||
|
|
||||||
#include <assert.h>
|
|
||||||
|
|
||||||
int qext_rank;
|
int qext_rank;
|
||||||
mpq_t *qext_coefficient;
|
mpq_t *qext_coefficient;
|
||||||
|
|
||||||
@ -11,20 +9,13 @@ mpq_t *qext_coefficient;
|
|||||||
#define LOOP(i,n) for(int i = 0; i < (n); i++)
|
#define LOOP(i,n) for(int i = 0; i < (n); i++)
|
||||||
#define RANK(x) ((x)->type->rank)
|
#define RANK(x) ((x)->type->rank)
|
||||||
|
|
||||||
struct qext_type *QT_TRIVIAL = &(struct qext_type) {
|
const int qext_type_trivial_coeffs[] = {-1, 1};
|
||||||
.rank = 1,
|
struct qext_type qext_type_trivial_v = {1, qext_type_trivial_coeffs, NULL};
|
||||||
.integer_coeffs = (const int[]) {-1, 1}
|
struct qext_type *QT_TRIVIAL = &qext_type_trivial_v;
|
||||||
};
|
|
||||||
|
|
||||||
struct qext_type *QT_SQRT5 = &(struct qext_type) {
|
const int qext_type_sqrt5_coeffs[] = {-5, 0, 1};
|
||||||
.rank = 2,
|
struct qext_type qext_type_sqrt5_v = {2, qext_type_sqrt5_coeffs, NULL};
|
||||||
.integer_coeffs = (const int[]) {-5, 0, 1}
|
struct qext_type *QT_SQRT5 = &qext_type_sqrt5_v;
|
||||||
};
|
|
||||||
|
|
||||||
struct qext_type *QT_GAUSS_SQRT5 = &(struct qext_type) {
|
|
||||||
.rank = 4,
|
|
||||||
.integer_coeffs = (const int[]) {36, 0, -8, 0, 1}
|
|
||||||
};
|
|
||||||
|
|
||||||
static void qext_init_type(struct qext_type *type)
|
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);
|
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)
|
void qext_div(qext_number out, qext_number x, qext_number y)
|
||||||
{
|
{
|
||||||
int rk = RANK(x);
|
// todo: implement this by solving a linear system
|
||||||
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]);
|
|
||||||
}
|
}
|
||||||
|
|||||||
6
qext.h
6
qext.h
@ -5,8 +5,8 @@
|
|||||||
|
|
||||||
struct qext_type {
|
struct qext_type {
|
||||||
int rank;
|
int rank;
|
||||||
const int *integer_coeffs; // actually the rank+1 coefficients of the polynomial
|
const int *integer_coeffs;
|
||||||
mpq_t *coeffs; // components of a^rank
|
mpq_t *coeffs;
|
||||||
};
|
};
|
||||||
|
|
||||||
struct qext_number_internal {
|
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_TRIVIAL;
|
||||||
extern struct qext_type *QT_SQRT5;
|
extern struct qext_type *QT_SQRT5;
|
||||||
extern struct qext_type *QT_GAUSS_SQRT5;
|
|
||||||
|
|
||||||
struct qext_type *qext_newtype(int rank, const int *coeffs);
|
struct qext_type *qext_newtype(int rank, const int *coeffs);
|
||||||
void qext_init(qext_number x, struct qext_type *type);
|
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_sub(qext_number result, qext_number x, qext_number y);
|
||||||
void qext_neg(qext_number result, qext_number x);
|
void qext_neg(qext_number result, qext_number x);
|
||||||
void qext_mul(qext_number out, qext_number x, qext_number y);
|
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);
|
void qext_div(qext_number out, qext_number x, qext_number y);
|
||||||
int qext_cmp(qext_number x, qext_number y);
|
int qext_cmp(qext_number x, qext_number y);
|
||||||
void qext_print(qext_number x);
|
void qext_print(qext_number x);
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user