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