# lib/algcurves/src/Siegel
#
# algcurves[Siegel]
#
# Author: Bernard Deconinck (2001).
# Modifications by Mark van Hoeij
# This procedure takes as input a Riemann matrix B.
# Every such matrix has a decomposition as
# B=real part+i*T*transpose(T). Then an ellipsoid
# is determined by (1/2)T*transpose(T)=constant. For
# many applications it is beneficial to have this ellipsoid
# be as round as possible. An optimal algorithm to achieve this
# is known in one dimension, but not beyond that. The algorithm
# below, due to Siegel, comes up with a Riemann matrix, which
# corresponds to a much more spherical ellipsoid; it is not
# known in how far this construction is optimal.
# The output is a list with two entries: the first one is the
# new Riemann matrix. The second one gives the symplectic
# matrix which gives the transformation between the old and new
# Riemann matrices.
macro(
sortlattice = `RiemannTheta/sortlattice`,
Siegel = `algcurves/Siegel`,
rowdim=LinearAlgebra:-RowDimension,
transpose=LinearAlgebra:-Transpose,
cholesky=LinearAlgebra:-LUDecomposition,
col=LinearAlgebra:-Column,
linsolve=LinearAlgebra:-LinearSolve,
mu=LinearAlgebra:-Multiply,
diag=LinearAlgebra:-DiagonalMatrix
):
Siegel:=proc(B::Matrix(square))
description "Author: Bernard Deconinck, changes by Mark van Hoeij";
option remember,system;
local i,g,id,t,b,zip,ro,Y,T,L,U,X,x,t1,t2,t3,t4;
b := evalf(B);
if indets(b)<>{} then
error "entries should be complex constants"
fi;
g:=rowdim(b);
id := Matrix(g,'shape'='identity');
t := Matrix(2*g, 'shape'='identity');
zip := Matrix(g,g);
ro:=proc(a)
local i;
Matrix([seq(map(round,i),i=convert(a,listlist))]);
end:
do
b:=0.5*(b+transpose(b));
Y:=map(Im,b);
T:=transpose(cholesky(Y,'method'='Cholesky'));
L:=transpose(Matrix(sortlattice(
[seq(convert(col(T,i),list),i=1..g)])));
U:=ro(linsolve(T,L));
b:=mu(mu(transpose(U), b), U);
t:=mu(Matrix([[transpose(U),zip],[zip,U^(-1)]]),t);
X:=map(Re,b);
x:=ro(X);
b:=b-x;
t:=mu(Matrix( [[id,-x],[zip,id]] ), t);
if abs(b[1,1]) >= 1 then
break
fi;
t1:=diag(<0 ,seq(1,i=1..g-1)>);
t2:=diag(<-1,seq(0,i=1..g-1)>);
t3:=diag(<1 ,seq(0,i=1..g-1)>);
t4:=diag(<0 ,seq(1,i=1..g-1)>);
t:=mu(Matrix([[t1,t2],[t3,t4]]),t);
b:=mu(mu(t1,b)+t2, linsolve(mu(t3,b)+t4,id))
od;
eval( [b,t] )
end:
#savelib('Siegel'):