delete all files not useful to complex calculation
This commit is contained in:
Florian Stecker committed
1 parent
1a85a5ef14
commit
3309c37955
23 files changed
-1531
No files matched your search
-41
@@ -1,41 +0,0 @@
|
||||
set size ratio 1.3
|
||||
set grid
|
||||
set xrange [1/2.5:2.5]
|
||||
set yrange [0.2:5]
|
||||
set samples 2000
|
||||
set parametric
|
||||
|
||||
set log x
|
||||
set log y
|
||||
|
||||
set terminal pngcairo size 1000, 1200
|
||||
set output "barbot-symmetric.png"
|
||||
|
||||
# solves the equation rho(s)q^2 - (s+1)q + 1 = 0
|
||||
# (s^2 + cs + 1)q^2 - (s+1)q + 1 = q^2 s^2 + (cq - 1)qs + (q^2 - q + 1)
|
||||
c = 2*cos(4*pi/5)
|
||||
rho(s) = s**2 + c*s + 1
|
||||
# qp(s) = (s+1+sqrt((s+1)**2 - 4*rho(s)))/(2*rho(s))
|
||||
# qm(s) = (s+1-sqrt((s+1)**2 - 4*rho(s)))/(2*rho(s))
|
||||
# d(q) = q**2*(c*q-1)**2 - 4*(q**2-q+1)*q**2
|
||||
# sp(q) = (-(c*q-1)*q + sqrt(d(q)))/(2*q**2)
|
||||
# sm(q) = (-(c*q-1)*q - sqrt(d(q)))/(2*q**2)
|
||||
qp(s) = (s+1+sqrt((s+1)**2 - 4*rho(s)))/(2*sqrt(rho(s)))
|
||||
qm(s) = (s+1-sqrt((s+1)**2 - 4*rho(s)))/(2*sqrt(rho(s)))
|
||||
d(q) = q**2*(c*q-1)**2 - 4*(q**2-q+1)*q**2
|
||||
sp(q) = (-(c*q-1)*q + sqrt(d(q)))/(2*q**2)
|
||||
sm(q) = (-(c*q-1)*q - sqrt(d(q)))/(2*q**2)
|
||||
|
||||
plot "output/barbot_map_5000_frequencies" \
|
||||
using ($1/100.0):($2/100.0*sqrt(rho($1/100.0))):(($3==0)?1/0:$4) w p pt 5 ps 0.5 lc palette t '', \
|
||||
"output/barbot_map_500000" \
|
||||
using ($1/10.0) :($2/10.0*sqrt(rho($1/10.0))) :($3*0.7+0.3) w p pt 7 ps variable lc 3 t '', \
|
||||
sm(t), d(t) > 0 ? t*sqrt(rho(sm(t))) : 1/0 w l lw 2 lc 7 t "", \
|
||||
sp(t), d(t) > 0 ? t*sqrt(rho(sp(t))) : 1/0 w l lw 2 lc 7 t "", \
|
||||
t, qm(t) w l lw 2 lc 7 t "", \
|
||||
t, qp(t) w l lw 2 lc 7 t ""
|
||||
|
||||
# plot "output/barbot_map_50000_zoom_b" using ($1/1000.0):($2/5000):(($3==0)?1/0:$4) w p pt 5 ps 0.5 lc palette t ''
|
||||
|
||||
# pause mouse keypress
|
||||
# reread
|
||||
Binary file not shown.
Binary file not shown.
@@ -1,31 +0,0 @@
|
||||
#!/bin/bash
|
||||
|
||||
trap 'exit 130' INT
|
||||
|
||||
wordlength=30
|
||||
sdenom=1
|
||||
sstart=1
|
||||
send=1
|
||||
qdenom=100
|
||||
qstart=1
|
||||
qend=200 # 1/sqrt(2) = 0.7071...
|
||||
|
||||
words="$(./billiard_words $wordlength | awk '{print $1}')"
|
||||
#words="cbabcabacabcacbcab cabacabcacbcabcbab cabcacbcabcbabcaba"
|
||||
#words="abcabc abcb cbabcacbcacbab"
|
||||
#words="abcabc abcbcabcbc"
|
||||
#words="abcabc bcbab bcbabcac"
|
||||
|
||||
for s in $(seq $sstart $send); do
|
||||
for q in $(seq $qstart $qend); do
|
||||
i=0
|
||||
echo -n "$s/$sdenom $q/$qdenom "
|
||||
# MAXIMUM=only ./special_element $s/$sdenom $q/$qdenom $words
|
||||
# MAXIMUM=no ./special_element $s/$sdenom $q/$qdenom abcb
|
||||
MAXIMUM=no ./special_element $s/$sdenom $q/$qdenom $words | while read line; do
|
||||
echo -n "$line "
|
||||
((i=i+1))
|
||||
done
|
||||
echo
|
||||
done
|
||||
done
|
||||
@@ -1,62 +0,0 @@
|
||||
import Data.List
|
||||
import Data.Ord
|
||||
import Text.Printf
|
||||
import System.Environment
|
||||
|
||||
main = do
|
||||
argv <- getArgs
|
||||
listWordsUpToLength $ read $ argv !! 0
|
||||
|
||||
listWordsUpToLength :: Int -> IO ()
|
||||
listWordsUpToLength n = do
|
||||
putStr $ unlines [printf "%s %d/%d %f"
|
||||
w
|
||||
(p `div` gcd p q)
|
||||
(q `div` gcd p q)
|
||||
(atan (sqrt 3 / (2*q_/p_ + 1))) |
|
||||
((p,q),w) <- wordlist (n `div` 2, n `div` 2),
|
||||
let p_ = fromIntegral p :: Double,
|
||||
let q_ = fromIntegral q :: Double,
|
||||
length w <= n,
|
||||
let x = 2*q + p,
|
||||
let y = 2*p + q]
|
||||
|
||||
-- (sqrt 3 / 2 * fromIntegral p / (fromIntegral q + fromIntegral p / 2) :: Double) |
|
||||
-- (slopeWord "bca" (orthogonalSlope (p,q))) |
|
||||
|
||||
wordlist :: (Int,Int) -> [((Int,Int),String)]
|
||||
wordlist (pmax,qmax) = nub $
|
||||
sortBy (comparing sl)
|
||||
[((p `div` gcd p q, q `div` gcd p q), slopeWord "bca" (p,q)) |
|
||||
p <- [0..pmax],
|
||||
q <- [0..qmax],
|
||||
q /= 0] -- use p /= 0 || q /= 0 for more symmetric output
|
||||
where
|
||||
sl ((p,q),_) = fromIntegral p / fromIntegral q
|
||||
|
||||
-- letters: reflection along e_1, reflection along e_2, other one; p,q >= 0
|
||||
-- the "slope" (p,q) means the Euclidean vector q*e_1 + p*e_2, where e_1,e_2 are at a 60 degree angle
|
||||
-- in Euclidean coordinates this is (q + p/2, sqrt(3)/2 * p)
|
||||
slopeWord :: [Char] -> (Int,Int) -> String
|
||||
slopeWord [x,y,z] (p,q)
|
||||
| p > q = slopeWord [y,x,z] (q,p)
|
||||
| otherwise = concat $ map word $ zipWith step list (tail list)
|
||||
where
|
||||
p_ = p `div` gcd p q
|
||||
q_ = q `div` gcd p q
|
||||
xmax = if (p_-q_) `mod` 3 == 0 then q_ else 3*q_ :: Int
|
||||
list = [(x,(x*p) `div` q) | x <- [0..xmax]]
|
||||
step (x1,y1) (x2,y2) = ((x1-y1) `mod` 3, y2-y1)
|
||||
word (0,0) = [z,x]
|
||||
word (1,0) = [y,z]
|
||||
word (2,0) = [x,y]
|
||||
word (0,1) = [z,y,x,y]
|
||||
word (1,1) = [y,x,z,x]
|
||||
word (2,1) = [x,z,y,z]
|
||||
|
||||
-- assuming p, q >= 0
|
||||
orthogonalSlope :: (Int, Int) -> (Int, Int)
|
||||
orthogonalSlope (p,q)
|
||||
| p > q = (p-q, p+2*q)
|
||||
| p < q = (q+2*p, q-p)
|
||||
| otherwise = (1,0)
|
||||
@@ -1,57 +0,0 @@
|
||||
#if(!exists("logt")) logt = log(1.80)
|
||||
if(!exists("n")) n = 263
|
||||
if(!exists("logt")) logt = log(1)
|
||||
if(!exists("logs")) logs = log(1)
|
||||
|
||||
#logt = 0.01*n
|
||||
logt = log(1000000000)
|
||||
|
||||
file = sprintf("< ./singular_values 713698 %f %f", exp(logs), exp(logt))
|
||||
#file = sprintf("< ./singular_values 1621 %f %f", exp(logs), exp(logt))
|
||||
#outfile = sprintf("cdf/cdf_hires_%05d.png", n)
|
||||
outfile = sprintf("cdf/cdf_hires_limit.png")
|
||||
|
||||
set log x
|
||||
set zeroaxis
|
||||
set samples 1000
|
||||
set size square
|
||||
set xrange [0.5:2]
|
||||
set yrange [0:500000]
|
||||
#set yrange [0:1000]
|
||||
set trange [0:30]
|
||||
set grid
|
||||
set parametric
|
||||
|
||||
set terminal pngcairo enhanced size 1024, 1024
|
||||
set output outfile
|
||||
print sprintf("n = %d, t = %.2f", n, exp(logt))
|
||||
|
||||
# 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 6:7 w p pt 7 ps 0.5 lc 1 t columnheader,
|
||||
# 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 file using 8:3 w steps lw 2 lc 1 t sprintf("t = %.2f", exp(logt))
|
||||
|
||||
#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 ""
|
||||
|
||||
#n=n+1
|
||||
#if(n < 1000) reread
|
||||
|
||||
# pause mouse keypress
|
||||
# if(MOUSE_KEY == 60) logt=logt-0.02
|
||||
# if(MOUSE_KEY == 62) logt=logt+0.02
|
||||
# if(MOUSE_KEY == 44) logs=logs-0.02
|
||||
# if(MOUSE_KEY == 46) logs=logs+0.02
|
||||
# if(MOUSE_KEY != 113) reread
|
||||
@@ -1 +0,0 @@
|
||||
ffmpeg -f image2 -framerate 20 -i test%03d.png -s 1024x1024 -c:v libvpx-vp9 -lossless 1 test.webm
|
||||
@@ -1,353 +0,0 @@
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <gsl/gsl_math.h>
|
||||
#include <gsl/gsl_eigen.h>
|
||||
#include <gsl/gsl_blas.h>
|
||||
#include <gsl/gsl_linalg.h>
|
||||
#include <gsl/gsl_complex_math.h>
|
||||
#include <memory.h>
|
||||
#include <math.h>
|
||||
|
||||
#include "linalg.h"
|
||||
|
||||
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
|
||||
#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
|
||||
|
||||
/*********************************************** temporary storage ********************************************************/
|
||||
|
||||
workspace_t *workspace_alloc(int n)
|
||||
{
|
||||
workspace_t *result = (workspace_t*)malloc(sizeof(workspace_t));
|
||||
result->n = n;
|
||||
result->work_nonsymmv = gsl_eigen_nonsymmv_alloc(n);
|
||||
result->work_symmv = gsl_eigen_symmv_alloc(n);
|
||||
result->work_sv = gsl_vector_alloc(n);
|
||||
result->eval_complex = gsl_vector_complex_alloc(n);
|
||||
result->evec_complex = gsl_matrix_complex_alloc(n, n);
|
||||
result->eval_real = gsl_vector_alloc(n);
|
||||
result->evec_real = gsl_matrix_alloc(n, n);
|
||||
result->permutation = gsl_permutation_alloc(n);
|
||||
|
||||
result->tmp_mat = malloc(MAX_TEMP_MATRICES*sizeof(gsl_matrix*));
|
||||
for(int i = 0; i < MAX_TEMP_MATRICES; i++)
|
||||
result->tmp_mat[i] = gsl_matrix_alloc(3, 3);
|
||||
result->tmp_mat_used = 0;
|
||||
result->tmp_vec = malloc(MAX_TEMP_MATRICES*sizeof(gsl_vector*));
|
||||
for(int i = 0; i < MAX_TEMP_MATRICES; i++)
|
||||
result->tmp_vec[i] = gsl_vector_alloc(3);
|
||||
result->tmp_vec_used = 0;
|
||||
return result;
|
||||
}
|
||||
|
||||
void workspace_free(workspace_t *workspace)
|
||||
{
|
||||
gsl_eigen_nonsymmv_free(workspace->work_nonsymmv);
|
||||
gsl_eigen_symmv_free(workspace->work_symmv);
|
||||
gsl_vector_free(workspace->work_sv);
|
||||
gsl_vector_complex_free(workspace->eval_complex);
|
||||
gsl_matrix_complex_free(workspace->evec_complex);
|
||||
gsl_vector_free(workspace->eval_real);
|
||||
gsl_matrix_free(workspace->evec_real);
|
||||
gsl_permutation_free(workspace->permutation);
|
||||
|
||||
for(int i = 0; i < MAX_TEMP_MATRICES; i++)
|
||||
gsl_matrix_free(workspace->tmp_mat[i]);
|
||||
free(workspace->tmp_mat);
|
||||
for(int i = 0; i < MAX_TEMP_VECTORS; i++)
|
||||
gsl_vector_free(workspace->tmp_vec[i]);
|
||||
free(workspace->tmp_vec);
|
||||
}
|
||||
|
||||
/************************************************** basic operations ********************************************************/
|
||||
|
||||
void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws)
|
||||
{
|
||||
int s;
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
|
||||
gsl_matrix_memcpy(tmp, in);
|
||||
gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
|
||||
gsl_linalg_LU_invert(tmp, ws->permutation, out);
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
}
|
||||
|
||||
void solve(gsl_matrix *A, gsl_vector *b, gsl_vector *result, workspace_t *ws)
|
||||
{
|
||||
int s;
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
|
||||
gsl_matrix_memcpy(tmp, A);
|
||||
gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
|
||||
gsl_linalg_LU_solve(tmp, ws->permutation, b, result);
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
}
|
||||
|
||||
void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
|
||||
invert(conjugator, out, ws); // use out to temporarily store inverse conjugator
|
||||
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, in, out, 0.0, tmp); // in * conjugator^{-1}
|
||||
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, conjugator, tmp, 0.0, out);
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
}
|
||||
|
||||
void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out)
|
||||
{
|
||||
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, out);
|
||||
}
|
||||
|
||||
void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, tmp);
|
||||
gsl_matrix_memcpy(a, tmp);
|
||||
releaseTempMatrices(ws, 1);
|
||||
}
|
||||
|
||||
void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, tmp);
|
||||
gsl_matrix_memcpy(b, tmp);
|
||||
releaseTempMatrices(ws, 1);
|
||||
}
|
||||
|
||||
void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...)
|
||||
{
|
||||
va_list args;
|
||||
va_start(args, n);
|
||||
|
||||
gsl_matrix_set_identity(out);
|
||||
|
||||
for(int i = 0; i < n; i++) {
|
||||
gsl_matrix *cur = va_arg(args, gsl_matrix *);
|
||||
multiply_right(out, cur, ws);
|
||||
}
|
||||
|
||||
va_end(args);
|
||||
}
|
||||
|
||||
void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
|
||||
gsl_matrix_memcpy(tmp, g);
|
||||
gsl_linalg_SV_decomp(tmp, ws->evec_real, ws->eval_real, ws->work_sv);
|
||||
|
||||
for(int i = 0; i < ws->n - 1; i++)
|
||||
mu[i] = log(gsl_vector_get(ws->eval_real, i) / gsl_vector_get(ws->eval_real, i+1));
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
}
|
||||
|
||||
void initialize(gsl_matrix *g, double *data, int x, int y)
|
||||
{
|
||||
gsl_matrix_view view = gsl_matrix_view_array(data, x, y);
|
||||
gsl_matrix_memcpy(g, &view.matrix);
|
||||
}
|
||||
|
||||
void rotation_matrix(gsl_matrix *g, double *vector)
|
||||
{
|
||||
double normalized[3];
|
||||
double norm = sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]);
|
||||
for(int i = 0; i < 3; i++)
|
||||
normalized[i] = vector[i] / norm;
|
||||
gsl_matrix_set_identity(g);
|
||||
gsl_matrix_set(g, 0, 0, cos(norm));
|
||||
gsl_matrix_set(g, 0, 1, -sin(norm) * normalized[2]);
|
||||
gsl_matrix_set(g, 0, 2, +sin(norm) * normalized[1]);
|
||||
gsl_matrix_set(g, 1, 0, +sin(norm) * normalized[2]);
|
||||
gsl_matrix_set(g, 1, 1, cos(norm));
|
||||
gsl_matrix_set(g, 1, 2, -sin(norm) * normalized[0]);
|
||||
gsl_matrix_set(g, 2, 0, -sin(norm) * normalized[1]);
|
||||
gsl_matrix_set(g, 2, 1, +sin(norm) * normalized[0]);
|
||||
gsl_matrix_set(g, 2, 2, cos(norm));
|
||||
for(int i = 0; i < 3; i++)
|
||||
for(int j = 0; j < 3; j++)
|
||||
g->data[i * g->tda + j] += (1 - cos(norm)) * normalized[i] * normalized[j];
|
||||
}
|
||||
|
||||
double trace(gsl_matrix *g)
|
||||
{
|
||||
return gsl_matrix_get(g, 0, 0) + gsl_matrix_get(g, 1, 1) + gsl_matrix_get(g, 2, 2);
|
||||
}
|
||||
|
||||
double determinant(gsl_matrix *g, workspace_t *ws)
|
||||
{
|
||||
int s;
|
||||
double result;
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
|
||||
gsl_matrix_memcpy(tmp, g);
|
||||
gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
|
||||
result = gsl_linalg_LU_det(tmp, s);
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *g_ = getTempMatrix(ws);
|
||||
int success = 0;
|
||||
|
||||
gsl_matrix_memcpy(g_, g);
|
||||
gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
|
||||
int r = gsl_eigen_nonsymmv(g_, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
|
||||
ERROR(r, "gsl_eigen_nonsymmv failed!\n");
|
||||
|
||||
gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
|
||||
|
||||
int real = 1;
|
||||
for(int i = 0; i < ws->n; i++)
|
||||
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
|
||||
real = 0;
|
||||
|
||||
if(!real)
|
||||
goto eigenvectors_out;
|
||||
|
||||
for(int i = 0; i < ws->n; i++)
|
||||
for(int j = 0; j < ws->n; j++)
|
||||
gsl_matrix_set(evec_real, i, j, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j)));
|
||||
|
||||
success = 1;
|
||||
|
||||
eigenvectors_out:
|
||||
releaseTempMatrices(ws, 1);
|
||||
return success;
|
||||
}
|
||||
|
||||
// only fills in the real eigenvectors and returns their count
|
||||
int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *g_ = getTempMatrix(ws);
|
||||
|
||||
gsl_matrix_memcpy(g_, g);
|
||||
gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv);
|
||||
int r = gsl_eigen_nonsymmv(g_, 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 = 0;
|
||||
|
||||
for(int i = 0; i < ws->n; i++) {
|
||||
if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) == 0) {// real
|
||||
if(evec_real) {
|
||||
for(int j = 0; j < ws->n; j++)
|
||||
gsl_matrix_set(evec_real, j, real, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, j, i)));
|
||||
}
|
||||
real++;
|
||||
}
|
||||
}
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
return real;
|
||||
}
|
||||
|
||||
void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *g_ = getTempMatrix(ws);
|
||||
|
||||
gsl_matrix_memcpy(g_, g);
|
||||
int r = gsl_eigen_symmv (g_, eval, evec, ws->work_symmv);
|
||||
ERROR(r, "gsl_eigen_symmv failed!\n");
|
||||
|
||||
gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
}
|
||||
|
||||
// returns number of positive directions and matrix transforming TO diagonal basis
|
||||
int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *A_ = getTempMatrix(ws);
|
||||
|
||||
gsl_matrix_memcpy(A_, A);
|
||||
int r = gsl_eigen_symmv (A_, ws->eval_real, cob, ws->work_symmv);
|
||||
ERROR(r, "gsl_eigen_symmv failed!\n");
|
||||
|
||||
gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC);
|
||||
|
||||
gsl_matrix_transpose(cob);
|
||||
|
||||
int positive = 0;
|
||||
|
||||
for(int i = 0; i < ws->n; i++) {
|
||||
if(gsl_vector_get(ws->eval_real, i) > 0)
|
||||
positive++;
|
||||
|
||||
for(int j = 0; j < ws->n; j++)
|
||||
*gsl_matrix_ptr(cob, i, j) *= sqrt(fabs(gsl_vector_get(ws->eval_real, i)));
|
||||
}
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
return positive;
|
||||
}
|
||||
|
||||
// computes a matrix in SL(3, R) which projectively transforms (e1, e2, e3, e1+e2+e3) to the 4 given vectors
|
||||
void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
gsl_vector *coeff = getTempVector(ws);
|
||||
int s;
|
||||
double det, scale;
|
||||
|
||||
for(int i = 0; i < 3; i++)
|
||||
for(int j = 0; j < 3; j++)
|
||||
gsl_matrix_set(tmp, i, j, gsl_vector_get(vertices[j], i));
|
||||
|
||||
gsl_linalg_LU_decomp(tmp, ws->permutation, &s);
|
||||
gsl_linalg_LU_solve(tmp, ws->permutation, vertices[3], coeff);
|
||||
det = gsl_linalg_LU_det(tmp, s);
|
||||
|
||||
for(int i = 0; i < 3; i++)
|
||||
det *= gsl_vector_get(coeff, i);
|
||||
scale = 1/cbrt(det);
|
||||
|
||||
for(int i = 0; i < 3; i++)
|
||||
for(int j = 0; j < 3; j++)
|
||||
gsl_matrix_set(result, i, j, scale*gsl_vector_get(vertices[j], i)*gsl_vector_get(coeff, j));
|
||||
|
||||
releaseTempMatrices(ws, 1);
|
||||
releaseTempVectors(ws, 1);
|
||||
}
|
||||
|
||||
void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws)
|
||||
{
|
||||
gsl_matrix *tmp = getTempMatrix(ws);
|
||||
gsl_matrix *rot_basis = getTempMatrix(ws);
|
||||
|
||||
gsl_matrix_memcpy(tmp, rotation);
|
||||
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");
|
||||
|
||||
double arg, minarg = 5; // greater than pi
|
||||
int minidx;
|
||||
for(int i = 0; i < 3; i++) {
|
||||
arg = gsl_complex_arg(gsl_vector_complex_get(ws->eval_complex, i));
|
||||
if(abs(arg) < minarg)
|
||||
{
|
||||
minidx = i;
|
||||
minarg = abs(arg);
|
||||
}
|
||||
}
|
||||
ERROR(FCMP(minarg, 0.0) != 0, "rotation_frame() failed! No eigenvalue was 1.\n");
|
||||
|
||||
for(int i = 0; i < 3; i++) {
|
||||
gsl_complex x = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+1)%3);
|
||||
gsl_complex y = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+2)%3);
|
||||
gsl_complex z = gsl_matrix_complex_get(ws->evec_complex, i, minidx);
|
||||
gsl_matrix_set(result, i, 0, GSL_REAL(x)+GSL_REAL(y));
|
||||
gsl_matrix_set(result, i, 1, GSL_IMAG(x)-GSL_IMAG(y));
|
||||
gsl_matrix_set(result, i, 2, GSL_REAL(z));
|
||||
}
|
||||
|
||||
releaseTempMatrices(ws, 2);
|
||||
}
|
||||
@@ -1,98 +0,0 @@
|
||||
#ifndef LINALG_H
|
||||
#define LINALG_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <stdarg.h>
|
||||
#include <gsl/gsl_math.h>
|
||||
#include <gsl/gsl_eigen.h>
|
||||
#include <gsl/gsl_blas.h>
|
||||
#include <gsl/gsl_linalg.h>
|
||||
#include <memory.h>
|
||||
|
||||
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
|
||||
#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
|
||||
|
||||
#define MAX_TEMP_MATRICES 10
|
||||
#define MAX_TEMP_VECTORS 10
|
||||
|
||||
typedef struct _workspace {
|
||||
int n;
|
||||
gsl_eigen_nonsymmv_workspace *work_nonsymmv;
|
||||
gsl_eigen_symmv_workspace *work_symmv;
|
||||
gsl_vector *work_sv;
|
||||
gsl_vector_complex *eval_complex;
|
||||
gsl_matrix_complex *evec_complex;
|
||||
gsl_vector *eval_real;
|
||||
gsl_matrix *evec_real;
|
||||
gsl_permutation *permutation;
|
||||
|
||||
gsl_matrix **tmp_mat;
|
||||
int tmp_mat_used;
|
||||
gsl_vector **tmp_vec;
|
||||
int tmp_vec_used;
|
||||
} workspace_t;
|
||||
|
||||
workspace_t *workspace_alloc(int n);
|
||||
void workspace_free(workspace_t *workspace);
|
||||
void solve(gsl_matrix *A, gsl_vector *b, gsl_vector *result, workspace_t *ws);
|
||||
void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws);
|
||||
void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws);
|
||||
void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out);
|
||||
void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws);
|
||||
void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws);
|
||||
void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...);
|
||||
void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
|
||||
void initialize(gsl_matrix *g, double *data, int x, int y);
|
||||
void rotation_matrix(gsl_matrix *g, double *vector);
|
||||
int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
|
||||
double trace(gsl_matrix *g);
|
||||
double determinant(gsl_matrix *g, workspace_t *ws);
|
||||
int eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
|
||||
int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws);
|
||||
void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws);
|
||||
int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws);
|
||||
void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws);
|
||||
void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws);
|
||||
|
||||
// matrix allocation stuff
|
||||
|
||||
static gsl_matrix **getTempMatrices(workspace_t *ws, int n)
|
||||
{
|
||||
ERROR(ws->tmp_mat_used + n > MAX_TEMP_MATRICES, "Ran out of temporary matrices. Consider increasing MAX_TEMP_MATRICES\n");
|
||||
int index = ws->tmp_mat_used;
|
||||
ws->tmp_mat_used += n;
|
||||
return ws->tmp_mat + index;
|
||||
}
|
||||
|
||||
static gsl_matrix *getTempMatrix(workspace_t *ws)
|
||||
{
|
||||
return *getTempMatrices(ws, 1);
|
||||
}
|
||||
|
||||
static void releaseTempMatrices(workspace_t *ws, int n)
|
||||
{
|
||||
ERROR(ws->tmp_mat_used - n < 0, "Released more matrices then in use\n");
|
||||
ws->tmp_mat_used -= n;
|
||||
}
|
||||
|
||||
static gsl_vector **getTempVectors(workspace_t *ws, int n)
|
||||
{
|
||||
ERROR(ws->tmp_vec_used + n > MAX_TEMP_VECTORS, "Ran out of temporary vectors. Consider increasing MAX_TEMP_VECTORS\n");
|
||||
int index = ws->tmp_vec_used;
|
||||
ws->tmp_vec_used += n;
|
||||
return ws->tmp_vec + index;
|
||||
}
|
||||
|
||||
static gsl_vector *getTempVector(workspace_t *ws)
|
||||
{
|
||||
return *getTempVectors(ws, 1);
|
||||
}
|
||||
|
||||
static void releaseTempVectors(workspace_t *ws, int n)
|
||||
{
|
||||
ERROR(ws->tmp_vec_used - n < 0, "Released more vectors then in use\n");
|
||||
ws->tmp_vec_used -= n;
|
||||
}
|
||||
|
||||
#endif
|
||||
@@ -1,18 +0,0 @@
|
||||
set log x
|
||||
set y2tics
|
||||
set xrange [exp(-1):1]
|
||||
set yrange [1.5:2]
|
||||
set y2range [1.98:2.1]
|
||||
set grid
|
||||
|
||||
# set terminal pngcairo enhanced size 1500,1000
|
||||
# set output "output/max_slope.png"
|
||||
|
||||
plot "output/max_slope_1621.dat" using 1:3 w lp pt 7 ps 0.6 lw 2 t "1621 elements", \
|
||||
"output/max_slope_24428.dat" using 1:3 w lp pt 7 ps 0.6 lw 2 t "24428 elements", \
|
||||
"output/max_slope_94252.dat" using 1:3 w lp lw 2 pt 7 ps 0.6 t "94252 elements", \
|
||||
"output/max_slope_713698.dat" using 1:3 w lp lw 2 pt 7 ps 0.6 t "713698 elements", \
|
||||
"output/max_slope_1621.dat" using 1:2 w p pt 7 ax x1y2 t "parameter"
|
||||
|
||||
pause mouse keypress
|
||||
if(MOUSE_KEY != 113) reread
|
||||
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
@@ -1,32 +0,0 @@
|
||||
if(!exists("logt")) logt = log(1)
|
||||
if(!exists("logs")) logs = log(1)
|
||||
|
||||
#file = sprintf("< ./singular_values 713698 %f %f", exp(logs), exp(logt))
|
||||
file = sprintf("< ./singular_values 1621 %f %f", exp(logs), exp(logt))
|
||||
|
||||
set zeroaxis
|
||||
set samples 1000
|
||||
set size square
|
||||
set xrange [0:30]
|
||||
set yrange [0:30]
|
||||
set trange [0:5]
|
||||
set grid
|
||||
set parametric
|
||||
|
||||
plot file using 8:9 w p pt 7 ps 1 lc 1 t sprintf("t = %.2f", exp(logt))
|
||||
|
||||
#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 ""
|
||||
|
||||
#n=n+1
|
||||
#if(n < 1000) reread
|
||||
|
||||
pause mouse keypress
|
||||
if(MOUSE_KEY == 60) logt=logt-0.02
|
||||
if(MOUSE_KEY == 62) logt=logt+0.02
|
||||
if(MOUSE_KEY == 44) logs=logs-0.02
|
||||
if(MOUSE_KEY == 46) logs=logs+0.02
|
||||
if(MOUSE_KEY != 113) reread
|
||||
@@ -1,23 +0,0 @@
|
||||
if(!exists("i")) i = 0
|
||||
|
||||
file = sprintf("< ./singular_values %f 1.78", exp((i-50)*0.02))
|
||||
|
||||
set samples 1000
|
||||
set size square
|
||||
set xrange [0:30]
|
||||
set yrange [0:30]
|
||||
set trange [0:30]
|
||||
set grid
|
||||
set parametric
|
||||
|
||||
set terminal pngcairo enhanced size 1024,1024
|
||||
img = sprintf("output/animation/test%03d.png", i);
|
||||
print sprintf("write %s", img)
|
||||
set output img
|
||||
|
||||
plot file using 6:7 w p pt 7 ps 0.5 lc 1 t columnheader, \
|
||||
t,2*t w l lw 2 t "", \
|
||||
t,t/2 w l lw 2 t ""
|
||||
|
||||
i=i+1
|
||||
if(i <= 100) reread
|
||||
@@ -1,642 +0,0 @@
|
||||
#include "coxeter.h"
|
||||
//#include "linalg.h"
|
||||
#include "mat.h"
|
||||
|
||||
//#include <gsl/gsl_poly.h>
|
||||
#include <mps/mps.h>
|
||||
#include <mpi.h>
|
||||
#include <sys/stat.h>
|
||||
#include <sys/mman.h>
|
||||
#include <fcntl.h>
|
||||
#include <errno.h>
|
||||
#include <string.h>
|
||||
#include <unistd.h>
|
||||
|
||||
#define MIN(x,y) ((x)<(y)?(x):(y))
|
||||
#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
|
||||
#define DEBUG(msg, ...) do { print_time(); fprintf(stderr, msg, ##__VA_ARGS__); } while (0);
|
||||
//#define DEBUG(msg, ...)
|
||||
|
||||
#define TDIV 10
|
||||
#define TFROM 1
|
||||
#define TTO 9
|
||||
#define MDIV 10
|
||||
#define MFROM 1
|
||||
#define MTO 9
|
||||
#define JOBNR(t,m) (((t)-TFROM) + ((m)-MFROM)*(TTO-TFROM+1))
|
||||
#define NJOBS ((TTO-TFROM+1)*(MTO-MFROM+1))
|
||||
#define FLUSH_INTERVAL 100
|
||||
|
||||
enum message_tag {
|
||||
JOB_ORDER,
|
||||
JOB_RESULT,
|
||||
JOB_SHUTDOWN,
|
||||
};
|
||||
|
||||
struct job {
|
||||
int tparam, mparam;
|
||||
int done;
|
||||
double max_slope;
|
||||
double time;
|
||||
};
|
||||
|
||||
struct result {
|
||||
mpq_t tr;
|
||||
mpq_t trinv;
|
||||
};
|
||||
|
||||
struct global_data {
|
||||
int n;
|
||||
group_t *group;
|
||||
mat* matrices;
|
||||
struct result *invariants;
|
||||
struct result **distinct_invariants;
|
||||
mps_context *solver;
|
||||
};
|
||||
|
||||
|
||||
struct timespec starttime;
|
||||
char processor_name[MPI_MAX_PROCESSOR_NAME];
|
||||
int world_rank;
|
||||
int world_size;
|
||||
MPI_Datatype job_datatype;
|
||||
|
||||
void print_time()
|
||||
{
|
||||
double diff;
|
||||
struct timespec current;
|
||||
|
||||
clock_gettime(CLOCK_REALTIME, ¤t);
|
||||
|
||||
diff = (current.tv_sec - starttime.tv_sec) + (current.tv_nsec - starttime.tv_nsec)*1e-9;
|
||||
|
||||
fprintf(stderr, "[%04d %.3f] ", world_rank, diff);
|
||||
}
|
||||
|
||||
static struct global_data allocate_global_data(int n)
|
||||
{
|
||||
struct global_data result;
|
||||
result.n = n;
|
||||
result.matrices = malloc(n*sizeof(mat));
|
||||
for(int i = 0; i < n; i++)
|
||||
mat_init(result.matrices[i], 3);
|
||||
result.invariants = malloc(n*sizeof(struct result));
|
||||
result.distinct_invariants = malloc(n*sizeof(struct result*));
|
||||
for(int i = 0; i < n; i++) {
|
||||
mpq_init(result.invariants[i].tr);
|
||||
mpq_init(result.invariants[i].trinv);
|
||||
result.distinct_invariants[i] = &result.invariants[i];
|
||||
}
|
||||
result.solver = mps_context_new();
|
||||
mps_context_set_output_prec(result.solver, 20); // relative precision
|
||||
mps_context_set_output_goal(result.solver, MPS_OUTPUT_GOAL_APPROXIMATE);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
void free_global_data(struct global_data dat)
|
||||
{
|
||||
for(int i = 0; i < dat.n; i++)
|
||||
mat_clear(dat.matrices[i]);
|
||||
free(dat.matrices);
|
||||
for(int i = 0; i < dat.n; i++) {
|
||||
mpq_clear(dat.invariants[i].tr);
|
||||
mpq_clear(dat.invariants[i].trinv);
|
||||
}
|
||||
free(dat.invariants);
|
||||
free(dat.distinct_invariants);
|
||||
mps_context_free(dat.solver);
|
||||
}
|
||||
|
||||
static int compare_result(const void *a_, const void *b_)
|
||||
{
|
||||
int d = 0;
|
||||
|
||||
struct result **a = (struct result **)a_;
|
||||
struct result **b = (struct result **)b_;
|
||||
|
||||
d = mpq_cmp((*a)->tr,(*b)->tr);
|
||||
if(d == 0)
|
||||
d = mpq_cmp((*a)->trinv, (*b)->trinv);
|
||||
|
||||
return d;
|
||||
}
|
||||
|
||||
int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, double *eigenvalues)
|
||||
{
|
||||
mpq_t coeff1, coeff2, zero;
|
||||
cplx_t *roots;
|
||||
double radii[3];
|
||||
double *radii_p[3];
|
||||
mps_monomial_poly *poly;
|
||||
mps_boolean errors;
|
||||
int result = 0;
|
||||
|
||||
mpq_inits(coeff1, coeff2, zero, NULL);
|
||||
mpq_set(coeff1, trinv);
|
||||
mpq_sub(coeff2, zero, tr);
|
||||
|
||||
poly = mps_monomial_poly_new(solv, 3);
|
||||
mps_monomial_poly_set_coefficient_int(solv, poly, 0, -1, 0);
|
||||
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_int(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[i] = cplx_Re(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;
|
||||
}
|
||||
|
||||
void continued_fraction_approximation(mpq_t out, double in, int level)
|
||||
{
|
||||
mpq_t tmp;
|
||||
|
||||
if(in < 0) {
|
||||
mpq_init(tmp);
|
||||
mpq_set_ui(tmp, 0, 1);
|
||||
continued_fraction_approximation(out, -in, level);
|
||||
mpq_sub(out, tmp, out);
|
||||
mpq_clear(tmp);
|
||||
return;
|
||||
}
|
||||
|
||||
if(level == 0) {
|
||||
mpq_set_si(out, (signed long int)round(in), 1); // floor(in)
|
||||
} else {
|
||||
continued_fraction_approximation(out, 1/(in - floor(in)), level - 1);
|
||||
mpq_init(tmp);
|
||||
mpq_set_ui(tmp, 1, 1);
|
||||
mpq_div(out, tmp, out); // out -> 1/out
|
||||
mpq_set_si(tmp, (signed long int)in, 1); // floor(in)
|
||||
mpq_add(out, out, tmp);
|
||||
mpq_clear(tmp);
|
||||
}
|
||||
}
|
||||
|
||||
void quartic(mpq_t out, mpq_t in, int a, int b, int c, int d, int e)
|
||||
{
|
||||
mpq_t tmp;
|
||||
mpq_init(tmp);
|
||||
|
||||
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);
|
||||
|
||||
mpq_clear(tmp);
|
||||
}
|
||||
|
||||
// this version is only for the (4,4,4) group
|
||||
void initialize_triangle_generators(mat_workspace *ws, mat *gen, mpq_t m, mpq_t t)
|
||||
{
|
||||
mpq_t s,sinv,q,x,y;
|
||||
mpq_t zero, one, two;
|
||||
mpq_t tmp;
|
||||
|
||||
mpq_inits(s,sinv,q,x,y,zero,one,two,tmp,NULL);
|
||||
mpq_set_ui(zero, 0, 1);
|
||||
mpq_set_ui(one, 1, 1);
|
||||
mpq_set_ui(two, 2, 1);
|
||||
|
||||
// s = (1-m^2)/2m
|
||||
mpq_mul(s, m, m);
|
||||
mpq_sub(s, one, s);
|
||||
mpq_div(s, s, m);
|
||||
mpq_div(s, s, two);
|
||||
mpq_div(sinv, one, s);
|
||||
|
||||
// q = (1+m^2)/(1-m^2) = 2/(1-m^2) - 1
|
||||
mpq_mul(q, m, m);
|
||||
mpq_sub(q, one, q);
|
||||
mpq_div(q, two, q);
|
||||
mpq_sub(q, q, one);
|
||||
|
||||
// x = -tq, y = -q/t
|
||||
mpq_mul(x, q, t);
|
||||
mpq_sub(x, zero, x);
|
||||
mpq_div(y, q, t);
|
||||
mpq_sub(y, zero, y);
|
||||
|
||||
// q^2 = xy = 1 + 1/s^2
|
||||
// [ -s s*y 0]
|
||||
// [ -s*x s*x*y - 1/s 0]
|
||||
// [ -s*y s*y^2 - x 1]
|
||||
LOOP(i,3) {
|
||||
mat_zero(gen[i]);
|
||||
mpq_sub(tmp, zero, s);
|
||||
mat_set(gen[i%3], i%3, i%3, tmp);
|
||||
mpq_mul(tmp, s, y);
|
||||
mat_set(gen[i%3], i%3, (i+1)%3, tmp);
|
||||
mpq_mul(tmp, s, x);
|
||||
mpq_sub(tmp, zero, tmp);
|
||||
mat_set(gen[i%3], (i+1)%3, i%3, tmp);
|
||||
mpq_mul(tmp, s, x);
|
||||
mpq_mul(tmp, tmp, y);
|
||||
mpq_sub(tmp, tmp, sinv);
|
||||
mat_set(gen[i%3], (i+1)%3, (i+1)%3, tmp);
|
||||
mpq_mul(tmp, s, y);
|
||||
mpq_sub(tmp, zero, tmp);
|
||||
mat_set(gen[i%3], (i+2)%3, i%3, tmp);
|
||||
mpq_mul(tmp, s, y);
|
||||
mpq_mul(tmp, tmp, y);
|
||||
mpq_sub(tmp, tmp, x);
|
||||
mat_set(gen[i%3], (i+2)%3, (i+1)%3, tmp);
|
||||
mat_set(gen[i%3], (i+2)%3, (i+2)%3, one);
|
||||
}
|
||||
|
||||
LOOP(i,3) mat_pseudoinverse(ws, gen[i+3], gen[i]);
|
||||
|
||||
// debug output
|
||||
/*
|
||||
gmp_printf("m = %Qd, s = %Qd, t = %Qd, q = %Qd, x = %Qd, y = %Qd\n", m, s, t, q, x, y);
|
||||
mat_print(gen[0]);
|
||||
mat_print(gen[1]);
|
||||
mat_print(gen[2]);
|
||||
*/
|
||||
|
||||
mpq_inits(s,sinv,q,x,y,zero,one,two,tmp,NULL);
|
||||
}
|
||||
|
||||
char *print_word(groupelement_t *g, char *str)
|
||||
{
|
||||
int i = g->length - 1;
|
||||
|
||||
str[g->length] = 0;
|
||||
while(g->parent) {
|
||||
str[i--] = 'a' + g->letter;
|
||||
g = g->parent;
|
||||
}
|
||||
|
||||
return str;
|
||||
}
|
||||
|
||||
void enumerate(group_t *group, mat *matrices, mpq_t m, mpq_t t)
|
||||
{
|
||||
mat_workspace *ws;
|
||||
mat tmp;
|
||||
mat gen[6];
|
||||
char buf[100], buf2[100], buf3[100];
|
||||
|
||||
// allocate stuff
|
||||
ws = mat_workspace_init(3);
|
||||
for(int i = 0; i < 6; i++)
|
||||
mat_init(gen[i], 3);
|
||||
mat_init(tmp, 3);
|
||||
|
||||
initialize_triangle_generators(ws, gen, m, t);
|
||||
|
||||
mat_identity(matrices[0]);
|
||||
for(int i = 1; i < group->size; i++) {
|
||||
if(group->elements[i].length % 2 != 0)
|
||||
continue;
|
||||
if(!group->elements[i].inverse)
|
||||
continue;
|
||||
|
||||
int parent = group->elements[i].parent->id;
|
||||
int grandparent = group->elements[i].parent->parent->id;
|
||||
int letter;
|
||||
|
||||
if(group->elements[parent].letter == 1 && group->elements[i].letter == 2)
|
||||
letter = 0; // p = bc
|
||||
else if(group->elements[parent].letter == 2 && group->elements[i].letter == 0)
|
||||
letter = 1; // q = ca
|
||||
else if(group->elements[parent].letter == 0 && group->elements[i].letter == 1)
|
||||
letter = 2; // r = ab
|
||||
if(group->elements[parent].letter == 2 && group->elements[i].letter == 1)
|
||||
letter = 3; // p^{-1} = cb
|
||||
else if(group->elements[parent].letter == 0 && group->elements[i].letter == 2)
|
||||
letter = 4; // q^{-1} = ac
|
||||
else if(group->elements[parent].letter == 1 && group->elements[i].letter == 0)
|
||||
letter = 5; // r^{-1} = ba
|
||||
|
||||
mat_multiply(ws, matrices[i], matrices[grandparent], gen[letter]);
|
||||
}
|
||||
|
||||
// free stuff
|
||||
for(int i = 0; i < 6; i++)
|
||||
mat_clear(gen[i]);
|
||||
mat_clear(tmp);
|
||||
mat_workspace_clear(ws);
|
||||
}
|
||||
|
||||
|
||||
double compute_max_slope(struct global_data dat, mpq_t t, mpq_t m)
|
||||
{
|
||||
// mpq_set_ui(t, ttick, 100);
|
||||
// mpq_set_ui(m, mtick, 100); // 414/1000 ~ sqrt(2)-1 <-> s=1
|
||||
// s = (1-mpq_get_d(m)*mpq_get_d(m))/(2*mpq_get_d(m));
|
||||
|
||||
int n = 0;
|
||||
int nmax = dat.n;
|
||||
int nuniq;
|
||||
double max_slope;
|
||||
int retval;
|
||||
double evs[3];
|
||||
|
||||
group_t *group = dat.group;
|
||||
mat *matrices = dat.matrices;
|
||||
struct result *invariants = dat.invariants;
|
||||
struct result **distinct_invariants = dat.distinct_invariants;
|
||||
mps_context *solver = dat.solver;
|
||||
|
||||
// DEBUG("Compute matrices\n");
|
||||
enumerate(group, matrices, m, t);
|
||||
|
||||
// DEBUG("Compute traces\n");
|
||||
n = 0;
|
||||
for(int i = 0; i < nmax; i++) {
|
||||
if(group->elements[i].length % 2 != 0 || !group->elements[i].inverse)
|
||||
continue;
|
||||
|
||||
mat_trace(invariants[i].tr, matrices[i]);
|
||||
mat_trace(invariants[i].trinv, matrices[group->elements[i].inverse->id]);
|
||||
|
||||
distinct_invariants[n++] = &invariants[i];
|
||||
}
|
||||
|
||||
// DEBUG("Get unique traces\n");
|
||||
|
||||
qsort(distinct_invariants, n, sizeof(struct result*), compare_result);
|
||||
|
||||
nuniq = 0;
|
||||
for(int i = 0; i < n; i++) {
|
||||
if(i == 0 || compare_result(&distinct_invariants[i], &distinct_invariants[nuniq-1]) != 0)
|
||||
distinct_invariants[nuniq++] = distinct_invariants[i];
|
||||
}
|
||||
|
||||
max_slope = 0;
|
||||
int max_slope_index;
|
||||
|
||||
// DEBUG("Solve characteristic polynomials\n");
|
||||
for(int i = 0; i < nuniq; i++) {
|
||||
retval = solve_characteristic_polynomial(solver, distinct_invariants[i]->tr, distinct_invariants[i]->trinv, evs);
|
||||
if(retval == 1) {
|
||||
fprintf(stderr, "Error! Could not solve polynomial.\n");
|
||||
continue;
|
||||
} else if(retval == 2) {
|
||||
continue;
|
||||
}
|
||||
|
||||
if(fabs(evs[0]) < fabs(evs[1]))
|
||||
SWAP(double, evs[0], evs[1]);
|
||||
if(fabs(evs[1]) < fabs(evs[2]))
|
||||
SWAP(double, evs[1], evs[2]);
|
||||
if(fabs(evs[0]) < fabs(evs[1]))
|
||||
SWAP(double, evs[0], evs[1]);
|
||||
|
||||
double x = log(fabs(evs[0]));
|
||||
double y = -log(fabs(evs[2]));
|
||||
|
||||
if(y/x > max_slope && (x > 0.1 || y > 0.1)) {
|
||||
max_slope_index = distinct_invariants[i] - invariants;
|
||||
max_slope = y/x;
|
||||
}
|
||||
|
||||
// gmp_printf("%Qd %Qd %f %f %f\n", distinct_invariants[i]->tr, distinct_invariants[i]->trinv, x, y, y/x);
|
||||
}
|
||||
|
||||
return max_slope;
|
||||
}
|
||||
|
||||
void write_results_and_end(struct job *jobs, const char *outfile)
|
||||
{
|
||||
DEBUG("writing output and shutting down\n");
|
||||
|
||||
FILE *f = fopen(outfile, "w");
|
||||
for(int i = 0; i < NJOBS; i++)
|
||||
fprintf(f, "%d/%d %d/%d %.10f %.10f %.10f %.3f\n",
|
||||
jobs[i].tparam, TDIV, jobs[i].mparam, MDIV,
|
||||
(double)jobs[i].tparam/TDIV, (double)jobs[i].mparam/MDIV, jobs[i].max_slope,
|
||||
jobs[i].time);
|
||||
fclose(f);
|
||||
|
||||
for(int i = 1; i < world_size; i++)
|
||||
MPI_Send(NULL, 0, job_datatype, i, JOB_SHUTDOWN, MPI_COMM_WORLD);
|
||||
|
||||
}
|
||||
|
||||
void run_master_process(int nmax, const char *restart, const char *outfile)
|
||||
{
|
||||
int total_jobs = NJOBS;
|
||||
int completed = 0;
|
||||
int queue_jobs = MIN(total_jobs, 2*world_size);
|
||||
struct job current_job;
|
||||
MPI_Status status;
|
||||
FILE *f;
|
||||
int continuing = 1;
|
||||
int restartf;
|
||||
struct job *alljobs;
|
||||
struct job *current;
|
||||
|
||||
restartf = open(restart, O_RDWR);
|
||||
if(restartf == -1 && errno == ENOENT) {
|
||||
restartf = open(restart, O_RDWR | O_CREAT, 0666);
|
||||
continuing = 0;
|
||||
}
|
||||
if(restartf == -1) {
|
||||
DEBUG("error opening restart file: %s\n", strerror(errno));
|
||||
exit(1);
|
||||
}
|
||||
ftruncate(restartf, total_jobs*sizeof(struct job));
|
||||
alljobs = (struct job*) mmap(0, total_jobs*sizeof(struct job), PROT_READ | PROT_WRITE, MAP_SHARED, restartf, 0);
|
||||
if(alljobs == MAP_FAILED) {
|
||||
DEBUG("error mapping restart file: %s\n", strerror(errno));
|
||||
exit(1);
|
||||
}
|
||||
|
||||
if(continuing) {
|
||||
for(int i = 0; i < total_jobs; i++)
|
||||
if(alljobs[i].done)
|
||||
completed++;
|
||||
} else {
|
||||
for(int tparam = TFROM; tparam <= TTO; tparam++) {
|
||||
for(int mparam = MFROM; mparam <= MTO; mparam++) {
|
||||
alljobs[JOBNR(tparam,mparam)].tparam = tparam;
|
||||
alljobs[JOBNR(tparam,mparam)].mparam = mparam;
|
||||
alljobs[JOBNR(tparam,mparam)].done = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fsync(restartf);
|
||||
|
||||
if(continuing) {
|
||||
DEBUG("continuing from restart file, %d/%d jobs completed, %d nodes\n", completed, total_jobs, world_size);
|
||||
} else {
|
||||
DEBUG("starting from scratch, %d jobs, %d nodes\n", total_jobs, world_size);
|
||||
}
|
||||
|
||||
if(completed >= total_jobs)
|
||||
{
|
||||
write_results_and_end(alljobs, outfile);
|
||||
goto cleanup;
|
||||
}
|
||||
|
||||
// assign initial jobs
|
||||
current = alljobs-1;
|
||||
for(int i = 0; i < 2*world_size; i++) {
|
||||
do {
|
||||
current++;
|
||||
} while(current < alljobs + total_jobs && current->done);
|
||||
if(current >= alljobs + total_jobs) // all jobs are assigned
|
||||
break;
|
||||
MPI_Send(current, 1, job_datatype, 1 + i%(world_size-1), JOB_ORDER, MPI_COMM_WORLD);
|
||||
}
|
||||
|
||||
while(1) {
|
||||
MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
|
||||
if(status.MPI_TAG == JOB_RESULT) {
|
||||
MPI_Recv(¤t_job, 1, job_datatype, MPI_ANY_SOURCE, JOB_RESULT, MPI_COMM_WORLD, &status);
|
||||
completed++;
|
||||
|
||||
DEBUG("job (%d,%d) completed by instance %d in %f seconds, result = %.3f, %d/%d done\n",
|
||||
current_job.tparam, current_job.mparam,
|
||||
status.MPI_SOURCE, current_job.time, current_job.max_slope, completed, total_jobs);
|
||||
|
||||
int nr = JOBNR(current_job.tparam, current_job.mparam);
|
||||
memcpy(&alljobs[nr], ¤t_job, sizeof(struct job));
|
||||
alljobs[nr].done = 1;
|
||||
|
||||
if(completed % FLUSH_INTERVAL == 0)
|
||||
fsync(restartf);
|
||||
|
||||
// find the next unassigned job
|
||||
do {
|
||||
current++;
|
||||
} while(current < alljobs + total_jobs && current->done);
|
||||
|
||||
if(current < alljobs + total_jobs) {
|
||||
MPI_Send(current, 1, job_datatype, status.MPI_SOURCE, JOB_ORDER, MPI_COMM_WORLD);
|
||||
}
|
||||
|
||||
if(completed >= total_jobs) {
|
||||
write_results_and_end(alljobs, outfile);
|
||||
goto cleanup;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
cleanup:
|
||||
|
||||
munmap(alljobs, total_jobs*sizeof(struct job));
|
||||
close(restartf);
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
int name_len;
|
||||
|
||||
MPI_Status status;
|
||||
|
||||
mpq_t m, t;
|
||||
double s;
|
||||
struct job current_job;
|
||||
int nmax;
|
||||
double max_slope;
|
||||
struct global_data dat;
|
||||
double jobtime;
|
||||
|
||||
clock_gettime(CLOCK_REALTIME, &starttime);
|
||||
|
||||
if(argc < 4) {
|
||||
fprintf(stderr, "Usage: mpirun -n <nr> --hostfile <hostfile> %s <number of elements> <restartfile> <outfile>\n", argv[0]);
|
||||
return 0;
|
||||
}
|
||||
nmax = atoi(argv[1]);
|
||||
|
||||
MPI_Init(NULL, NULL);
|
||||
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
|
||||
MPI_Get_processor_name(processor_name, &name_len);
|
||||
|
||||
// DEBUG("instance %d/%d started on %s\n", world_rank, world_size, processor_name);
|
||||
|
||||
int blocklengths[2] = {3, 2};
|
||||
MPI_Datatype types[2] = {MPI_INT, MPI_DOUBLE};
|
||||
MPI_Aint displacements[2] = {(size_t)&((struct job*)0)->tparam, (size_t)&((struct job*)0)->max_slope};
|
||||
MPI_Type_create_struct(2, blocklengths, displacements, types, &job_datatype);
|
||||
MPI_Type_commit(&job_datatype);
|
||||
|
||||
if(world_rank == 0) { // master processor
|
||||
run_master_process(nmax, argv[2], argv[3]);
|
||||
MPI_Finalize();
|
||||
return 0;
|
||||
}
|
||||
|
||||
// DEBUG("Allocate & generate group\n");
|
||||
mpq_inits(m, t, NULL);
|
||||
dat = allocate_global_data(nmax);
|
||||
dat.group = coxeter_init_triangle(4, 4, 4, nmax);
|
||||
|
||||
// fprintf(stderr, "max word length = %d\n", dat.group->elements[nmax-1].length);
|
||||
|
||||
while(1) {
|
||||
MPI_Probe(0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
|
||||
// MPI_Recv(¤t_job, 1, job_datatype, 0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
|
||||
if(status.MPI_TAG == JOB_SHUTDOWN) {
|
||||
// DEBUG("instance %d shutting down\n", world_rank);
|
||||
break;
|
||||
}
|
||||
else if(status.MPI_TAG == JOB_ORDER) {
|
||||
MPI_Recv(¤t_job, 1, job_datatype, 0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
|
||||
DEBUG("instance %d starting order (%d,%d)\n", world_rank, current_job.tparam, current_job.mparam);
|
||||
|
||||
jobtime = -MPI_Wtime();
|
||||
|
||||
// do the actual work
|
||||
mpq_set_ui(t, current_job.tparam, TDIV);
|
||||
mpq_set_ui(m, current_job.mparam, MDIV);
|
||||
s = (1-mpq_get_d(m)*mpq_get_d(m))/(2*mpq_get_d(m));
|
||||
|
||||
max_slope = compute_max_slope(dat, t, m);
|
||||
|
||||
jobtime += MPI_Wtime();
|
||||
|
||||
// fprintf(stdout, "%.5f %.5f %.5f %f\n",
|
||||
// mpq_get_d(t), mpq_get_d(m), s, max_slope);
|
||||
current_job.max_slope = max_slope;
|
||||
current_job.time = jobtime;
|
||||
|
||||
DEBUG("instance %d finished order (%d,%d) in %f seconds\n", world_rank, current_job.tparam, current_job.mparam, jobtime);
|
||||
|
||||
MPI_Send(¤t_job, 1, job_datatype, 0, JOB_RESULT, MPI_COMM_WORLD);
|
||||
}
|
||||
}
|
||||
|
||||
// DEBUG("Clean up\n");
|
||||
coxeter_clear(dat.group);
|
||||
free_global_data(dat);
|
||||
mpq_clears(m, t, NULL);
|
||||
|
||||
MPI_Type_free(&job_datatype);
|
||||
MPI_Finalize();
|
||||
}
|
||||
@@ -1,173 +0,0 @@
|
||||
#include "coxeter.h"
|
||||
#include "linalg.h"
|
||||
#include "mat.h"
|
||||
#include "enumerate_triangle_group.h"
|
||||
|
||||
#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
|
||||
#define DEBUG(msg, ...)
|
||||
|
||||
double mpq_log(mpq_t m_op)
|
||||
{
|
||||
static double logB = log(ULONG_MAX);
|
||||
|
||||
// Undefined logs (should probably return NAN in second case?)
|
||||
if (mpz_get_ui(mpq_numref(m_op)) == 0 || mpz_sgn(mpq_numref(m_op)) < 0)
|
||||
return -INFINITY;
|
||||
|
||||
// Log of numerator
|
||||
double lognum = log(mpq_numref(m_op)->_mp_d[abs(mpq_numref(m_op)->_mp_size) - 1]);
|
||||
lognum += (abs(mpq_numref(m_op)->_mp_size)-1) * logB;
|
||||
|
||||
// Subtract log of denominator, if it exists
|
||||
if (abs(mpq_denref(m_op)->_mp_size) > 0)
|
||||
{
|
||||
lognum -= log(mpq_denref(m_op)->_mp_d[abs(mpq_denref(m_op)->_mp_size)-1]);
|
||||
lognum -= (abs(mpq_denref(m_op)->_mp_size)-1) * logB;
|
||||
}
|
||||
return lognum;
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
mpq_t m, t, s, q, tmp, tmp2;
|
||||
mat_workspace *ws;
|
||||
mat gen[6];
|
||||
mps_context *solver;
|
||||
mps_monomial_poly *poly;
|
||||
mat element, inverse;
|
||||
int letter1, letter2, letter;
|
||||
mpq_t tr, trinv;
|
||||
double x, y, slope;
|
||||
int retval;
|
||||
double evs[3];
|
||||
char buf[100];
|
||||
double max_slope = 0;
|
||||
int max_slope_index = 0;
|
||||
double min_slope = INFINITY;
|
||||
int min_slope_index = 0;
|
||||
char *env;
|
||||
int mode;
|
||||
|
||||
if(argc < 2) {
|
||||
fprintf(stderr,
|
||||
"Usage: %s <s> <q> <word1> <word2> ...\n"
|
||||
"Computes jordan slopes of a list of group elements for a fixed representation.\n"
|
||||
"s,q: representation in the Hitchin component, given as rational numbers, e.g. 2/7\n"
|
||||
"word1, word2, ...: elements in the triangle rotation group, given as reflection group words\n"
|
||||
"output: word - jordan slope pairs\n"
|
||||
"+max slope index, max slope value, max slope word, min slope index, min slope value, min slope word\n"
|
||||
"controlled by environment variable MAXIMUM=no/yes/only, default yes\n",
|
||||
argv[0]);
|
||||
exit(0);
|
||||
}
|
||||
|
||||
mpq_inits(m, t, s, q, tmp, tmp2, tr, trinv, NULL);
|
||||
ws = mat_workspace_init(3);
|
||||
for(int i = 0; i < 6; i++)
|
||||
mat_init(gen[i], 3);
|
||||
mat_init(element, 3);
|
||||
mat_init(inverse, 3);
|
||||
|
||||
solver = mps_context_new();
|
||||
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);
|
||||
|
||||
mpq_set_str(s, argv[1], 10);
|
||||
mpq_set_str(q, argv[2], 10);
|
||||
|
||||
env = getenv("MAXIMUM");
|
||||
if(!env || strcmp(env, "yes") == 0) {
|
||||
mode = 1; // yes
|
||||
} else if(strcmp(env, "no") == 0) {
|
||||
mode = 0; // no
|
||||
} else if(strcmp(env, "only") == 0) {
|
||||
mode = 2; // only
|
||||
}
|
||||
|
||||
for(int w = 0; w < argc - 3; w++) {
|
||||
initialize_triangle_generators(ws, gen, 6, 4, 3, s, q);
|
||||
|
||||
mat_identity(element);
|
||||
mat_identity(inverse);
|
||||
for(int k = 0; k < strlen(argv[w+3]); k+=2) {
|
||||
letter1 = argv[w+3][k] - 'a';
|
||||
letter2 = argv[w+3][k+1] - 'a';
|
||||
|
||||
if(letter1 == 1 && letter2 == 2)
|
||||
letter = 0; // p = bc
|
||||
else if(letter1 == 2 && letter2 == 0)
|
||||
letter = 1; // q = ca
|
||||
else if(letter1 == 0 && letter2 == 1)
|
||||
letter = 2; // r = ab
|
||||
else if(letter1 == 2 && letter2 == 1)
|
||||
letter = 3; // p^{-1} = cb
|
||||
else if(letter1 == 0 && letter2 == 2)
|
||||
letter = 4; // q^{-1} = ac
|
||||
else if(letter1 == 1 && letter2 == 0)
|
||||
letter = 5; // r^{-1} = ba
|
||||
|
||||
mat_multiply(ws, element, element, gen[letter]);
|
||||
mat_multiply(ws, inverse, gen[(letter+3)%6], inverse);
|
||||
}
|
||||
|
||||
mat_trace(tr, element);
|
||||
mat_trace(trinv, inverse);
|
||||
|
||||
retval = solve_characteristic_polynomial(solver, poly, tr, trinv, evs);
|
||||
if(retval == 1) {
|
||||
fprintf(stderr, "Error! Could not solve polynomial.\n");
|
||||
return 1;
|
||||
}
|
||||
|
||||
if(fabs(evs[0]) < fabs(evs[1]))
|
||||
SWAP(double, evs[0], evs[1]);
|
||||
if(fabs(evs[1]) < fabs(evs[2]))
|
||||
SWAP(double, evs[1], evs[2]);
|
||||
if(fabs(evs[0]) < fabs(evs[1]))
|
||||
SWAP(double, evs[0], evs[1]);
|
||||
|
||||
x = log(fabs(evs[0]));
|
||||
y = -log(fabs(evs[2]));
|
||||
|
||||
if(x > DBL_MAX || y > DBL_MAX) {
|
||||
mpq_abs(tmp, tr);
|
||||
mpq_abs(tmp2, trinv);
|
||||
slope = mpq_log(tmp)/mpq_log(tmp2);
|
||||
} else {
|
||||
slope = y/x;
|
||||
}
|
||||
|
||||
if(slope < 1)
|
||||
slope = 1/slope;
|
||||
|
||||
if(slope > max_slope) {
|
||||
max_slope = slope;
|
||||
max_slope_index = w;
|
||||
}
|
||||
|
||||
if(slope < min_slope) {
|
||||
min_slope = slope;
|
||||
min_slope_index = w;
|
||||
}
|
||||
|
||||
if(mode != 2) {
|
||||
// gmp_printf("%s %.9f %Qd %Qd\n", argv[w+3], slope, tr, trinv);
|
||||
gmp_printf("%s %.9f %.9f %.9f\n", argv[w+3], slope, x, y);
|
||||
}
|
||||
}
|
||||
|
||||
if(mode != 0)
|
||||
printf("%d %.9f %s %d %.9f %s\n",
|
||||
max_slope_index, max_slope, argv[max_slope_index+3],
|
||||
min_slope_index, min_slope, argv[min_slope_index+3]);
|
||||
fflush(stdout);
|
||||
|
||||
mpq_clears(m, t, s, q, tmp, tmp2, tr, trinv, NULL);
|
||||
mat_workspace_clear(ws);
|
||||
for(int i = 0; i < 6; i++)
|
||||
mat_clear(gen[i]);
|
||||
mat_clear(element);
|
||||
mat_clear(inverse);
|
||||
mps_context_free(solver);
|
||||
}
|
||||
Reference in new issue
Block a user