(*:Name: NumberTheory`SiegelTheta` *)

(*:Title: Numerical Evaluation of Siegel's Theta Function *)

(*:Author: Jerry B. Keiper *)

(*:Summary:
This package provides for the numerical evaluation of
Siegel's theta function.
*)

(*:Context: NumberTheory`SiegelTheta`  *)

(*:Package Version: 1.0 *)

(* :Copyright: Copyright 1994-1999,  Wolfram Research, Inc.
*)

(* :History:
	Originally by Jerry B. Keiper, June 1994.
*)

(*:Keywords: theta functions, Siegel's theta function *)

(*:Source:
	Siegel, Topics in complex function theory, Vol II,
	Wiley-Interscience, p 163.
*)

(*:Mathematica Version: 2.0 *)

(*:Limitation:
	Multi-dimensional summation for high dimensions is inherently
	slow.  The sum in one dimension can be dealt with by looking
	at how the exponent Pi I (t . z . t + 2 t . s) changes.  But
	the trick doesn't seem to work for the other dimensions because
	there is coupling.
*)

BeginPackage["NumberTheory`SiegelTheta`"]

Unprotect[SiegelTheta];

SiegelTheta::usage = 
"SiegelTheta[Z, s] gives Siegel's theta function.  Z must be a square
complex matrix, Z == X + I Y, where Y is positive definite, and s must
be a complex vector compatible with Z."

Begin["`Private`"]

SiegelTheta[z_ /; MatrixQ[z, NumericQ], s_ /; VectorQ[s, NumericQ]] :=
    Module[{ans, len, p},
	ans /; (len = Length[s]; {len, len} === Dimensions[z] &&
		positivedefiniteQ[z, Im[z]] &&
		(p = Min[Precision[z], Precision[s]]) < Infinity &&
		NumberQ[ans = SiegelTheta0[z, s, 0, Table[0,{len}], 1, len, p]])
	]

SiegelTheta::npd =
"The imaginary part of the matrix `1` is not positive definite."

positivedefiniteQ[z_, y_] :=
	If[Det[y] > 0, True, Message[SiegelTheta::npd, z]; False] &&
	(Length[y] == 1 || positivedefiniteQ[z, Drop[#, -1]& /@ Drop[y, -1]]);

SiegelTheta0[z_, s_, ss_, t_, len_, len_, p_] :=
    Block[{oldsum = 0, sum, term, csp, csm, lsp=1, lsm=1, cq, lq},
	sum = ss + (term = E^(Pi I N[t . z . t + 2 t . s, p]));
	csp =  E^(Pi I N[z[[-1]].t + t.(Transpose[z][[-1]])+2 s[[-1]], p]);
	csm =  1/csp;
	lq = E^(Pi I N[z[[-1,-1]], p]);
	cq = lq^2;
	While[sum - oldsum != 0,
		oldsum = sum;
		term *= lq;
		lq *= cq;
		lsp *= csp;
		lsm *= csm;
		sum += term*(lsp + lsm);
		];
	sum
	];

SiegelTheta0[z_, s_, ss_, t_, i_, len_, p_] :=
    Block[{oldsum = 0, sum, tt = t},
	sum = SiegelTheta0[z, s, ss, tt, i+1, len, p];
	While[sum - oldsum != 0,
		oldsum = sum;
		tt[[i]]++;
		sum = SiegelTheta0[z, s, sum, tt, i+1, len, p];
		];
	oldsum = 0;
	tt = t;
	While[sum - oldsum != 0,
		oldsum = sum;
		tt[[i]]--;
		sum = SiegelTheta0[z, s, sum, tt, i+1, len, p];
		];
	sum
	];

End[] (* `Private` *)

Protect[SiegelTheta];

EndPackage[] (* NumberTheory`SiegelTheta` *)

(*

This is for brute-force checking:

test[n_] :=
    Module[{z, s, t, tt, k},
	z = Table[Random[Complex], {n}, {n}] + I IdentityMatrix[n];
	s = Table[Random[Complex], {n}];
	tt = Table[t[k], {k, n}];
	Sum[E^(Pi I N[tt . z . tt + 2 tt . s]),
		Evaluate[Sequence @@ Table[{t[k], -5, 5}, {k, n}]]]/
	    SiegelTheta[z,s] - 1
	]

*)
