# 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'):