From 015b391cc072a152e3214f0bc763a2b3fbc36136 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Sat, 26 Mar 2022 09:42:16 -0500 Subject: [PATCH] draw x^alpha curves; no want to refactor limit curve --- Makefile | 11 ++- draw.c | 216 ++++++++++++++++++++++++++++++------------------- exp_equation.c | 183 +++++++++++++++++++++++++++++++++++++++++ exp_equation.h | 14 ++++ main.c | 4 +- main.h | 8 ++ 6 files changed, 345 insertions(+), 91 deletions(-) create mode 100644 exp_equation.c create mode 100644 exp_equation.h diff --git a/Makefile b/Makefile index 4f13fc4..fb64c4a 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -HEADERS=triangle.h linalg.h queue.h initcairo.h main.h +HEADERS=triangle.h linalg.h queue.h initcairo.h main.h exp_equation.h SPECIAL_OPTIONS=-O0 -g -D_DEBUG #SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline @@ -11,8 +11,8 @@ OPTIONS=$(GENERAL_OPTIONS) $(CAIRO_OPTIONS) $(SPECIAL_OPTIONS) all: limit_set -limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o - gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o -lm -lgsl -lcblas -lcairo -lX11 +limit_set: limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.o + gcc $(OPTIONS) -o limit_set limit_set.o linalg.o triangle.o initcairo.o draw.o main.o exp_equation.o -lm -lgsl -lcblas -lcairo -lX11 linalg.o: linalg.c $(HEADERS) gcc $(OPTIONS) -c linalg.c @@ -32,5 +32,8 @@ draw.o: draw.c $(HEADERS) main.o: main.c $(HEADERS) gcc $(OPTIONS) -c main.c +exp_equation.o: exp_equation.c $(HEADERS) + gcc $(OPTIONS) -c exp_equation.c + clean: - rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o + rm -f limit_set linalg.o triangle.o limit_set.o draw.o main.o exp_equation.o diff --git a/draw.c b/draw.c index 6092222..ddf9658 100644 --- a/draw.c +++ b/draw.c @@ -1,4 +1,5 @@ #include "main.h" +#include "exp_equation.h" #define FMOD(x,y) (fmod(x,y) < 0 ? fmod(x,y) + y : fmod(x,y)) #define ANGLE_DIFF(x,y) (FMOD((x)-(y), 2*M_PI)) @@ -145,6 +146,44 @@ int intersect_line_and_conic(DrawingContext *ctx, vector_t line, gsl_matrix *con releaseTempMatrices(ctx->ws, 1); } +// intersect the line given by the covector "line" with the orbit of "orbit_point" +// by the one--parameter subgroup of SL(3,R) which contains the element "loxodromic" +// in an eigenbasis of "loxodromic", this corresponds +int intersect_line_and_loxodromic_orbit(DrawingContext *ctx, vector_t line, gsl_matrix *frame, double *logeigenvalues, vector_t start, vector_t *out) +{ + vector_t line_in_frame = apply_transpose(frame, line); + vector_t start_in_frame = apply_pseudoinverse(frame, start); + + vector_t a, x; + LOOP(i) a.x[i] = logeigenvalues[i]; + LOOP(i) x.x[i] = line_in_frame.x[i]*start_in_frame.x[i]; + + double t[2]; + vector_t v[2]; + int n1, n2; + + n1 = solve_linear_exp(a, x, t); + for(int i = 0; i < n1; i++) { + LOOP(j) v[i].x[j] = exp(a.x[j]*t[i]) * start_in_frame.x[j]; + out[i] = apply(frame, v[i]); + } + + x.x[1] *= -1; + n2 = solve_linear_exp(a, x, t); + for(int i = 0; i < n2; i++) { + LOOP(j) v[i].x[j] = exp(a.x[j]*t[i]) * start_in_frame.x[j]; + v[i].x[1] *= -1; + out[i+n1] = apply(frame, v[i]); + } + + if(n1+n2 > 2) { + fprintf(stderr, "more than 2 solutions in intersect_line_and_loxodromic_orbit()!\n"); + exit(1); + } + + return n1+n2; +} + // should be three collinear vectors! double halfCR(vector_t x, vector_t y, vector_t z) { @@ -500,6 +539,48 @@ void drawRotationOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, vector_t sta releaseTempVectors(ctx->ws, 2); } +void drawLoxodromicOrbitFrame(DrawingContext *ctx, gsl_matrix *frame, double *logeigenvalues, vector_t start) +{ + vector_t start_in_frame = apply_pseudoinverse(frame, start); + int iterations = 500; + double stepsize = 0.02; + vector_t x, w; + point_t p; + cairo_t *C = ctx->cairo; + double t; + + int previous_inside = 0; + for(int k = 0; k <= iterations; k++) { + // 0 = repelling fixed point, iterations/2 = attracting fixed point + if(k == 0 || k == iterations) { + w.x[0] = 0.0; w.x[1] = 0.0; w.x[2] = 1.0; + } else if(k == iterations/2) { + w.x[0] = 1.0; w.x[1] = 0.0; w.x[2] = 0.0; + } else if(k < iterations/2) { + t = (k-(double)iterations/4.0)*stepsize; + LOOP(i) w.x[i] = start_in_frame.x[i] * exp(logeigenvalues[i]*t); + w.x[1] *= -1; + } else { + t = ((double)iterations*3.0/4.0-k)*stepsize; + LOOP(i) w.x[i] = start_in_frame.x[i] * exp(logeigenvalues[i]*t); + } + + x = apply(frame, w); + p = vectorToPoint(ctx, x); + + if(isInsideBB(ctx, p)) { + if(!previous_inside) + cairo_move_to(C, p.x, p.y); + else + cairo_line_to(C, p.x, p.y); + previous_inside = 1; + } else { + previous_inside = 0; + } + } + + cairo_stroke(C); +} void drawConic(DrawingContext *ctx, gsl_matrix *form) { @@ -892,92 +973,22 @@ void drawBoxes(DrawingContext *ctx) cairo_set_line_width(C, 2.0/ctx->dim->scalefactor); cairo_set_source_rgb(C, 0.6, 0.6, 0.6); -// drawRotationOrbit(ctx, "ab", p[0][0]); -// drawRotationOrbit(ctx, "bc", p[0][0]); -// drawRotationOrbit(ctx, "ca", p[0][0]); + LOOP(i) LOOP(j) gsl_matrix_set(tmp, i, j, p[0][j].x[i]); + double evs[3]; + wordEigenvalues(ctx, "abc", evs); + LOOP(i) evs[i] = log(fabs(evs[i])); + drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][0]); + drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][2]); -// drawRotationOrbit(ctx, "ab", p[0][2]); -// drawRotationOrbit(ctx, "bc", p[0][2]); -// drawRotationOrbit(ctx, "ca", p[0][2]); - -// gsl_matrix_set_zero(order3); -// LOOP(i) gsl_matrix_set(order3, (i+1)%3, i, 1); -// rotation_frame(order3, frame, ctx->ws); -// drawRotationOrbitFrame(ctx, frame, p[0][0]); -// drawRotationOrbitFrame(ctx, frame, p[0][2]); - - vector_t line1; - vector_t line2; - double t; - int positives; - - /* - // conic 1 - line1 = cross(p[0][0],p[0][1]); - line2 = cross(p[1][0],p[1][1]); - degenerate_conic(line1, line2, tmp); - line1 = cross(p[0][0], p[1][0]); - line2 = cross(p[0][0], p[1][0]); - degenerate_conic(line1, line2, tmp2); - t = - conic_value(tmp, p[2][0]) / conic_value(tmp2, p[2][0]); - LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, - gsl_matrix_get(tmp, i, j) + - gsl_matrix_get(tmp2, i, j)*t); - positives = diagonalize_symmetric_form(conic, tmp, ctx->ws); - if(positives == 2) - LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; - drawConic(ctx, conic); - */ - - // conic 2 - line1 = cross(p[0][0],p[0][1]); - line2 = cross(p[2][0],p[2][1]); - degenerate_conic(line1, line2, tmp); - line1 = cross(p[0][0], p[2][0]); - line2 = cross(p[0][0], p[2][0]); - degenerate_conic(line1, line2, tmp2); - t = - conic_value(tmp, p[1][0]) / conic_value(tmp2, p[1][0]); - LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, - gsl_matrix_get(tmp, i, j) + - gsl_matrix_get(tmp2, i, j)*t); - positives = diagonalize_symmetric_form(conic, tmp, ctx->ws); - if(positives == 2) - LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; - drawConic(ctx, conic); - - /* - // conic 3 - line1 = cross(p[0][2],p[0][1]); - line2 = cross(p[1][2],p[1][1]); - degenerate_conic(line1, line2, tmp); - line1 = cross(p[0][2], p[1][2]); - line2 = cross(p[0][2], p[1][2]); - degenerate_conic(line1, line2, tmp2); - t = - conic_value(tmp, p[2][2]) / conic_value(tmp2, p[2][2]); - LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, - gsl_matrix_get(tmp, i, j) + - gsl_matrix_get(tmp2, i, j)*t); - positives = diagonalize_symmetric_form(conic, tmp, ctx->ws); - if(positives == 2) - LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; - drawConic(ctx, conic); - */ - - // conic 4 - line1 = cross(p[0][2],p[0][1]); - line2 = cross(p[2][2],p[2][1]); - degenerate_conic(line1, line2, tmp); - line1 = cross(p[0][2], p[2][2]); - line2 = cross(p[0][2], p[2][2]); - degenerate_conic(line1, line2, tmp2); - t = - conic_value(tmp, p[1][2]) / conic_value(tmp2, p[1][2]); - LOOP(i) LOOP(j) gsl_matrix_set(conic, i, j, - gsl_matrix_get(tmp, i, j) + - gsl_matrix_get(tmp2, i, j)*t); - positives = diagonalize_symmetric_form(conic, tmp, ctx->ws); - if(positives == 2) - LOOP(i) LOOP(j) *gsl_matrix_ptr(conic, i, j) *= -1; - drawConic(ctx, conic); + vector_t x; + for(int i = 0; i < ctx->n_group_elements; i++) { + LOOP(j) x.x[j] = ctx->limit_curve[12*i+3+j]; + x = apply_pseudoinverse(tmp, x); + printf("%f\n", + pow(fabs(x.x[0]), evs[1]-evs[2]) * + pow(fabs(x.x[1]), evs[2]-evs[0]) * + pow(fabs(x.x[2]), evs[0]-evs[1])); + } cairo_restore(C); releaseTempMatrices(ctx->ws, 9); @@ -1101,6 +1112,39 @@ void drawBoxes2(DrawingContext *ctx) perimeter2 += log(value); approx_perimeter2 += log(value); + cairo_set_source_rgb(C, 0, 0.6, 0.1); + + LOOP(i) LOOP(j) gsl_matrix_set(tmp, i, j, p[0][j].x[i]); + double evs[3]; + wordEigenvalues(ctx, "abc", evs); + LOOP(i) evs[i] = log(fabs(evs[i])); + drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][0]); + drawLoxodromicOrbitFrame(ctx, tmp, evs, p[1][2]); + + vector_t intersection[2]; + double hC0, hC1; + + cairo_set_source_rgb(C, 1, 0, 0); + intersect_line_and_loxodromic_orbit(ctx, cross(a,b), tmp, evs, p[1][2], intersection); + hC0 = halfCR(a, b, intersection[0]); + hC1 = halfCR(a, b, intersection[1]); + if(hC0 > hC1) + drawVector(ctx, intersection[0]); + else + drawVector(ctx, intersection[1]); + approx_perimeter1 += log(hC0 > hC1 ? hC0 : hC1); + + cairo_set_source_rgb(C, 0, 0, 1); + intersect_line_and_loxodromic_orbit(ctx, cross(a,b), tmp, evs, p[1][0], intersection); + hC0 = halfCR(b, a, intersection[0]); + hC1 = halfCR(b, a, intersection[1]); + if(hC0 > hC1) + drawVector(ctx, intersection[0]); + else + drawVector(ctx, intersection[1]); + approx_perimeter2 += log(hC0 > hC1 ? hC0 : hC1); + + /* vector_t conic_intersection[2]; double hC0, hC1; @@ -1150,10 +1194,12 @@ void drawBoxes2(DrawingContext *ctx) snprintf(ctx->extra_text, 1000, "perimeter1 = %f (%f), perimeter2 = %f (%f)", perimeter1, approx_perimeter1, perimeter2, approx_perimeter2); + */ printf("%f %f %f %f %f %f\n", t, s, perimeter1, perimeter2, approx_perimeter1, approx_perimeter2); + cairo_restore(C); releaseTempMatrices(ctx->ws, 7); releaseTempVectors(ctx->ws, 2); diff --git a/exp_equation.c b/exp_equation.c new file mode 100644 index 0000000..c256073 --- /dev/null +++ b/exp_equation.c @@ -0,0 +1,183 @@ +#include "main.h" +#include "exp_equation.h" + +#include +#include +#include +#include + +#define EPSILON 1e-9 +#define LOOP(i) for(int i = 0; i < 3; i++) + +struct solve_exp_plus_exp_params +{ + double alpha; + double beta; + double x; +}; + +static double solve_exp_plus_exp_f(double t, void *_params) +{ + struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params; + +// return exp(params->alpha * t) + exp(params->beta * t) - params->x; + return log(exp(params->alpha * t) + exp(params->beta * t)) - log(params->x); +} + +static double solve_exp_plus_exp_df(double t, void *_params) +{ + struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params; + +// return params->alpha * exp(params->alpha * t) + params->beta * exp(params->beta * t); + return params->alpha + (params->beta - params->alpha) / (1 + exp((params->alpha-params->beta)*t)); +} + +static void solve_exp_plus_exp_fdf(double t, void *params, double *f, double *df) +{ + *f = solve_exp_plus_exp_f(t, params); + *df = solve_exp_plus_exp_df(t, params); +} + +// solve the equation exp(alpha t) + exp(beta t) = x for t +int solve_exp_plus_exp(double alpha, double beta, double x, double *t) +{ + if(alpha <= 0 && beta >= 0 || alpha >= 0 && beta <= 0) { + double critical = + pow(-beta/alpha, alpha/(alpha-beta)) + + pow(-beta/alpha, beta/(alpha-beta)); + + if(x < critical) + return 0; + else if (x < critical + EPSILON) { + t[0] = log(-beta/alpha)/(alpha-beta); + return 1; + } + + // Newton this + gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton); + struct solve_exp_plus_exp_params params; + params.alpha = alpha; + params.beta = beta; + params.x = x; + gsl_function_fdf FDF; + FDF.f = &solve_exp_plus_exp_f; + FDF.df = &solve_exp_plus_exp_df; + FDF.fdf = &solve_exp_plus_exp_fdf; + FDF.params = (void *)¶ms; + int status; + double root, lastroot; + + for(int r = 0; r < 2; r++) { + root = r == 0 ? log(x)/beta : log(x)/alpha; + gsl_root_fdfsolver_set(solver, &FDF, root); + for(int i = 0; i < 100; i++) { + gsl_root_fdfsolver_iterate(solver); + lastroot = root; + root = gsl_root_fdfsolver_root(solver); + status = gsl_root_test_delta(root, lastroot, 0, 1e-9); + + if(status == GSL_SUCCESS) + break; + +// printf("iteration %d, root %f\n", i, root); + } + t[r] = root; + } + + gsl_root_fdfsolver_free(solver); + + return 2; + } else { + // Newton with start value 0 + gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton); + struct solve_exp_plus_exp_params params; + params.alpha = alpha; + params.beta = beta; + params.x = x; + gsl_function_fdf FDF; + FDF.f = &solve_exp_plus_exp_f; + FDF.df = &solve_exp_plus_exp_df; + FDF.fdf = &solve_exp_plus_exp_fdf; + FDF.params = (void *)¶ms; + int status; + double root, lastroot; + + root = 0; + gsl_root_fdfsolver_set(solver, &FDF, root); + for(int i = 0; i < 100; i++) { + gsl_root_fdfsolver_iterate(solver); + lastroot = root; + root = gsl_root_fdfsolver_root(solver); + status = gsl_root_test_delta(root, lastroot, 0, 1e-9); + +// printf("iteration %d, root %f\n", i, root); + + if(status == GSL_SUCCESS) + break; + } + t[0] = root; + + gsl_root_fdfsolver_free(solver); + + return 1; + } +} + +// solve the equation x1 exp(a1 t) + x2 exp(a2 t) + x3 exp(a3 t) = 0 +int solve_linear_exp(vector_t a, vector_t x, double *t) +{ + if(x.x[0] > 0 && x.x[1] > 0 && x.x[2] > 0 || + x.x[0] < 0 && x.x[1] < 0 && x.x[2] < 0) + return 0; + + // ensure that y[0] < 0 and y[1], y[2] > 0 + int j; + vector_t y, b; + for(j = 0; j < 3; j++) + if(x.x[(j+1)%3] * x.x[(j+2)%3] > 0) + break; + LOOP(i) y.x[i] = x.x[(i+j)%3]; + if(y.x[0] > 0) + LOOP(i) y.x[i] *= -1; + LOOP(i) b.x[i] = a.x[(i+j)%3]; + + double T = (log(y.x[1]) - log(y.x[2])) / (b.x[2] - b.x[1]); + double rhs = - y.x[0] * + pow(y.x[1], (b.x[0]-b.x[2])/(b.x[2]-b.x[1])) * + pow(y.x[2], (b.x[0]-b.x[1])/(b.x[1]-b.x[2])); + + int n = solve_exp_plus_exp(b.x[1] - b.x[0], b.x[2] - b.x[0], rhs, t); + + for(int i = 0; i < n; i++) + t[i] += T; + + return n; +} + +/* +int main(int argc, char *argv[]) +{ + +// int n = solve_exp_plus_exp(atof(argv[1]), atof(argv[2]), atof(argv[3]), result); + + double result[2]; + vector_t a, x; + a.x[0] = atof(argv[1]); + a.x[1] = atof(argv[2]); + a.x[2] = atof(argv[3]); + x.x[0] = atof(argv[4]); + x.x[1] = atof(argv[5]); + x.x[2] = atof(argv[6]); + int n = solve_linear_exp(a, x, result); + + if(n == 0) + printf("0 results found\n"); + else if(n == 1) + printf("1 result found: %.9f\n", result[0]); + else if(n == 2) + printf("2 results found: %.9f and %.9f\n", result[0], result[1]); + else + printf("%d results found\n", n); + return 0; +} +*/ diff --git a/exp_equation.h b/exp_equation.h new file mode 100644 index 0000000..97e336a --- /dev/null +++ b/exp_equation.h @@ -0,0 +1,14 @@ +#ifndef EXP_EQUATION_H +#define EXP_EQUATION_H + +#include "main.h" + +#include +#include +#include +#include + +int solve_exp_plus_exp(double alpha, double beta, double x, double *t); +int solve_linear_exp(vector_t a, vector_t x, double *t); + +#endif diff --git a/main.c b/main.c index 4a3f0df..67c1e17 100644 --- a/main.c +++ b/main.c @@ -49,8 +49,8 @@ void setupContext(DrawingContext *ctx, int argc, char *argv[]) ctx->marking2.y = -0.11873; ctx->marking3.x = -0.73679; ctx->marking3.y = -0.21873; - ctx->distance_parameter1 = 0.95; - ctx->distance_parameter2 = 0.95; + ctx->distance_parameter1 = 0.45; + ctx->distance_parameter2 = 0.2; ctx->show_coxeter_orbit = 0; ctx->extra_text = malloc(1000*sizeof(char)); memset(ctx->extra_text, 0, 1000*sizeof(char)); diff --git a/main.h b/main.h index 5889879..b0fbd45 100644 --- a/main.h +++ b/main.h @@ -22,6 +22,14 @@ typedef struct { double y; } point_t; +typdef struct { + double x; + double y; + double angle; + vector_t fp[3]; + groupelement_t *g; +} limit_curve_t; + typedef struct { // infos about the screen to draw on cairo_t *cairo;