From c566a8741e74babc0ce4805c2285a1f7afdea0d1 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 1 Jul 2020 21:26:42 -0500 Subject: [PATCH] switch to rational numbers + add variable type matrix library --- Makefile | 15 ++- mat.c | 158 +++++++++++++++++++++++++ mat.h | 66 +++++++++++ singular_values.c | 281 +++++++++++++++++++------------------------- singular_values.plt | 30 ++++- 5 files changed, 379 insertions(+), 171 deletions(-) create mode 100644 mat.c create mode 100644 mat.h diff --git a/Makefile b/Makefile index fa862ff..f71aaf0 100644 --- a/Makefile +++ b/Makefile @@ -1,16 +1,16 @@ -HEADERS=triangle.h linalg.h queue.h +HEADERS=triangle.h linalg.h queue.h mat.h #SPECIAL_OPTIONS=-O0 -g -D_DEBUG -#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline -SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline +SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline +#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline #SPECIAL_OPTIONS= OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS) all: singular_values -singular_values: singular_values.o triangle.o linalg.o - gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o -lm -lgsl -lcblas +singular_values: singular_values.o triangle.o linalg.o mat.o + gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o mat.o -lm -lgsl -lcblas -lgmp singular_values.o: singular_values.c $(HEADERS) gcc $(OPTIONS) -c singular_values.c @@ -21,5 +21,8 @@ linalg.o: linalg.c $(HEADERS) triangle.o: triangle.c $(HEADERS) gcc $(OPTIONS) -c triangle.c +mat.o: mat.c $(HEADERS) + gcc $(OPTIONS) -c mat.c + clean: - rm -f singular_values triangle.o linalg.o singular_values.o + rm -f singular_values triangle.o linalg.o singular_values.o mat.o diff --git a/mat.c b/mat.c new file mode 100644 index 0000000..7ed0b5a --- /dev/null +++ b/mat.c @@ -0,0 +1,158 @@ +#include "mat.h" + +mat_workspace *mat_workspace_init(int n) +{ + mat_workspace *ws = (mat_workspace*)malloc(sizeof(mat_workspace)); + mat_init(ws->tmp_mat, n); + INIT(ws->tmp_num); + INIT(ws->tmp_num2); + return ws; +} + +void mat_workspace_clear(mat_workspace *ws) +{ + mat_clear(ws->tmp_mat); + CLEAR(ws->tmp_num); + CLEAR(ws->tmp_num2); + free(ws); +} + +void mat_init(mat m, int n) +{ + m->n = n; + m->x = malloc(n*n*sizeof(NUMBER)); + LOOP(i,n) LOOP(j,n) INIT(m->x[i+j*n]); +} + +void mat_get(NUMBER out, mat m, int i, int j) +{ + SET(out, M(m,i,j)); +} + +void mat_set(mat m, int i, int j, NUMBER x) +{ + SET(M(m,i,j), x); +} + +NUMBER *mat_ref(mat m, int i, int j) +{ + return &M(m,i,j); +} + +void mat_zero(mat m) +{ + LOOP(i,m->n) LOOP(j,m->n) SET_ZERO(M(m,i,j)); +} + +void mat_identity(mat m) +{ + LOOP(i,m->n) LOOP(j,m->n) { + if(i == j) + SET_ONE(M(m,i,j)); + else + SET_ZERO(M(m,i,j)); + } +} + +void mat_copy(mat to, mat from) +{ + LOOP(i,from->n) LOOP(j,from->n) SET(M(to,i,j), M(from,i,j)); +} + +void mat_clear(mat m) +{ + LOOP(i,m->n) LOOP(j,m->n) CLEAR(M(m,i,j)); + free(m->x); +} + +int mat_same(mat m1, mat m2) +{ + return m1 == m2; +} + +static void mat_multiply_outofplace(mat_workspace *ws, mat out, mat in1, mat in2) +{ + int n = out->n; + NUMBER *tmp = &(ws->tmp_num); + + LOOP(i,n) LOOP(j,n) { + SET_ZERO(M(out,i,j)); + LOOP(k,n) { + MULTIPLY(*tmp, M(in1,i,k), M(in2,k,j)); + ADD(M(out,i,j), M(out,i,j), *tmp); + } + } +} + +void mat_multiply(mat_workspace *ws, mat out, mat in1, mat in2) +{ + if(mat_same(out, in1) || mat_same(out, in2)) { + mat_multiply_outofplace(ws, ws->tmp_mat, in1, in2); + mat_copy(out, ws->tmp_mat); + } else { + mat_multiply_outofplace(ws, out, in1, in2); + } +} + +void mat_det(mat_workspace *ws, NUMBER out, mat in) +{ +// let's just assume n = 3 for the moment + NUMBER *tmp = &(ws->tmp_num); + int n = 3; + + SET_ZERO(out); + + LOOP(i,n) { + MULTIPLY(*tmp, M(in,0,i), M(in,1,(i+1)%3)); + MULTIPLY(*tmp, *tmp, M(in,2,(i+2)%3)); + ADD(out, out, *tmp); + + MULTIPLY(*tmp, M(in,0,(i+2)%3), M(in,1,(i+1)%3)); + MULTIPLY(*tmp, *tmp, M(in,2,i)); + SUB(out, out, *tmp); + } +} + +static void mat_pseudoinverse_outofplace(mat_workspace *ws, mat out, mat in) +{ + // let's just assume n = 3 for the moment + NUMBER *tmp = &(ws->tmp_num); + NUMBER *tmp2 = &(ws->tmp_num2); + int n = 3; + + LOOP(i,n) LOOP(j,n) { + MULTIPLY(*tmp, M(in,(i+1)%3,(j+1)%3), M(in,(i+2)%3,(j+2)%3)); + MULTIPLY(*tmp2, M(in,(i+1)%3,(j+2)%3), M(in,(i+2)%3,(j+1)%3)); + SUB(M(out,j,i), *tmp, *tmp2); + } +} + +void mat_pseudoinverse(mat_workspace *ws, mat out, mat in) +{ + if(mat_same(out, in)) { + mat_pseudoinverse_outofplace(ws, ws->tmp_mat, in); + mat_copy(out, ws->tmp_mat); + } else { + mat_pseudoinverse_outofplace(ws, out, in); + } +} + +void mat_trace(NUMBER out, mat in) +{ + SET_ZERO(out); + ADD(out, out, M(in,0,0)); + ADD(out, out, M(in,1,1)); + ADD(out, out, M(in,2,2)); +} + +void mat_print(mat in) +{ + int n = in->n; + + LOOP(i,n) { + LOOP(j,n) { + PRINT(M(in,i,j)); + fputc(j == n - 1 ? '\n' : ' ', stdout); + } + } +} diff --git a/mat.h b/mat.h new file mode 100644 index 0000000..06b067b --- /dev/null +++ b/mat.h @@ -0,0 +1,66 @@ +#ifndef MAT_H +#define MAT_H + +#include +#include + +#define LOOP(i,n) for(int i = 0; i < (n); i++) + +// library for matrix computations in variable rings (based on GMP types) + +/* + needed features: + x multiply matrices + - inverse + - pseudoinverse + x set + - eigenvalues + */ + +#define NUMBER mpq_t +#define INIT mpq_init +#define CLEAR mpq_clear +#define SET mpq_set +#define SET_ZERO(x) mpq_set_ui(x,0,1) +#define SET_ONE(x) mpq_set_ui(x,1,1) +#define ADD mpq_add +#define SUB mpq_sub +#define MULTIPLY mpq_mul +#define DIV mpq_div +#define PRINT(x) gmp_printf("%Qd", x) + +#define M(m,i,j) ((m)->x[(i)+(m)->n*(j)]) + +struct _mat{ + int n; + NUMBER *x; +} ; + +typedef struct _mat mat[1]; + +typedef struct _mat_workspace { + mat tmp_mat; + NUMBER tmp_num; + NUMBER tmp_num2; +} mat_workspace; + +mat_workspace *mat_workspace_init(int n); +void mat_workspace_clear(mat_workspace *ws); +void mat_init(mat m, int n); +void mat_get(NUMBER out, mat m, int i, int j); +void mat_set(mat m, int i, int j, NUMBER x); +NUMBER *mat_ref(mat m, int i, int j); +void mat_zero(mat m); +void mat_identity(mat m); +void mat_copy(mat to, mat from); +void mat_clear(mat m); +int mat_same(mat m1, mat m2); +static void mat_multiply_outofplace(mat_workspace *ws, mat out, mat in1, mat in2); +void mat_multiply(mat_workspace *ws, mat out, mat in1, mat in2); +void mat_det(mat_workspace *ws, NUMBER out, mat in); +static void mat_pseudoinverse_outofplace(mat_workspace *ws, mat out, mat in); +void mat_pseudoinverse(mat_workspace *ws, mat out, mat in); +void mat_trace(NUMBER out, mat in); +void mat_print(mat in); + +#endif diff --git a/singular_values.c b/singular_values.c index 44e9f76..bc06c61 100644 --- a/singular_values.c +++ b/singular_values.c @@ -1,60 +1,69 @@ #include "triangle.h" #include "linalg.h" +#include "mat.h" +#include + +//#define MAX_ELEMENTS 2800000 +//#define MAX_ELEMENTS 720000 #define MAX_ELEMENTS 14000 //#define DRAW_PICTURE 1 -void initialize_triangle_generators(workspace_t *ws, gsl_matrix **gen, double a1, double a2, double a3, double s, double t) +void mpq_quartic(mpq_t out, mpq_t in, int a, int b, int c, int d, int e) { - gsl_matrix **tmp; - tmp = getTempMatrices(ws, 6); + mpq_t tmp; + mpq_init(tmp); - double rho1sqrt = sqrt(s*(s+2*cos(2*M_PI*a1))+1); - double rho2sqrt = sqrt(s*(s+2*cos(2*M_PI*a2))+1); - double rho3sqrt = sqrt(s*(s+2*cos(2*M_PI*a3))+1); + mpq_set_si(out, a, 1); + mpq_mul(out, out, in); + mpq_set_si(tmp, b, 1); + mpq_add(out, out, tmp); + mpq_mul(out, out, in); + mpq_set_si(tmp, c, 1); + mpq_add(out, out, tmp); + mpq_mul(out, out, in); + mpq_set_si(tmp, d, 1); + mpq_add(out, out, tmp); + mpq_mul(out, out, in); + mpq_set_si(tmp, e, 1); + mpq_add(out, out, tmp); - for(int i = 0; i < 6; i++) - gsl_matrix_set_zero(tmp[i]); + mpq_clear(tmp); +} - gsl_matrix_set(tmp[0], 0, 0, 1); - gsl_matrix_set(tmp[0], 1, 0, -rho3sqrt/t); - gsl_matrix_set(tmp[0], 2, 0, -rho2sqrt*t); - gsl_matrix_set(tmp[0], 1, 1, -1); - gsl_matrix_set(tmp[0], 2, 2, -1); +void initialize_triangle_generators(mat_workspace *ws, mat *gen, mpq_t s, mpq_t t) +{ + mpq_set_ui(*mat_ref(gen[0], 0, 0), 0, 1); + mpq_set_ui(*mat_ref(gen[0], 0, 1), 0, 1); + mpq_set_ui(*mat_ref(gen[0], 0, 2), 1, 1); + mpq_set_ui(*mat_ref(gen[0], 1, 0), 1, 1); + mpq_set_ui(*mat_ref(gen[0], 1, 1), 0, 1); + mpq_set_ui(*mat_ref(gen[0], 1, 2), 0, 1); + mpq_set_ui(*mat_ref(gen[0], 2, 0), 0, 1); + mpq_set_ui(*mat_ref(gen[0], 2, 1), 1, 1); + mpq_set_ui(*mat_ref(gen[0], 2, 2), 0, 1); - gsl_matrix_set(tmp[1], 0, 0, -1); - gsl_matrix_set(tmp[1], 0, 1, -rho3sqrt*t); - gsl_matrix_set(tmp[1], 1, 1, 1); - gsl_matrix_set(tmp[1], 2, 1, -rho1sqrt/t); - gsl_matrix_set(tmp[1], 2, 2, -1); + mpq_set_ui(*mat_ref(gen[1], 0, 0), 1, 1); + mpq_set_ui(*mat_ref(gen[1], 1, 0), 0, 1); + mpq_set_ui(*mat_ref(gen[1], 2, 0), 0, 1); + mpq_quartic(*mat_ref(gen[1], 0, 1), t, 0, 0, 1, -1, 2); + mpq_quartic(*mat_ref(gen[1], 1, 1), t, 0, 0, -1, 2, -2); + mpq_quartic(*mat_ref(gen[1], 2, 1), t, 0, 0, 1, -3, 3); + mpq_quartic(*mat_ref(gen[1], 0, 2), t, 0, 0, 1, 0, 3); + mpq_quartic(*mat_ref(gen[1], 1, 2), t, 0, 0, -1, 1, -1); + mpq_quartic(*mat_ref(gen[1], 2, 2), t, 0, 0, 1, -2, 1); - gsl_matrix_set(tmp[2], 0, 0, -1); - gsl_matrix_set(tmp[2], 1, 1, -1); - gsl_matrix_set(tmp[2], 0, 2, -rho2sqrt/t); - gsl_matrix_set(tmp[2], 1, 2, -rho1sqrt*t); - gsl_matrix_set(tmp[2], 2, 2, 1); + mat_pseudoinverse(ws, gen[3], gen[0]); // p^{-1} + mat_pseudoinverse(ws, gen[4], gen[1]); // q^{-1} + mat_multiply(ws, gen[2], gen[4], gen[3]); // r = q^{-1}p^{-1} + mat_pseudoinverse(ws, gen[5], gen[2]); - gsl_matrix_set(tmp[3], 0, 0, 1); - gsl_matrix_set(tmp[3], 1, 1, s); - gsl_matrix_set(tmp[3], 2, 2, 1/s); - - gsl_matrix_set(tmp[4], 0, 0, 1/s); - gsl_matrix_set(tmp[4], 1, 1, 1); - gsl_matrix_set(tmp[4], 2, 2, s); - - gsl_matrix_set(tmp[5], 0, 0, s); - gsl_matrix_set(tmp[5], 1, 1, 1/s); - gsl_matrix_set(tmp[5], 2, 2, 1); - - multiply_many(ws, gen[0], 3, tmp[2], tmp[3], tmp[1]); - multiply_many(ws, gen[1], 3, tmp[0], tmp[4], tmp[2]); - multiply_many(ws, gen[2], 3, tmp[1], tmp[5], tmp[0]); - - invert(gen[0], gen[3], ws); - invert(gen[1], gen[4], ws); - invert(gen[2], gen[5], ws); - - releaseTempMatrices(ws, 6); +// mat_print(gen[0]); +// mat_print(gen[1]); +// mat_print(gen[2]); +// mat_print(gen[3]); +// mat_print(gen[4]); +// mat_print(gen[5]); } char *print_word(groupelement_t *g, char *str) @@ -73,138 +82,92 @@ char *print_word(groupelement_t *g, char *str) int main(int argc, char *argv[]) { groupelement_t *group; - workspace_t *ws; - gsl_matrix **matrices; - gsl_matrix *gen[6]; - gsl_matrix *tmp, *tmp2, *coxeter; - double mu[6], tr, trinv; - char buf[100]; - - /* - if(argc < 2) { - fprintf(stderr, "Too few arguments!\n"); - return 1; - } - */ + mat_workspace *ws; + mat *matrices; + mat tmp; + mat gen[6]; + char buf[100], buf2[100], buf3[100]; + mpq_t s,t; + mpq_t det, tr, trinv; // allocate stuff group = malloc(MAX_ELEMENTS*sizeof(groupelement_t)); - ws = workspace_alloc(3); - matrices = malloc(MAX_ELEMENTS*sizeof(gsl_matrix*)); + ws = mat_workspace_init(3); + matrices = malloc(MAX_ELEMENTS*sizeof(mat)); for(int i = 0; i < MAX_ELEMENTS; i++) - matrices[i] = gsl_matrix_alloc(3, 3); + mat_init(matrices[i], 3); for(int i = 0; i < 6; i++) - gen[i] = gsl_matrix_alloc(3, 3); - tmp = gsl_matrix_alloc(3, 3); - tmp2 = gsl_matrix_alloc(3, 3); - coxeter = gsl_matrix_alloc(3, 3); + mat_init(gen[i], 3); + mat_init(tmp, 3); -#ifdef DRAW_PICTURE - int width = 800; - int height = 800; - printf("P3\n%d %d\n255\n", width, height); + mpq_inits(s, t, det, tr, trinv, NULL); - for(int y = 0; y < height; y++) { - for(int x = 0; x < width; x++) { - double s = exp(((double)x/width-0.5)*4); - double t = exp(((double)y/height-0.5)*4); -#else - double s = atof(argv[1]); - double t = atof(argv[2]); -#endif + mpq_set_ui(s,1,1); + double t_ = atof(argv[1]); + mpq_set_ui(t,(int)(t_*100),100); + mpq_canonicalize(s); + mpq_canonicalize(t); - // the real action - generate_triangle_group(group, MAX_ELEMENTS, 5, 6, 7); - initialize_triangle_generators(ws, gen, 1.0/5.0, 1.0/6.0, 1.0/7.0, s, t); + // the real action + generate_triangle_group(group, MAX_ELEMENTS, 3, 3, 4); + initialize_triangle_generators(ws, gen, s, t); - gsl_matrix_set_identity(matrices[0]); - for(int i = 1; i < MAX_ELEMENTS; i++) { - if(group[i].length % 2 != 0) - continue; + mat_identity(matrices[0]); + for(int i = 1; i < MAX_ELEMENTS; i++) { + if(group[i].length % 2 != 0) + continue; + if(!group[i].inverse) + continue; - int parent = group[i].parent->id; - int grandparent = group[i].parent->parent->id; - int letter; + int parent = group[i].parent->id; + int grandparent = group[i].parent->parent->id; + int letter; - if(group[parent].letter == 0 && group[i].letter == 1) - letter = 3; // p = ab - else if(group[parent].letter == 1 && group[i].letter == 2) - letter = 4; // q = bc - else if(group[parent].letter == 2 && group[i].letter == 0) - letter = 5; // r = ca - else if(group[parent].letter == 1 && group[i].letter == 0) - letter = 0; // p^{-1} = ba - else if(group[parent].letter == 2 && group[i].letter == 1) - letter = 1; // q^{-1} = cb - else if(group[parent].letter == 0 && group[i].letter == 2) - letter = 2; // r^{-1} = ac + if(group[parent].letter == 1 && group[i].letter == 2) + letter = 0; // p = bc + else if(group[parent].letter == 2 && group[i].letter == 0) + letter = 1; // q = ca + else if(group[parent].letter == 0 && group[i].letter == 1) + letter = 2; // r = ab + if(group[parent].letter == 2 && group[i].letter == 1) + letter = 3; // p^{-1} = cb + else if(group[parent].letter == 0 && group[i].letter == 2) + letter = 4; // q^{-1} = ac + else if(group[parent].letter == 1 && group[i].letter == 0) + letter = 5; // r^{-1} = ba + + mat_multiply(ws, matrices[i], matrices[grandparent], gen[letter]); + } + + for(int i = 0; i < MAX_ELEMENTS; i++) { + if(group[i].length % 2 != 0) + continue; + if(!group[i].inverse) + continue; + + mat_trace(tr, matrices[i]); + mat_trace(trinv, matrices[group[i].inverse->id]); + + double lambda1, lambda2, lambda3; +// int realevs = gsl_poly_solve_cubic(-mpq_get_d(tr), mpq_get_d(trinv), -1, &lambda3, &lambda2, &lambda1); +// if(realevs != 3) +// continue; +// if(lambda1 < 0 || lambda2 < 0 || lambda3 < 0) +// continue; + +// gmp_printf("%d %d %s %Qd %Qd\n", i, group[i].length, print_word(&group[i], buf), tr, trinv); + printf("%d %d %s %f %f %f\n", i, group[i].length, print_word(&group[i], buf), log(mpq_get_d(tr)), log(mpq_get_d(trinv))); - multiply(matrices[grandparent], gen[letter], matrices[i]); - } - - double excentricity; - double max_excentricity = 0; - int max_excentricity_index = 0; - - for(int i = 0; i < MAX_ELEMENTS; i++) { - if(group[i].length % 2 != 0) - continue; - - // jordan projection - gsl_matrix_memcpy(tmp, matrices[i]); - gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); - int r = gsl_eigen_nonsymmv(tmp, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); - ERROR(r, "gsl_eigen_nonsymmv failed!\n"); - gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC); - - int real_evs = 0; - - for(int j = 0; j < 3; j++) - if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, j)), 0) == 0) - real_evs++; - - if(real_evs != 3) - continue; - - mu[0] = fabs(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 0)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 1))); - mu[1] = fabs(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 1)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 2))); - - excentricity = log(mu[1])/log(mu[0]); - if(excentricity > max_excentricity + 1e-3) { - max_excentricity = excentricity; - max_excentricity_index = i; - } - -#ifndef DRAW_PICTURE - printf("%d %f %f %f %f 0x0000ff %d %s\n", group[i].length, mu[0], mu[1], tr, trinv, i, print_word(&group[i], buf)); -#endif - } - -#ifdef DRAW_PICTURE - int color = (int)((max_excentricity-1)/(max_excentricity+4)*255); - if(color < 0) - color = 0; - if(color > 255) - color = 255; - printf("%d %d %d\n", color, color, color); - } - fprintf(stderr, "y = %d\n", y); } -#else - fprintf(stderr, "max = %f %s\n", max_excentricity, print_word(&group[max_excentricity_index], buf)); -#endif // free stuff for(int i = 0; i < MAX_ELEMENTS; i++) - gsl_matrix_free(matrices[i]); + mat_clear(matrices[i]); for(int i = 0; i < 6; i++) - gsl_matrix_free(gen[i]); - gsl_matrix_free(tmp); - gsl_matrix_free(tmp2); - gsl_matrix_free(coxeter); - free(matrices); - free(group); - workspace_free(ws); + mat_clear(gen[i]); + mat_clear(tmp); + mpq_clears(s, t, det, tr, trinv, NULL); + mat_workspace_clear(ws); return 0; } diff --git a/singular_values.plt b/singular_values.plt index d6775cc..ce251a7 100644 --- a/singular_values.plt +++ b/singular_values.plt @@ -1,20 +1,38 @@ -if(!exists("logt")) logt = log(1.5) +if(!exists("logt")) logt = log(7) if(!exists("logs")) logs = log(1.0) -file = sprintf("< ./singular_values %f %f", exp(logs), exp(logt)) -title = sprintf("s = %f, t = %f", exp(logs), exp(logt)) +file = sprintf("< ./singular_values %f", exp(logt)) +#title = sprintf("s = %f, t = %f", exp(logs), exp(logt)) +title = sprintf("t = %.3f", floor(exp(logt)*100)/100.0) # print title set zeroaxis +set samples 1000 set size square -set xrange [0:25] -set yrange [0:25] +set xrange [0:20] +set yrange [0:20] +set trange [0:20] set grid set parametric -plot file using (log($2)):(log($3)) w p pt 7 ps 0.5 lc 1 t title # plot file using 2:3 w p pt 7 ps 0.5 lc 1 t title +#tr(a,b) = exp((2*a+b)/3) + exp((b-a)/3) + exp(-(a+2*b)/3) +#trinv(a,b) = exp(-(2*a+b)/3) + exp((a-b)/3) + exp((a+2*b)/3) + +tr(a,b) = exp(a) + exp(b-a) + exp(-b) +trinv(a,b) = exp(-a) + exp(a-b) + exp(b) + +plot file using 4:5 w p pt 7 ps 0.5 lc 1 t title, \ + log(tr(t,t*2)),log(trinv(t,2*t)) w l lw 2 t "", \ + log(tr(t,t/2)),log(trinv(t,t/2)) w l lw 2 t "" + +#plot for[i=-10:10] log(tr(t,t*exp(log(2)*i/10.0))),log(trinv(t,t*exp(log(2)*i/10.0))) w l lw 2 t "" + +#plot for[i=-10:10] t,log(tr(t,t*exp(log(2)*i/10.0)))-t w l lw 2 t "" + +##plot for[i=20:20] t,log(tr(1/t,exp(2*log(2)*i/20.0-log(2)))) w l lw 2 t "" + pause mouse keypress if(MOUSE_KEY == 60) logt=logt-0.02 if(MOUSE_KEY == 62) logt=logt+0.02