# lib/src/RiemannTheta
#
# RiemannTheta
#
# Author: Bernard Deconinck  (2001).
# Modifications by Mark van Hoeij.
 
macro(
	bil = LinearAlgebra:-BilinearForm,
	dot = LinearAlgebra:-DotProduct,

	RE = `RiemannTheta/RE`,
	IM = `RiemannTheta/IM`,
	boundingellipsoid = `RiemannTheta/boundingellipsoid`,
	findvectors = `RiemannTheta/findvectors`,
	sortlattice = `RiemannTheta/sortlattice`,
	theta = RiemannTheta,
	doit = `RiemannTheta/doit`,
	make_proc = `RiemannTheta/make_proc`,
	finitesum = `RiemannTheta/finitesum`
):

RE := proc(b)
	local a;
	a:=expand(b);
	if type(a,`+`) then
		map(Re,a)
	else
		Re(a)
	fi
end:
IM := proc(b)
	local a;
	a:=expand(b);
	if type(a,`+`) then
		map(Im,a)
	else
		Im(a)
	fi
end:

`evalf/RiemannTheta` := proc()
	theta(args, `if`(type(args[-1],positive),NULL,Float(10,-Digits)))
end:
	

# Recursive procedure to find all integer vectors inside an ellipsoid. 
boundingellipsoid := proc(g, ch, shift, R)
local hc;
    hc := LinearAlgebra:-MatrixInverse(ch);
    {findvectors(g, ch, hc, Vector(shift), R, [])}
end:

findvectors := proc(g, ch, hc, shift, R, staart)
local a,b,k,newg,newch,newhc,i;
    a := ceil(evalf(-R/ch[g, g] + shift[g]));
    b := floor(evalf(R/ch[g, g] + shift[g]));
    k := a .. b;
    if g = 1 and a <= b then seq([i, op(staart)], i = k)
    elif 1 < g and a <= b then
	newg := g - 1;
	newch := LinearAlgebra:-SubMatrix(ch, 1 .. newg, 1 .. newg);
	newhc := LinearAlgebra:-SubMatrix(hc, 1 .. newg, 1 .. newg);
	seq(findvectors(newg, newch, newhc,
	    LinearAlgebra:-SubVector(shift, 1 .. newg) - (i - shift[g])*
	    LinearAlgebra:-MatrixVectorMultiply(newhc, LinearAlgebra:-SubVector(
	    LinearAlgebra:-Column(ch, g), 1 .. newg)),
	      sqrt(R^2 - ch[g, g]^2*(i - shift[g])^2), [i, op(staart)]),
	      i = k)
    else NULL
    fi
end:

# Sort lattice output w.r.t. length.
sortlattice := proc(L)
local i;
    sort(lattice(L), (s, t) -> evalb(add(i^2, i = s) <= add(i^2, i = t)))
end:


make_proc := proc(vv, g, omegar, omegai, S, dirs)
local V,n,S1,S2,v,i,ar,ai,ad,j,re,im,va,iv;
	V := [op(indets(vv,'name'))];
	n := nops(V);
	S1 :=	seq( Re(V[i]) = re[i], i = 1..n),
		seq( Im(V[i]) = im[i], i = 1..n),
		seq(    V[i]  = va[i], i = 1..n);
	S2 :=	seq( re[i] = Re(V[i]), i = 1..n),
		seq( im[i] = Im(V[i]), i = 1..n),
		seq( va[i] =    V[i] , i = 1..n);

	v := subs(S1, vv);
	i := Vector([seq(iv[j],j = 1..g)]);

	ar := i-v[1];
	ar := convert(
	   -.5*bil(ar, ar, omegar, conjugate = false), 
	   horner, [seq(iv[j],j = 1..g), op(indets(v,'name'))]);
	ai := convert( - .5*bil(i, i, omegai, conjugate = false)
	   + dot(i, v[2], conjugate = false),
	   horner, [seq(iv[j],j = 1..g), op(indets(v,'name'))]);

	ad := convert(mul(evalf(2*Pi*I)*add(iv[j]*i[j],j=1..g), i=dirs),
	   horner, [seq(iv[j],j = 1..g)]);

	S2 := {seq(`if`(has([ar,ai,ad],lhs(j)),j,NULL),j=[S2])};
	finitesum([ar,ai,ad,v[1]],S2,[seq(iv[j],j=1..g)],S,g,Digits)
end:


finitesum := proc(ar, S2, iv, S, g, dig)
local a,b,intshift,T,i,ed,e,j;
	Digits := dig;
	ed := 2;
	a := subs(S2,ar);
	if indets(a,'name') minus {op(iv)} <> {} then
		'procname(args)'
	else
		Digits := dig - 4;
		intshift := map(round,a[4]):
		a := a[1..3];
		T := 0;
		for i in S do
			b:=subs(seq(iv[j]=i[j]+intshift[j],j=1..g),a);
			e:=ilog10(b[2]);
			if e>ed then
				Digits:=Digits+e-ed;
				ed:=e;
				b:=subs(seq(iv[j]=i[j]+intshift[j],j=1..g),a)
			fi;
			T:=T+b[3]*`evalf/exp`(b[1])*
				(`evalf/cos`(b[2]) + `evalf/sin`(b[2])*I)
		od;
		T
	fi
end:


# Given a Riemann matrix B of dimension g, and a complex vector z of
# dimension g, calculate theta(z,B), the theta function associated with B, 
# evaluated at the point z, with an absolute error of epsilon.
# Note: this procedure does not require the use of Siegel to
# "well-condition" the Riemann matrix. Of course, doing this does result
# in faster computations.
# The output is a list of two elements [a,b]. Then theta=exp(a)*b. The
# exponential growth factor (first entry a) usually doesn't matter for
# applications, because one often takes quotients of Theta functions
# so we factor it out and return it seperately

theta := proc(zz, BB::Matrix, dirs, epsilon)
local B,z,i, old_dig, il;
    option `Copyright (c) 2001 Waterloo Maple Inc. All rights reserved.\
  Authors: B. Deconinck and M. van Hoeij.`;
    if nargs=2 or type(dirs,positive) then
	#
	# nargs < 2 --> error
	# nargs = 2 or argument dirs is
	# missing   --> add the argument dirs = []
	#
	return procname(zz,BB,[],args[3..-1])
    elif nargs<4 then
	#
	# epsilon missing --> unevaluated.
	# Note: evalf applied on this unevaluated expression
	# will insert epsilon so then the procedure should
	# evaluate.
	#
	return 'procname(args)'
    elif not type(dirs,list) then
	error dirs, "should be a list of directional derivatives"
    elif not type(epsilon,positive) then
	error "epsilon should be positive"
    fi;
    old_dig:=Digits;
    il := ilog10(epsilon);
    Digits := max(10, 10 - il, Digits+5);
    B := evalf(BB);
    z := evalf(2*I*Pi)*evalf(Vector(zz));
    if indets([B,dirs]) = {} then
	i:='output'='list';
	z:=doit(z,B,dirs,epsilon,Digits,op({args[5..-1]} minus {i}));
	if not member(i,[args]) then
		z:=exp(z[1])*z[2]
	fi;
	evalf(z, max(old_dig, 2-il))
    else
	'procname(args)'
    fi
end:

doit := proc(z,B,dirs,epsilon,dig)
local g,omega,omegar,omegai,omegarin,ch,v,r,nd,normks,k,normomegarin,
normx,d,bd,g2,x,R,i,y,ev,shift,intshift,leftshift,S,j,dch;
	Digits:=dig;
	g := LinearAlgebra:-RowDimension(B);
	if g <> LinearAlgebra:-Dimension(z) then error
	    "Dimensions of first and second argument are incompatible"
	fi;
	# omega = -2*I*Pi*B. Now B and omega should be symmetric because
	# for the Cholesky decomposition. Therefore, in the next
	# command we only work with the symmetric part of B, presumably
	# this only results in a Riemann matrix which is more accurate.
	# We switch to the normalization used in Matthias Heil's thesis.
	# Most of what follows uses the notation from this thesis. 
	omega := evalf(-I*Pi)*(B + LinearAlgebra:-Transpose(B));

	# omega := convert(omega,listlist):
	# omegar := Matrix([seq(map(RE,i),i=omega)]);
	# omegai := Matrix([seq(map(IM,i),i=omega)]);
	#
	# The following is a bit cleaner:
	#
	omegar := rtable(map(Re,omega),datatype=float);
	omegai := rtable(map(Im,omega),datatype=float);

	omegarin := LinearAlgebra:-MatrixInverse(omegar);
	ch := evalf(sqrt(1/2))*LinearAlgebra:-Transpose(
	    LinearAlgebra:-LUDecomposition(omegar, 'method' = 'Cholesky'));
	# We need the shortest lattice point from the lattice defined
	# by the columns of ch. A good estimate (in many cases, the best) is 
	# provided by the LLL algorithm, which is implemented in the lattice 
	# program. 
	v := sortlattice(
	    [seq(convert(LinearAlgebra:-Column(ch, i), list), i=1..g)])[1];
	r := sqrt(add(v[i]^2, i = 1 .. g));

	nd:=nops(dirs);
	normks:=[seq(evalf(sqrt(add(abs(k[i])^2,i=1..g))),k=dirs)];
     if nd>0 then
	for k in normks do if k=0 then
	  error "directions need to be nonzero vectors"
	fi od;
	normomegarin:=convert(LinearAlgebra:-Eigenvalues(
	  omegarin),list);
	normomegarin:=max(op(map(Re,normomegarin)));
	if nargs=6 and type(args[nargs]*1.0,float)
	then
		if args[6]>1 then normx:=args[6] fi
	else
		normx:=5/normomegarin
	fi;
	d:=1/min(seq(ch[i,i],i=1..g));
     fi;
	bd:=evalf(sqrt(g + 2*nd + sqrt(g^2 + 8*nd))/2 + r/2);
	g2:=g/2;
	# Compute the scaling factor R for the ellipsoid.
	dch := abs(LinearAlgebra:-Determinant(ch));
	x:=subs(((dummy-r/2)^2)^(1/2)=(dummy-r/2),evalf(add(
	  (normomegarin*normx)^(nd-i)*2^(i*(nd-3/2))*d^i*GAMMA(g2+i/2,(dummy-r/2)^2)
	 ,i=0..nd)*Pi^g2/GAMMA(g2)/dch*convert(normks,`*`)-epsilon));
	if evalf(subs(dummy=bd,x)) <= 0
	then
		R := bd
        else
		i:=bd+1.0;
		while evalf(subs(dummy=i,x))>0 do i:=i+1 od:
		while i-bd > 0.01 do
		  y:=(bd+i)*0.5;
		  if evalf(subs(dummy=y,x))>0 then bd:=y else i:=y fi
		od:
		R:=i
	fi;
	ev := evalb( indets(z) = {} );
	x := map(RE,z);
	y := map(IM,z);
	# shift is almost what is required to complete the square in the
	# exponent of the theta function. 
	shift := LinearAlgebra:-MatrixVectorMultiply(omegarin, x);
	if ev then
	    intshift := map(round, shift);
	    leftshift := shift - intshift
	else
	    leftshift := Vector(g,-0.5)
	fi;
	# The integer vectors inside the ellipsoid.
	S := boundingellipsoid(g, ch, leftshift, R);
	if not ev then for i to g do
	    S := S union {seq(subsop(i=j[i]+1,j),j=S)}
	od fi;
	evalf([expand(.5*dot(x, shift, conjugate = false)),
	  make_proc([shift,y], g, omegar, omegai, S ,dirs)])
end:

#savelib('RE','IM','`evalf/RiemannTheta`','boundingellipsoid',\
 'findvectors','sortlattice','make_proc','finitesum','theta','doit'):