# HG changeset patch # User Nick Alexander # Date 1244180879 25200 # Node ID d6e4600f23539bc7c604c6545d6ef898eaa66ce4 # Parent 341d08e697367f1a6b3ff76e7d154c8c7843b7fc [mq]: riemann-theta-1.patch diff -r 341d08e69736 -r d6e4600f2353 sage/functions/all.py --- a/sage/functions/all.py Fri Jun 19 16:06:12 2009 -0700 +++ b/sage/functions/all.py Thu Jun 04 22:47:59 2009 -0700 @@ -55,4 +55,10 @@ from spike_function import spike_function +from riemann_theta import (elliptic_theta, elliptic_theta_prime, + siegel_theta, siegel_theta_with_char, + shimura_phi, shimura_phi_with_char) +from complex_lattice import (ComplexLattice) +from klein_p import KleinPFunction # XXX + from prime_pi import prime_pi diff -r 341d08e69736 -r d6e4600f2353 sage/functions/complex_lattice.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/complex_lattice.py Thu Jun 04 22:47:59 2009 -0700 @@ -0,0 +1,110 @@ +# -*- sage -*- + +import sage.matrix.matrix +from sage.rings.all import ZZ, QQ, ComplexField +from riemann_theta import invert +from sage.matrix.all import matrix +from sage.modules.all import vector + +class ComplexLattice: + r""" + sage: from sage.functions.complex_lattice import ComplexLattice + + sage: tau = matrix(CDF, 2, 2, [I, 0, 0, I]) + sage: ComplexLattice(tau) + Lattice spanned by columns of + [1.0*I 0 1.0 0] + [ 0 1.0*I 0 1.0] + sage: L = ComplexLattice(tau, 0*tau + 2); L + Lattice spanned by columns of + [1.0*I 0 2.0 0] + [ 0 1.0*I 0 2.0] + sage: L._tau + [0.5*I 0] + [ 0 0.5*I] + + sage: L = ComplexLattice(tau.augment(identity_matrix(2))); L + Lattice spanned by columns of + [1.0*I 0 1.0 0] + [ 0 1.0*I 0 1.0] + sage: ComplexLattice(L) + Lattice spanned by columns of + [1.0*I 0 1.0 0] + [ 0 1.0*I 0 1.0] + + sage: tau3 = matrix(ComplexField(100), 3, 6, [[I, 0, 0, 1, 1, 0], [0, 2*I, 0, 0, 1, 1], [0, 0, 3*I, 1, 1, 1]]) + sage: L = ComplexLattice(tau3); L + Lattice spanned by columns of + [1.0*I 0 0 1.0 1.0 0] + [ 0 2.0*I 0 0 1.0 1.0] + [ 0 0 3.0*I 1.0 1.0 1.0] + """ + + def __init__(self, *inputs): + if len(inputs) == 1 and isinstance(inputs[0], ComplexLattice): + w1, w2 = inputs[0]._w1, inputs[0]._w2 + elif len(inputs) == 1 and sage.matrix.matrix.is_Matrix(inputs[0]): + w1 = inputs[0] + if w1.ncols() == w1.nrows(): + w2 = 0*w1 + ZZ(1) + elif w1.ncols() == 2*w1.nrows(): + w1, w2 = w1.submatrix(0, 0, w1.nrows(), w1.nrows()), w1.submatrix(0, w1.nrows(), w1.nrows(), w1.nrows()) + else: + w2 = None # will raise an error later + elif len(inputs) == 2 and sage.matrix.matrix.is_Matrix(inputs[0]) and sage.matrix.matrix.is_Matrix(inputs[1]): + w1, w2 = inputs + + if not (sage.matrix.matrix.is_Matrix(w1) and + sage.matrix.matrix.is_Matrix(w2) and + w1.is_square() and + w2.is_square() and + w1.nrows() == w2.nrows()): + raise ValueError, "Must give a lattice to ComplexLattice" + + self._tau = invert(w2) * w1 + self._C = self._tau.column_space() + self._g = self._tau.nrows() + assert self._tau.ncols() == self._g + self._Omega = w1.augment(w2) + self._w1 = w1 + self._w2 = w2 + + def rational_torsion_iterator(self, n): + n = ZZ(n) + assert n >= 1 + if n == 1: + yield vector(QQ, [0]*self._g) + raise StopIteration + v = vector(QQ, range(n))/n + + from sage.combinat.all import CartesianProduct + for w in CartesianProduct(*([v]*(2*self._g))): + yield vector(QQ, w) + + def rational_torsion(self, n): + return list(self.rational_torsion_iterator(n)) + + def torsion_iterator(self, n): + for w in self.rational_torsion_iterator(n): + yield self._Omega*vector(w) + + def torsion(self, n): + return list(self.torsion_iterator(n)) + + def ngens(self): + return 2*self._g + + def gens(self): + return self._Omega.columns() + + def periods(self): + return self._Omega.columns() + + def reduced(self): + return ComplexLattice(self._tau) + + def __repr__(self, prec=10): + return "Lattice spanned by columns of \n%s" % repr(self._Omega.change_ring(ComplexField(prec))) + + def acted_on_by(self, M): + return ComplexLattice(self._Omega * M.transpose()) diff -r 341d08e69736 -r d6e4600f2353 sage/functions/klein_p.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/klein_p.py Thu Jun 04 22:47:59 2009 -0700 @@ -0,0 +1,381 @@ +from sage.rings.all import ZZ, QQ, CC, CDF, RR, RDF, ComplexField, PolynomialRing +from sage.misc.cachefunc import cached_method +from sage.functions.all import ComplexLattice +from sage.matrix.all import matrix +from sage.modules.all import vector +from sage.functions.riemann_theta import invert, RiemannTheta +from sage.misc.misc import verbose +from sage.misc.flatten import flatten + +def myvar(*args, **kwds): + import sage.symbolic.ring + return sage.symbolic.ring.var(*args) # ignore kwds for now + +def symplectic_blocks(A): + n = A.nrows() // 2 + Gamma = A.copy() + Gamma.subdivide(n, n) + a, b, c, d = Gamma.subdivision(0, 0), Gamma.subdivision(0, 1), Gamma.subdivision(1, 0), Gamma.subdivision(1, 1) + return a, b, c, d + +def act_on_char(v, M): + a, b, c, d = symplectic_blocks(M) + B = (ZZ(1)/2*a.transpose()*c) + bs = [ B[i, i] for i in range(B.nrows()) ] + C = (ZZ(1)/2*b.transpose()*d) + cs = [ C[i, i] for i in range(B.nrows()) ] + u = M.transpose()*vector(v) + vector(bs + cs) + return u + +class KleinPFunction: + r""" + EXAMPLES:: + sage: Omega = matrix(CDF, 2, 4, [[ -4.7405 - 2.1033*I, -2.2680 - 0.73692*I, -1.4017 - 0.45544*I, 0.86631 + 1.1924*I], [ -4.6218 + 2.6971*I, 0.50844 - 0.69981*I, -0.82267 + 1.1323*I, -1.3311 - 0.43250*I]]) + + sage: c = [(1/2, 0), (1/2, 0)] + sage: T = KleinPFunction(c, Omega) + + The Klein $p$-functions are meromorphic and periodic on the given + lattice: + + sage: u = Omega * vector([1/2, -3/5, 4/7, 0]) + sage: T(u) + -2.64054775371 + 0.354328332569*I + sage: [ T(u + Omega*v) for v in (QQ^4).basis() ] + [-2.64052681336 + 0.354348633828*I, + -2.64053911586 + 0.354315398323*I, + -2.64054775371 + 0.354328332569*I, + -2.64054775371 + 0.354328332569*I] + + sage: list_of_derivs = [(0, 0, 0), (0, 1, 0), (1, 1)] + sage: u = Omega * vector([-1/9, -1/5, 4/7, 2/3]) + sage: T.values_at_point(u, list_of_derivs) + (-88.0400127429 - 185.610596981*I, 61.5755815636 - 52.9059237104*I, -2.80859789488 - 5.44145586052*I) + + sage: matrix(4, 3, [ T.values_at_point(u + Omega*v, list_of_derivs) for v in (QQ^4).basis() ]) + [ -88.039294736 - 185.627174678*I 61.5810988926 - 52.9057292964*I -2.80850414111 - 5.44148602994*I] + [-88.0474112895 - 185.640416204*I 61.5850902871 - 52.9132624176*I -2.80868603625 - 5.44190092188*I] + [-88.0400127429 - 185.610596981*I 61.5755815636 - 52.9059237104*I -2.80859789488 - 5.44145586052*I] + [-88.0400127429 - 185.610596981*I 61.5755815636 - 52.9059237104*I -2.80859789488 - 5.44145586052*I] + + The Klein $p$-functions as theta functions: + + sage: Omega = matrix(ComplexField(200), 2, 4, [[ -4.7405 - 2.1033*I, -2.2680 - 0.73692*I, -1.4017 - 0.45544*I, 0.86631 + 1.1924*I], [ -4.6218 + 2.6971*I, 0.50844 - 0.69981*I, -0.82267 + 1.1323*I, -1.3311 - 0.43250*I]]) + sage: u = Omega * vector([-1/9, -1/5, 4/7, 2/3]) + + sage: c = [(1/2, 0), (1/2, 0)] + sage: T = KleinPFunction(c, Omega) + + sage: w2i = sage.functions.riemann_theta.invert(T._lattice._w2) + + sage: F0 = lambda(u): (siegel_theta(T._lattice._tau, w2i*(u + T._lattice._w1*T._r + T._lattice._w2*T._s)))^2 + sage: v0 = F0(u) + sage: vs0 = [ F0(u + Omega*v)/v0 for v in (QQ^4).basis() ]; vs0 + [-16487.492394930226727441992580127569369598216846240721124250 + 479.19016920898539453431220364948886116325256409988519312547*I, + -75.383912252238665796461778245260453306538506293222785652535 - 199.03073804363052892930228798600037139096287262204524583069*I, + 1.0000000000000000000000000000000000000000000000000000000000 - 9.3345229167917125607160960806701863608853782530817506996517e-61*I, + 1.0000000000000000000000000000000000000000000000000000000000 - 4.6672614583958562803580480403350931804426891265408753498258e-59*I] + + sage: F1 = lambda(u): (T(u, theta=2)) + sage: v1 = F1(u) + sage: vs1 = [ F1(u + Omega*v)/v1 for v in (QQ^4).basis() ]; vs1 + [-16488.280681303702097631410738992696906860863450048519917067 + 478.77207179245557192756461248613786287191381296578412135710*I, + -75.386885202833233902841944316055644133621700984262144343630 - 199.05177602257541797724575926285400131223493890326758688094*I, + 1.0000000000000000000000000000000000000000000000000000000000 - 3.1115076389305708535720320268900621202951260843605835665506e-61*I, + 1.0000000000000000000000000000000000000000000000000000000000 + 5.9118645139680846217868608510911180285607395602851087764461e-60*I] + + sage: v1 + -3.6933728385553929940573381180805251112154357177580838632828 + 2.7371941934098641736456943872613357930208072935953309441321*I + """ + def __init__(self, rs, lattice): + self._lattice = ComplexLattice(lattice) + self._z = self._lattice._tau + self._g = self._z.nrows() # XXX only works for g = 2 anyways! + self._R = self._z.base_ring() + II = self._R.gen() + PI = self._R.pi() + self._CS = self._z.column_space() + self._r, self._s = map(self._CS, rs) + + u0, u1 = myvar('u0, u1', ns=1) + u = vector([u0, u1]) + v = u * invert(self._lattice._w2).transpose() + + self._Z = self._z - self._z.apply_map(lambda t: t.conjugate()) + self._Z = invert(self._Z) + self._scale = (II * PI * (v * self._Z * v)) + + import sage.symbolic.function + self._func = sage.symbolic.function.function('theta') + self._theta = self._func(*v) + self._RT = RiemannTheta(self._z) + + # cache with memory one + self._last_u = None + self._last_rs_diff = None + self._last_subs_dict = None + + def __repr__(self): + return "phi[...]( _, \n%s" % self._lattice._Omega.change_ring(ComplexField(10)) + + def reduced(self): + return KleinPFunction((self._r, self._s), self._lattice.reduced()) + + def acted_on_by(self, G): + d = act_on_char(vector(list(self._r) + list(self._s)), G) + Omega = self._lattice._Omega * G.transpose() + r = d[:self._g] + s = d[self._g:] + return KleinPFunction((r, s), Omega) + + def _derivs_upto(self, N): + r""" + EXAMPLES:: + sage: tau = matrix(CDF, 2, 2, [I, 0, 0, I]) + sage: T = KleinPFunction([(1/2, 0), (1/2, 0)], tau) + + sage: T._derivs_upto(0) + [] + sage: T._derivs_upto(1) + [(0), (1)] + sage: T._derivs_upto(2) + [(0), (1), (0, 0), (0, 1), (1, 1)] + sage: T._derivs_upto(3) + [(0), (1), (0, 0), (0, 1), (1, 1), (0, 0, 0), (0, 0, 1), (0, 1, 1), (1, 1, 1)] + """ + L = [] + for n in range(1, N+1): + xs = [0]*n + L.append(vector(xs)) + for m in reversed(range(0, n)): + xs[m] = 1 + L.append(vector(xs)) + return L + + def _rational_function_for_derivative(self, derivs, theta=0): + r""" + Calculate the rational function giving $p_{derivs}$ in terms of theta + derivatives. + + sage: tau = matrix(CDF, 2, 2, [I, 0, 0, I]) + sage: T = KleinPFunction([(1/2, 0), (1/2, 0)], tau) + sage: T._rational_function_for_derivative((0, 0)) + (((-3.14159265359)*X^2 + X_0^2 - X*X_00)/X^2, {(0,): 'X_0', (0, 0): 'X_00'}) + sage: T._rational_function_for_derivative((0, 1)) + ((-X*X_01 + X_0*X_1)/X^2, {(0, 1): 'X_01', (0,): 'X_0', (1,): 'X_1'}) + sage: T._rational_function_for_derivative((1, 1, 1)) + (((-2.0)*X_1^3 + 3.0*X*X_1*X_11 - X^2*X_111)/X^3, {(1, 1, 1): 'X_111', (1, 1): 'X_11', (1,): 'X_1'}) + + The theta options works: + + sage: T._rational_function_for_derivative((0, 0), theta=0) + (((-3.14159265359)*X^2 + X_0^2 - X*X_00)/X^2, {(0,): 'X_0', (0, 0): 'X_00'}) + sage: T._rational_function_for_derivative((0, 0), theta=1) + (((-3.14159265359)*X^2 + X_0^2 - X*X_00)/X, {(0,): 'X_0', (0, 0): 'X_00'}) + sage: T._rational_function_for_derivative((0, 0), theta=2) + ((-3.14159265359)*X^2 + X_0^2 - X*X_00, {(0,): 'X_0', (0, 0): 'X_00'}) + sage: T._rational_function_for_derivative((0, 0), theta=3) + ((-3.14159265359)*X^3 + X*X_0^2 - X^2*X_00, {(0,): 'X_0', (0, 0): 'X_00'}) + """ + assert len(derivs) >= 2 + + u0, u1 = myvar('u0, u1', ns=1) + u = vector([u0, u1]) + v = u * invert(self._lattice._w2).transpose() + + total_symb = -1 * (self._theta.log() + self._scale) + for deriv in derivs: + total_symb = total_symb.derivative(u[deriv]) + + # since v references the vars of u, need to be careful substituting + V0, V1 = myvar('V0, V1', ns=1) + D1 = {u0:V0, u1:V1} + D2 = {V0:v[0], V1:v[1]} + def patch_up(t): + return t.subs(D1).subs(D2) + + empty_subs_dict = {} + for deriv in self._derivs_upto(len(derivs)): + f = self._func(u0, u1) + for ind in deriv: + f = f.derivative(u[ind]) + f = patch_up(f) # very important! + nm = 'X_' + ''.join(map(repr, deriv)) + + # print "Checking %s for %s, replacing with %s.." % (total_symb, f, nm) + if total_symb.has(f): + total_symb = total_symb.subs(f == myvar(nm, ns=1)) + empty_subs_dict[tuple(deriv)] = nm + # print "Built empty_subs_dict" + + # these ones are always present + X = myvar('X', ns=1) + total_symb = total_symb.subs(self._theta.log() == myvar('logX', ns=1)) # important that log is substitutated + total_symb = total_symb.subs(self._theta == X) # before the value + total_symb = total_symb.subs(u[0] == myvar('u0', ns=1)) + total_symb = total_symb.subs(u[1] == myvar('u1', ns=1)) + + names = map(repr, total_symb.variables()) + names.remove('X') + names = ['X'] + names # make 'X' first name + + S = PolynomialRing(self._R, len(names), names=names) # self._R[*names] + XinS = S.gen(0) + # one last bit of hijinks, to avoid computing GCDs over inexact rings + + # usually we scale numerater and denominator by X**(len(derivs)) + total_symb = (X**len(derivs) * total_symb) + total_poly = (total_symb.expand())._polynomial_(S) + if len(derivs) - theta < 0: + total_poly *= XinS**( theta - len(derivs) ) + elif len(derivs) - theta > 0: + total_poly /= XinS**( len(derivs) - theta ) + + assert total_poly in S.fraction_field() + return (total_poly, empty_subs_dict) + + def _subs_dict_for(self, empty_subs_dict, u, rs_diff=None): + u_ = self._CS(u) + if rs_diff is None: + rs_diff = [(0, 0), (0, 0)] + + v_ = self._CS(invert(self._lattice._w2) * u_) + + rdiff, sdiff = map(self._CS, rs_diff) # weird way of thinking of this + v__ = v_ + self._z*rdiff + sdiff + + v__ = v__ + self._z*self._r + self._s # XXX this is *not enough* if we're not taking second or higher derivatives! + + items = empty_subs_dict.items() + derivs = [] + for inds, name in items: + deriv = [ (QQ**2).gen(ind) for ind in inds ] + derivs.append(deriv) + + # do it! + verbose("Calculating %s derivatives.." % len(derivs), level=1) + vals = self._RT.derivatives_at_point(v__, derivs) + verbose("Calculating %s derivatives.. DONE" % len(derivs), level=1) + X = vals.pop(0) + logX = X.log() + + subs_dict = {} + for i in range(len(items)): + inds, name = items[i] + subs_dict[name] = vals[i] + + subs_dict['X'] = X + subs_dict['logX'] = logX + + self._last_subs_dict = subs_dict + return subs_dict + + def values_at_point(self, u, list_of_derivs, rs_diff=None, theta=0): + r""" + EXAMPLES:: + sage: tau = matrix(CDF, 2, 2, [I, 1, 1, 2*I]) + sage: u = tau*vector([1/5, 5/7]) + vector([1/3, 1/4]) + + sage: P = KleinPFunction([(1/2, 0), (1/2, 1/2)], tau) + sage: P.values_at_point(u, [(1, 1), (0, 0)]) + (-4.42960579531 + 1.07101018644*I, -21.4773820675 - 10.4793989819*I) + sage: P(u, (0, 0)) + -21.4773820675 - 10.4793989819*I + + Here's an identity satisifed by the Klein $p$-functions: + + sage: tau = 1 + matrix(CDF, 2, 2, [5*I, 1, 1, I]) + sage: u = tau*vector([1/5, 5/7]) + vector([1/3, 1/4]) + sage: P = KleinPFunction([(1/2, 0), (1/2, 1/2)], tau) + + sage: X112, X222, X12, X122, X22 = P.values_at_point(u, [(0, 0, 1), (1, 1, 1), (0, 1), (0, 1, 1), (1, 1)]) + sage: v = X112 - X222*X12 + X122*X22 + sage: v.abs() < 1e-10 + True + + sage: X112, X222, X12, X122, X22 = P.values_at_point(2*u, [(0, 0, 1), (1, 1, 1), (0, 1), (0, 1, 1), (1, 1)]) + sage: v = X112 - X222*X12 + X122*X22 + sage: v.abs() < 1e-10 + True + """ + empty_subs_items = [] + polys = [] + for deriv in list_of_derivs: + poly, esd = self._rational_function_for_derivative(tuple(deriv), theta=theta) + polys.append(poly) + empty_subs_items += esd.items() # merge all dicts + empty_subs_dict = dict(empty_subs_items) + + subs_dict = self._subs_dict_for(empty_subs_dict, u, rs_diff=rs_diff) + vals = [ self._R(poly.subs(**subs_dict)) for poly in polys ] + return vector(vals) + + def __call__(self, u, derivs=[0, 0], rs_diff=None, theta=0): + vals = self.values_at_point(u, [derivs], rs_diff=rs_diff, theta=theta) + return vals[0] + +def generators_of_SP4Z(scale=1): + r""" + """ + from sage.matrix.all import block_diagonal_matrix, identity_matrix + + Gs = [] + for A in [ matrix(2, 2, [0, 1, -1, 0]), matrix(2, 2, [1, 1, 0, 1]) ]: + G = block_diagonal_matrix([A, ~A.transpose()]) + Gs.append(G) + + # the form [1, b, 0, 1], b symmetric + for B in [ matrix(2, 2, [1, 0, 0, 0]), matrix(2, 2, [0, 0, 0, 1]), matrix(2, 2, [0, 1, 1, 0]) ]: + # B *= scale # for GammaTheta + assert B.is_symmetric() + G = identity_matrix(4, 4) + G.set_block(0, 2, B) + # assert G in GammaTheta(2) + Gs.append(G) + + G = matrix(4, 4, [0, 0, 1, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, -1, 0, 0]) + Gs.append(G) + + assert all(G.det() == 1 for G in Gs) + # assert all(is_symplectic(G) for G in Gs) + return Gs + +def assert_klein_p_automorphy(n=5): + r""" + sage: from sage.functions.klein_p import assert_klein_p_automorphy + sage: assert_klein_p_automorphy(n=5) + Asserting that klein_p is automorphic.. + Asserting that klein_p is automorphic.. DONE + """ + print "Asserting that klein_p is automorphic.." + R = ComplexField(100) + TPII = 2 * R.gen() * R.pi() + import sage.functions.riemann_theta_tests + it = sage.functions.riemann_theta_tests.good_random_siegel_theta_element_iterator(2, R=R) + cnt = 0 + + list_of_derivs = [ (0, 0), (0, 1, 1), (1, 0), (0, 1) ] + + c = [(ZZ(0), ZZ(1)/2), (ZZ(0), ZZ(1)/2)] + while cnt < n: + z, u = it.next() + verbose("z = \n%s" % z, level=2) + verbose("u = \n%s" % u, level=2) + + T = KleinPFunction(c, z) + v1 = T.values_at_point(u, list_of_derivs) + if not all(1e-2 < ge.abs() < 1e20 for ge in v1): + continue + + failed = False + for G in generators_of_SP4Z(scale=1): + T2 = T.acted_on_by(G) + v2 = T2.values_at_point(u, list_of_derivs) # , rs_diff=rs_diff) + + e2 = (v1 - v2).norm() < 1e-6 + if not e2: + print "\nv1=%s\nv2=%s\ne2=%s\nT._last_subs_dict=%s" % (v1, v2, e2, T._last_subs_dict) + failed = True + assert not failed + cnt += 1 + print "Asserting that klein_p is automorphic.. DONE" diff -r 341d08e69736 -r d6e4600f2353 sage/functions/riemann_theta.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/riemann_theta.py Thu Jun 04 22:47:59 2009 -0700 @@ -0,0 +1,1337 @@ +r""" +Computing Riemann Theta Functions +================================= + +This module follows the paper "Computing Riemann Theta Functions" +listed below, hereafter referred to as [CRTF]. It is a more or less +direct translation of the accompanying Maple code. + +Let `g` be a positive integer, the *genus* or *dimension* of the Riemann theta +function. Let `H_g` denote the Siegel upper half space of dimension `g`, that +is the space of symmetric complex matrices whose imaginary parts are positive +definite. When `g = 1`, this is just the complex upper half plane. + +The Riemann theta function `\theta : CC^g \times H_g \to CC` is defined by +the infinite series + +.. math:: + + \theta(z, \Omega) = \sum_{n \in ZZ^g} \exp \left( 2 \pi i \left( + n^{t} \cdot \Omega \cdot n + n^t \cdot z \right) \right) . + +It is holomorphic in both `z` and `\Omega`. + +We write `z = x + i y` and `\Omega = X + i Y`. We denote by `[ V ]` the vector +with integer components closest to `V` and by `[[ V ]]` we denote `V - [V]`. + +REFERENCES: + +- [CRTF] Computing Riemann Theta Functions. Bernard Deconinck, Matthias Heil, + Alexander Bobenko, Mark van Hoeij and Markus Schmies. Mathematics + of Computation 73 (2004) 1417-1442. The paper is available at + http://www.amath.washington.edu/~bernard/papers/pdfs/computingtheta.pdf + and accompanying Maple code at http://www.math.fsu.edu/~hoeij/RiemannTheta/ . + +- http://en.wikipedia.org/wiki/Theta_function + +- http://mathworld.wolfram.com/JacobiThetaFunctions.html + +AUTHORS: + +- Nick Alexander (2009-03) +""" + +from sage.ext.fast_callable import fast_callable +from sage.rings.all import RDF, CDF, RealField, ComplexField, ZZ, QQ, GF +from sage.misc.misc import verbose +from sage.misc.misc_c import prod # sum is built in +from sage.modules.free_module_element import vector +from sage.matrix.constructor import matrix, identity_matrix +from sage.parallel.decorate import parallel + +@parallel() +def _global_theta_for_parallel(characteristics, RT, u): + return RT.siegel_theta_with_char(characteristics, u) + +def finite_sum_without_derivatives(shift, intshift, x, y, X, Y, S, domain): + r""" + Calculate the finite sum of [CRTF] Theorems 2, 4, and 6 when no + directional derivatives are wanted. + + .. note:: + + It is anticipated that this loop could be tightened in Cython + at some point in the future, so it has been separated into a + global function for ease of importing from a Cython module. + """ + g = X.nrows() + from sage.rings.all import PolynomialRing + if g > 1: + R = PolynomialRing(domain, 'u', g) + else: + R = PolynomialRing(domain, 'u') + a_ = vector(R.gens()) + + epart_ = -ZZ(1)/2 * (a_ + intshift - shift) * X * (a_ + intshift - shift) + epart = fast_callable(epart_, domain=domain) + + sincospart_ = (-ZZ(1)/2 * (a_ + intshift) * Y * (a_ + intshift) + (a_ + intshift) * y) + sincospart = fast_callable(sincospart_, domain=domain) + + Treal = 0 + Timag = 0 + for i in S: + ep = epart(*i).exp() + sc = sincospart(*i) + c = sc.cos() + s = sc.sin() + Treal += ep * c + Timag += ep * s + + return Treal, Timag + + +def finite_sum_with_many_derivatives(shift, intshift, x, y, X, Y, S, list_of_derivs, domain): + r""" + Calculate the finite sum of [CRTF] Theorems 2, 4, and 6 when + *many* directional derivatives are wanted. + + .. note:: + + Breaking things up into real and imaginary parts lets us use + special ``fast_callable`` interpreters over ``RDF`` and + ``RealField``, which are very fast. + + .. note:: + + It is anticipated that this loop could be tightened in Cython + at some point in the future, so it has been separated into a + global function for ease of importing from a Cython module. + """ + g = X.nrows() + from sage.rings.all import PolynomialRing + if g > 1: + R = PolynomialRing(domain, 'u', g) + else: + R = PolynomialRing(domain, 'u') + a_ = vector(R.gens()) + + epart_ = -ZZ(1)/2 * (a_ + intshift - shift) * X * (a_ + intshift - shift) + epart = fast_callable(epart_, domain=domain) + + sincospart_ = (-ZZ(1)/2 * (a_ + intshift) * Y * (a_ + intshift) + (a_ + intshift) * y) + sincospart = fast_callable(sincospart_, domain=domain) + + derivrs = [] + derivis = [] + for derivs in list_of_derivs: + derivr_ = R(prod(deriv.apply_map(lambda t: t.real()) * (a_ + intshift) for deriv in derivs)) + derivr = fast_callable(derivr_, domain=domain) + derivi_ = R(prod(deriv.apply_map(lambda t: t.imag()) * (a_ + intshift) for deriv in derivs)) + derivi = fast_callable(derivi_, domain=domain) + derivrs.append(derivr) + derivis.append(derivi) + + derivrs = [ fast_callable(R(1), domain=domain) ] + derivrs + derivis = [ fast_callable(R(0), domain=domain) ] + derivis + + # doesn't seem to make much difference + # return finite_sum_with_many_derivatives_inner(derivrs, derivis, S, epart, sincospart) + + Treal = [0] * len(derivrs) + Timag = [0] * len(derivrs) + rang = range(0, len(derivrs)) + + for i in S: + ep = epart(*i).exp() + sc = sincospart(*i) + c = ep * sc.cos() + s = ep * sc.sin() + for k in rang: + dpr = derivrs[k](*i) + dpi = derivis[k](*i) + Treal[k] += dpr * c - dpi * s + Timag[k] += dpi * c + dpr * s + + return zip(Treal, Timag) + +def finite_sum_with_derivatives(shift, intshift, x, y, X, Y, S, derivs, domain): + r""" + Calculate the finite sum of [CRTF] Theorems 2, 4, and 6 when + directional derivatives are wanted. + + .. note:: + + Breaking things up into real and imaginary parts lets us use + special ``fast_callable`` interpreters over ``RDF`` and + ``RealField``, which are very fast. + + .. note:: + + It is anticipated that this loop could be tightened in Cython + at some point in the future, so it has been separated into a + global function for ease of importing from a Cython module. + """ + g = X.nrows() + from sage.rings.all import PolynomialRing + if g > 1: + R = PolynomialRing(domain, 'u', g) + else: + R = PolynomialRing(domain, 'u') + a_ = vector(R.gens()) + + epart_ = -ZZ(1)/2 * (a_ + intshift - shift) * X * (a_ + intshift - shift) + epart = fast_callable(epart_, domain=domain) + + sincospart_ = (-ZZ(1)/2 * (a_ + intshift) * Y * (a_ + intshift) + (a_ + intshift) * y) + sincospart = fast_callable(sincospart_, domain=domain) + + derivr_ = R(prod(deriv.apply_map(lambda t: t.real()) * (a_ + intshift) for deriv in derivs)) + derivr = fast_callable(derivr_, domain=domain) + derivi_ = R(prod(deriv.apply_map(lambda t: t.imag()) * (a_ + intshift) for deriv in derivs)) + derivi = fast_callable(derivi_, domain=domain) + + Treal = 0 + Timag = 0 + for i in S: + dpr = derivr(*i) + dpi = derivi(*i) + ep = epart(*i).exp() + sc = sincospart(*i) + + c = ep * sc.cos() + s = ep * sc.sin() + Treal += dpr * c - dpi * s + Timag += dpi * c + dpr * s + + return Treal, Timag + +def finite_sum(shift, intshift, x, y, X, Y, S, derivs, domain): + if derivs: + return finite_sum_with_derivatives(shift, intshift, x, y, X, Y, S, derivs, domain) + else: + return finite_sum_without_derivatives(shift, intshift, x, y, X, Y, S, domain) + +def invert(M): + # sage can't invert inexact matrices reliably, so we do this + n = M.ncols() + return M.augment(M.parent().identity_matrix()).echelon_form().matrix_from_columns(range(n, 2*n)) + +def upper_gamma_1(a, b): + c = a.gamma_inc(b) + verbose("Computing gamma_inc(%s, %s) = %s" % (a, b, c), level=3) + return c + +import sympy.mpmath +def upper_gamma_2(a, b): + a = a.real() + b = b.real() + aaa = a.str() # truncate=False, no_sci=2) + bbb = b.str() # truncate=False, no_sci=2) + # print aaa + # print bbb + aa = sympy.mpmath.mpf( aaa, prec=a.parent().prec() ) + bb = sympy.mpmath.mpf( bbb, prec=b.parent().prec() ) + # print aa, aa.prec() + # print bb + s = repr( sympy.mpmath.upper_gamma(aa, bb) ) + i = s.find("'") + j = s.find("'", i+1) + return a.parent()(s[i+1:j]) + +def upper_gamma_3(a, b): + aa = sympy.mpmath.mpf( a.str(), prec=a.parent().prec() ) + bb = sympy.mpmath.mpf( b.str(), prec=b.parent().prec() ) + return a.parent()(float(sympy.mpmath.upper_gamma(aa, bb))) + +upper_gamma = upper_gamma_1 + +class RiemannTheta: + def __init__(self, Omega): + self._max_points_in_one_evaluation = 25000 + self._g = ZZ(Omega.nrows()) + assert self._g == Omega.ncols() + + self._ring = Omega.base_ring() # let the coercion mechanism deduce where to compute + + self._Omega = -self._ring.pi() * self._ring.gen() * (Omega + Omega.transpose()).change_ring(self._ring) + self._X = self._Omega.apply_map(lambda t: t.real()) + self._Y = self._Omega.apply_map(lambda t: t.imag()) + self._T = (self._ring(1)/2).real().sqrt() * self._X.cholesky_decomposition().transpose() + verbose("Cholesky decomposition is:", level=2) + verbose(self._T, level=2) + + # cached + self._no_derivs_radius = None + self._gradient_radius = None + self._hessian_radius = None + + self._Xin = invert(self._X) + self._Tin = invert(self._T) + + def lattice(self): + return ComplexLattice(self._Omega) + + def genus(self): + return self._g + + def integer_points(self, z, R): + r""" + The set `S_R` of [CRTF], (7). + + `S_R` is the set of integral points in the ellipsoid: + + .. math:: + + \left\{ n \in ZZ^g : \pi ( n + [[ Y^{-1} y ]] )^{t} \cdot + Y \cdot ( n + [[ Y^{-1} y ]] ) < R^2 \right\} . + + Since `Y` is positive definite it has a Cholesky decomposition `Y = T^t + T`. Letting `\Lambda` be the lattice of vectors `v(n), n \in ZZ^g` of + the form `v(n) = \sqrt{\pi} T (n + [[ Y^{-1} n ]])`, we have that + + .. math:: + + S_R = \left\{ v(n) \in \Lambda : || v(n) || < R \right\} . + + .. warning:: + + At times we work over ``CDF`` and ``RDF``, which have very + low precision (53 bits). This could be a problem when + given ill-conditioned input. In general computing theta + functions with such ill-conditioned will not be possible, + so we do not concern ourselves with this case. + + EXAMPLES: + + One-dimensional examples:: + + sage: I = CDF.gen() + sage: Omega = matrix(CDF, 1, 1, [I]) + sage: from sage.functions.riemann_theta import RiemannTheta + sage: RiemannTheta(Omega).integer_points(vector([I]), 2) + [[-1], [0], [1]] + sage: RiemannTheta(Omega).integer_points(vector([I]), 8.5) + [[-4], [-3], [-2], [-1], [0], [1], [2], [3], [4]] + sage: RiemannTheta(Omega).integer_points(vector([20.5 + I]), 8.5) + [[-4], [-3], [-2], [-1], [0], [1], [2], [3], [4], [5]] + + Two-dimensional examples:: + + sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I]) + sage: RiemannTheta(Omega).integer_points(vector([I, 1]), 2) + [[-1, 0], [0, 0], [1, 0], [0, 1]] + sage: RiemannTheta(Omega).integer_points(2*vector([I, 1]), 2) + [[-1, 0], [0, 0], [1, 0], [0, 1]] + sage: RiemannTheta(Omega).integer_points(2*vector([I + 100, 1]), 2.2) + [[-1, 0], [0, 0], [1, 0], [-1, 1], [0, 1]] + sage: len( RiemannTheta(Omega/4).integer_points(2*vector([I + 100, 1]), 2.2) ) + 20 + + Now we work with higher precision:: + + sage: Omega = matrix(ComplexField(100), 2, 2, [1 + I, 2, 2, I]) + sage: RiemannTheta(Omega).integer_points(vector([I, 1]), 2) + [[-1, 0], [0, 0], [1, 0], [0, 1]] + sage: RiemannTheta(Omega).integer_points(2*vector([I, 1]), 2) + [[-1, 0], [0, 0], [1, 0], [0, 1]] + sage: RiemannTheta(Omega).integer_points(2*vector([I + 100, 1]), 2.2) + [[-1, 0], [0, 0], [1, 0], [-1, 1], [0, 1]] + sage: len( RiemannTheta(Omega/4).integer_points(2*vector([I + 100, 1]), 2.2) ) + 20 + """ + R = self._ring(R).real() + z = self._Omega.column_space()(z) + x = z.apply_map(lambda t: t.real()) + y = z.apply_map(lambda t: t.imag()) + shift = self._Xin * x + intshift = shift.apply_map(lambda t: ZZ(t.round())) + leftshift = shift - intshift + verbose("Shift %s, Integer part of shift %s, difference is %s" % (shift, intshift, leftshift), level=2) + + T = self._T + Tin = self._Tin + + # internal helper function + def find_integer_points(g, shift, R, start): + t = (-R / T[g, g] + shift[g]) + a = (-R / T[g, g] + shift[g]).ceil() + b = (R / T[g, g] + shift[g]).floor() + if b - a > self._max_points_in_one_evaluation: + raise ValueError, "Evaluating theta function required too many summands (at least %s)" % (b-a) + + if not a <= b: + return [] + if g == 0: + return [ [i] + start for i in range(a, b+1) ] + + newg = g - 1 + newTin = Tin.submatrix(nrows=newg+1, ncols=newg+1) + pts = [] + for i in range(a, b+1): + v1 = vector(shift[:newg+1]) + v2 = vector(T.column(g)[:newg+1]) + newshift = v1 - (i - shift[g]) * newTin * v2 + newR = (R**2 - T[g, g]**2*(i - shift[g])**2).sqrt() + newstart = [i] + start + pts += find_integer_points(newg, newshift, newR, newstart) + if len(pts) > self._max_points_in_one_evaluation: + raise ValueError, "Evaluating theta function required too many summands (at least %s)" % len(pts) + + return pts + + return find_integer_points(self._g-1, leftshift, R, []) + + def radius(self, derivs): + r""" + Calculate the radius `R` to compute the value of a theta + function to within `2^{-P + 1}`. + + `P` is the precision in bits desired, computed from the input + matrix. + + `R` is the radius of [CRTF] Theorems 2, 4, and 6. + + EXAMPLES:: + + sage: from sage.functions.riemann_theta import RiemannTheta + sage: Omega1 = matrix(CDF, 2, 2, [1 + I, 2, 2, I]) + sage: E1 = RiemannTheta(Omega1); R = E1.radius([]); R + 8.10732947638 + sage: len(E1.integer_points(Omega1.column_space()(0), R)) + 69 + + sage: Omega2 = matrix(ComplexField(100), 2, 2, [1 + I, 2, 2, I]) + sage: E2 = RiemannTheta(Omega2); R = E2.radius([]); R + 10.8461012965 + sage: len(E2.integer_points(Omega2.column_space()(0), R)) + 121 + + sage: Omega = matrix(ComplexField(100), 2, 2, [1 + I, 2, 2, I]) + sage: E = RiemannTheta(1/10 * Omega); R = E.radius([]); R + 10.3551680946 + sage: len(E.integer_points(Omega.column_space()(0), R)) + 1085 + + sage: Omega = matrix(ComplexField(300), 2, 2, [1 + I, 2, 2, I]) + sage: E = RiemannTheta(1/10 * Omega); R = E.radius([]); R + 17.6483193006 + sage: len(E.integer_points(Omega.column_space()(0), R)) + 3125 + """ + Pi = self._ring.pi() + I = self._ring.gen() + g = self._g + + num_derivs = len(derivs) + norm_derivs = [ deriv.norm() for deriv in derivs ] + assert all( norm_deriv > 1e-10 for norm_deriv in norm_derivs ), "Directions need to be non-zero vectors" + prod_norm_derivs = prod( norm_derivs ) + + U = self._T._pari_().qflll()._sage_() # find shortest vector in lattice + v = (U * self._T).column(0) + r = v.norm() + verbose("Smallest vector is %s of length %s" % (v, r), level=2) + + d = ZZ(1) / min( self._T[i, i] for i in range(0, g) ) + detT = self._T.det() + + def compute_R(RRR, range_factor): + # we need to compute log(gamma_inc) for values whose log are + # about -self._ring.prec(), so we work in a larger precision ring temporarily + # attempt to compute to within 2^{-prec + 1} + const = -(self._ring.prec() - 1) + (RRR(g/2).gamma() * detT).log() - (RRR.pi()**(g/2) * prod_norm_derivs).log() + const = const.real() + def func(R): + # find_root requires RDF's, so we return an RDF giving the + # difference between the log of the desired function to + # minimize and the desired epsilon + v = sum( RRR(5)**(num_derivs - i) * + RRR(2)**(i*(num_derivs - 3/2)) * RRR(d)**i * + upper_gamma( RRR((g + i)/2), (R - r/2)**2) + for i in range(0, num_derivs+1) ) + v = v.real().log() + return RDF(v - const) + + left = (RRR(g + 2*num_derivs + (RRR(g)**2 + 8*num_derivs).sqrt()).sqrt() + r)/2 + left = RDF(left.real()) + + from sage.numerical.optimize import find_root + try: + rt = find_root(func, left, range_factor*left) + except RuntimeError: + rt = -left + + R = max(left, rt) + if func(R) > 2**(-20): + # Could not find an accurate bound for the radius R + return None + R += ZZ(1)/100 # plus a slight epsilon, to deal with errors in root finding to low precision + return RDF(R) + + verbose("Trying to compute R with low precision..", level=1) + RRR1 = CDF + range_factor1 = 10 + R = compute_R(RRR1, range_factor1) + if R is not None: + verbose("Trying to compute R with low precision.. SUCCEEDED", level=1) + else: + verbose("Trying to compute R with low precision.. FAILED", level=1) + verbose("Trying to compute R with high precision..", level=1) + RRR2 = ComplexField((ZZ(3)/2 * self._ring.prec()).ceil()) + range_factor2 = ZZ(self._ring.prec()).isqrt() + R = compute_R(RRR2, range_factor2) + if R is not None: + verbose("Trying to compute R with high precision.. SUCCEEDED", level=1) + else: + verbose("Trying to compute R with high precision.. FAILED", level=1) + raise ValueError, "Could not find an accurate bound for the radius R" + + R = RDF(R) + verbose("Computed radius R = %s" % R, level=1) + return R + + def exp_and_osc_at_point(self, z, derivs=[]): + r""" + Calculate the exponential and oscillating parts of `\theta(z, + \Omega)`. + + OUTPUT: + + A pair `u, v` such that `\theta(z, \Omega) = e^u \cdot v`. + + EXAMPLES:: + + sage: from sage.functions.riemann_theta import RiemannTheta + sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I]) + sage: E = RiemannTheta(1/10 * Omega) + sage: E.exp_and_osc_at_point(Omega.column_space()(0)) + (0.0, 4.41764896374 + 0.409494157593*I) + sage: E.value_at_point(Omega.column_space()(0)) + 4.41764896374 + 0.409494157593*I + sage: mathematica.SiegelTheta(1/10*Omega, [0, 0]).N(40) # optional - mathematica + 4.417648963714461 + 0.40949415758330304*I + + sage: E.exp_and_osc_at_point(Omega.column_space()([1/2, 2+2/3*I])) + (13.962634016, 3.78266700319 + 0.201093599178*I) + sage: E.value_at_point(Omega.column_space()([1/2, 2+2/3*I])) + 4382208.29921 + 232966.327327*I + sage: v = mathematica.SiegelTheta(1/10*Omega, [1/2, 2+2/3*I]).N(40); v.Re()/10^6, v.Im() # optional - mathematica + (4.382208299206245, 232966.32732660917) + + To higher precision:: + + sage: v = RiemannTheta(1/10*matrix(ComplexField(200), 2, 2, [1 + I, 2, 2, I])).value_at_point(vector([1/2, 2+2/3*I])) + sage: v.real()/10^6, v.imag() + (4.3822082992062459038743076349829658635958787958520566763119, + 232966.32732660939410871639806475516727466443779468195756093) + sage: v = mathematica.SiegelTheta(1/10*matrix(2, 2, [1 + I, 2, 2, I]), [1/2, 2+2/3*I]).N(100); v.Re()/10^6, v.Im() # optional - mathematica + (4.3822082992062459038743076349829658635958787958520566763119247589002543494309123447131812897634579125, + 232966.3273266093941087163980647551672746644377946819575609732661864324204353400769930567694437801562) + """ + Pi = self._ring.pi() + I = self._ring.gen() + TPII = 2*Pi*I + CS = self._Omega.column_space() + derivs = [ TPII * CS(deriv) for deriv in derivs ] # XXX + z = TPII * CS(z) + + x = z.apply_map(lambda t: t.real()) + y = z.apply_map(lambda t: t.imag()) + shift = self._Xin * x + intshift = shift.apply_map(lambda t: ZZ(t.round())) + + if not derivs: + if self._no_derivs_radius is None: + self._no_derivs_radius = self.radius(derivs) + R = self._no_derivs_radius + else: + R = self.radius(derivs) + verbose("Radius of S_R is %s" % R, level=2) + S = self.integer_points(z, R) + verbose("S_R consists of %s points" % len(S), level=2) + verbose("S_R = %s" % S, level=3) + + u, v = finite_sum(shift, intshift, x, y, self._X, self._Y, S, derivs, x.base_ring()) + osc_part = u + I*v + + exp_part = ZZ(1)/2 * (x * shift) + verbose("exponential part, oscillating part = (%s, %s)" % (exp_part, osc_part), level=2) + return exp_part, osc_part + + def value_at_point(self, z, derivs=[]): + r""" + EXAMPLES: + sage: from sage.functions.riemann_theta import RiemannTheta + sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I]) + sage: E = RiemannTheta(1/10 * Omega) + sage: E.value_at_point(0) + 4.41764896374 + 0.409494157593*I + sage: E.value_at_point([1/3, 2/3 + 1/2*I]) + 1582.31272298 + 996.188182401*I + """ + exp_part, osc_part = self.exp_and_osc_at_point(z, derivs=derivs) + scale = 1 + if derivs: + scale = z.base_ring().gen()**(len(derivs)-1) # experimentally derived fudge factor + return scale * exp_part.exp() * osc_part + + def derivatives_at_point(self, z, list_of_derivs=[]): + r""" + XXX + sage: from sage.functions.riemann_theta import RiemannTheta + sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I]) + sage: E = RiemannTheta(1/10 * Omega) + sage: u = vector(CDF, [1/3, 2/3 + 1/2*I]) + sage: X, X0, X1, X00, X01, X11 = E.derivatives_at_point(u, [ [(1, 0)], [(0, 1)], [(1, 0), (1, 0)], [(1, 0), (0, 1)], [(0, 1), (0, 1)] ]) + sage: X + 1582.31272298 + 996.188182401*I + sage: E.value_at_point(u) + 1582.31272298 + 996.188182401*I + sage: vector([X0, X1]) + (6869.0440534 - 11042.144279*I, 61131.0804751 - 49270.0715471*I) + sage: E.gradient_at_point(u) + (6869.0440534 - 11042.144279*I, 61131.0804751 - 49270.0715471*I) + sage: matrix(2, 2, [X00, X01, X01, X11]) + [ 34666.7793688 - 139103.389251*I -481117.989651 + 81904.8136253*I] + [-481117.989651 + 81904.8136253*I -1460873.9145 - 2861525.00205*I] + sage: E.hessian_at_point(u) + [ 34666.7793688 - 139103.389251*I -481117.989651 + 81904.8136253*I] + [-481117.989651 + 81904.8136253*I -1460873.9145 - 2861525.00205*I] + + sage: X, X1, X00, X0, X11, X01 = E.derivatives_at_point(u, [ [(0, 1)], [(1, 0), (1, 0)], [(1, 0)], [(0, 1), (0, 1)], [(1, 0), (0, 1)] ]) + sage: X + 1582.31272298 + 996.188182401*I + sage: vector([X0, X1]) + (6869.0440534 - 11042.144279*I, 61131.0804751 - 49270.0715471*I) + sage: matrix(2, 2, [X00, X01, X01, X11]) + [ 34666.7793688 - 139103.389251*I -481117.989651 + 81904.8136253*I] + [-481117.989651 + 81904.8136253*I -1460873.9145 - 2861525.00205*I] + + sage: (X, ) = E.derivatives_at_point(u, []) + sage: X + 1582.31272298 + 996.188182401*I + """ + Pi = self._ring.pi() + I = self._ring.gen() + TPII = 2*Pi*I + + # scale derivatives appropriately + CS = self._Omega.column_space() + for i in range(len(list_of_derivs)): + derivs = list_of_derivs[i] + derivs = [ TPII * CS(deriv) for deriv in derivs ] + list_of_derivs[i] = derivs + # scale input point appropriately + z = TPII * self._Omega.column_space()(z) + + x = z.apply_map(lambda t: t.real()) + y = z.apply_map(lambda t: t.imag()) + shift = self._Xin * x + intshift = shift.apply_map(lambda t: ZZ(t.round())) + + # the radius should be the largest product of norms *with the most terms in the product* + if len(list_of_derivs) > 0: + temps = [] + for derivs in list_of_derivs: + norm_derivs = prod( deriv.norm() for deriv in derivs ) + temps.append((len(derivs), norm_derivs)) + max_terms, max_norm = max(temps) + for temp in temps: + # need both conditions + assert max_terms >= temp[0] + assert max_norm >= temp[1] + + # now actually find the largest + max_norm = 0 + max_elem = None + for derivs in list_of_derivs: + if len(derivs) < max_terms: + continue + norm_derivs = prod( deriv.norm() for deriv in derivs ) + if norm_derivs > max_norm: + max_norm = norm_derivs + max_elem = derivs + assert max_elem is not None + else: + max_elem = [] + + R = self.radius(max_elem) + + verbose("Radius of S_R is %s" % R, level=2) + S = self.integer_points(z, R) + verbose("S_R consists of %s points" % len(S), level=2) + verbose("S_R = %s" % S, level=3) + + uvs = finite_sum_with_many_derivatives(shift, intshift, x, y, self._X, self._Y, S, list_of_derivs, x.base_ring()) + vals = [] + for i in range(len(uvs)): # the first item is the value of the function at the point + u, v = uvs[i] + + scale = 1 + if i > 0: + derivs = list_of_derivs[i-1] # again, first item is value of function + scale = z.base_ring().gen()**(len(derivs)-1) # experimentally derived fudge factor + osc_part = u + I*v + exp_part = ZZ(1)/2 * (x * shift) + verbose("scale, exponential part, oscillating part = (%s, %s, %s)" % (scale, exp_part, osc_part), level=2) + vals.append(scale * exp_part.exp() * osc_part) + return vals + + def gradient_at_point(self, z): + r""" + EXAMPLES: + sage: from sage.functions.riemann_theta import RiemannTheta + sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I]) + sage: E = RiemannTheta(1/10 * Omega) + sage: E.gradient_at_point([1/3, 2/3 + 1/2*I]) + (6869.0440534 - 11042.144279*I, 61131.0804751 - 49270.0715471*I) + sage: E.gradient_at_point([1/3*I, 2/3*I + 1/2]) + (706078145.06 - 649539106.278*I, 2333263039.98 - 412406488.368*I) + """ + Pi = self._ring.pi() + I = self._ring.gen() + TPII = 2*Pi*I + v = [0]*self.genus() + list_of_derivs = [] + for i in range(0, self.genus()): + v[i] = TPII + list_of_derivs.append([ vector(self._ring, v) ]) + v[i] = 0 + z = TPII * self._Omega.column_space()(z) + + x = z.apply_map(lambda t: t.real()) + y = z.apply_map(lambda t: t.imag()) + shift = self._Xin * x + intshift = shift.apply_map(lambda t: ZZ(t.round())) + + # R = max( self.radius(derivs) for derivs in list_of_derivs ) # XXX + if self._gradient_radius is None: + self._gradient_radius = self.radius(list_of_derivs[-1]) # all derivatives have equal norm + R = self._gradient_radius + verbose("Radius of S_R is %s" % R, level=2) + S = self.integer_points(z, R) + verbose("S_R consists of %s points" % len(S), level=2) + verbose("S_R = %s" % S, level=3) + + uvs = finite_sum_with_many_derivatives(shift, intshift, x, y, self._X, self._Y, S, list_of_derivs, x.base_ring()) + vals = [] + for u, v in uvs[1:]: # the first value is the value of the function at the point + osc_part = u + I*v + exp_part = ZZ(1)/2 * (x * shift) + verbose("exponential part, oscillating part = (%s, %s)" % (exp_part, osc_part), level=2) + vals.append(exp_part.exp() * osc_part) + return vector(self._ring, vals) + + def hessian_at_point(self, z): + r""" + EXAMPLES: + sage: from sage.functions.riemann_theta import RiemannTheta + sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I]) + sage: E = RiemannTheta(1/10 * Omega) + sage: E.hessian_at_point([1/3, 2/3 + 1/2*I]) + [ 34666.7793688 - 139103.389251*I -481117.989651 + 81904.8136253*I] + [-481117.989651 + 81904.8136253*I -1460873.9145 - 2861525.00205*I] + sage: E.hessian_at_point([1/3*I, 2/3*I + 1/2]) + [-16895727870.4 - 11148594973.2*I -21658283152.2 - 41419586987.9*I] + [-21658283152.2 - 41419586987.9*I -3006810401.8 - 130403155418.0*I] + """ + Pi = self._ring.pi() + I = self._ring.gen() + TPII = 2*Pi*I + v1 = [0]*self.genus() + v2 = [0]*self.genus() + list_of_derivs = [] + for i in range(0, self.genus()): + v1[i] = TPII + for j in range(i, self.genus()): + v2[j] = TPII + list_of_derivs.append([ vector(self._ring, v1), vector(self._ring, v2) ]) + v2[j] = 0 + v1[i] = 0 + + z = TPII * self._Omega.column_space()(z) + x = z.apply_map(lambda t: t.real()) + y = z.apply_map(lambda t: t.imag()) + shift = self._Xin * x + intshift = shift.apply_map(lambda t: ZZ(t.round())) + + # R = max( self.radius(derivs) for derivs in list_of_derivs ) # XXX + if self._hessian_radius is None: + self._hessian_radius = self.radius(list_of_derivs[-1]) # all derivatives have equal norm + R = self._hessian_radius + verbose("Radius of S_R is %s" % R, level=2) + S = self.integer_points(z, R) + verbose("S_R consists of %s points" % len(S), level=2) + verbose("S_R = %s" % S, level=3) + + uvs = finite_sum_with_many_derivatives(shift, intshift, x, y, self._X, self._Y, S, list_of_derivs, x.base_ring()) + vals = [] + scale = z.base_ring().gen() # experimentally defined fudge factor + for u, v in uvs[1:]: # the first value is the value of the function at the point + osc_part = u + I*v + exp_part = ZZ(1)/2 * (x * shift) + verbose("exponential part, oscillating part = (%s, %s)" % (exp_part, osc_part), level=2) + vals.append(scale * exp_part.exp() * osc_part) + + vals.reverse() + + m = matrix(self._ring, self.genus(), self.genus()) + for i in range(0, self.genus()): + for j in range(i, self.genus()): + m[i, j] = vals.pop() + for i in range(1, self.genus()): + for j in range(0, i): + m[i, j] = m[j, i] + return m + + def siegel_theta(self, u): + return self.value_at_point(u) + + def siegel_theta_with_char(self, characteristics, u): + r, s = characteristics + CS = self._Omega.column_space() + r = CS(r) + s = CS(s) + u = CS(u) + + R = self._ring + II = R.gen() + PI = R.pi() + z = (self._Omega/(-2*PI*II)) + scale = (2*PI*II * ( ZZ(1)/2*r*z*r + r*(u + s) )).exp() + u_ = u + z*r + s + return scale * self.siegel_theta(u_) + + def shimura_phi_with_char(self, characteristics, u): + R = self._ring + II = R.gen() + PI = R.pi() + z = (self._Omega/(-2*PI*II)) + u = z.column_space()(u) + Z = z - z.apply_map(lambda t: t.conjugate()) + Z = invert(Z) + scale = (II * PI * (u * Z * u)).exp() + return scale * self.siegel_theta_with_char(characteristics, u) + + def shimura_phi(self, u): + CS = self._Omega.column_space() + characteristics = (CS(0), CS(0)) + return self.shimura_phi_with_char(characteristics, u) + + def theta_dict(self, u, list_of_characteristics=None): + r""" + sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6]) + sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2) + sage: D = RT.theta_dict(0, [ ((0, 0), (0, 0)), ((0, 1/2), (0, 1/2)) ]) + sage: D[((0, 0), (0, 0))] + 1.00171421242 + sage: D[((0, 1/2), (0, 1/2))].abs() < 1e-15 + True + + sage: RT.theta_dict([1/2, 1/5 + CDF(I)*0.234], [ ((0, 1/2), (0, 1/2)), ((0, 0), (0, 0)) ]) + {((0, 0), (0, 0)): 0.998788360598 + 0.00120506039896*I, ((0, 1/2), (0, 1/2)): -0.136735147549 - 0.117846448587*I} + + sage: RT.theta_vector([1/2, 1/5 + CDF(I)*0.234], [ ((0, 0), (0, 0)), ((0, 1/2), (0, 1/2)) ]) + (0.998788360598 + 0.00120506039896*I, -0.136735147549 - 0.117846448587*I) + sage: RT.theta_vector([1/2, 1/5 + CDF(I)*0.234], [ ((0, 1/2), (0, 1/2)), ((0, 0), (0, 0)) ]) + (-0.136735147549 - 0.117846448587*I, 0.998788360598 + 0.00120506039896*I) + """ + if list_of_characteristics is None: + list_of_characteristics = list(self._all_characteristics()) + L = list(_global_theta_for_parallel([ (c, self, u) for c in list_of_characteristics ])) + D = {} + for (args, kwds), val in L: + D[args[0]] = val + return D + + def even_theta_dict(self, u): + return self.theta_dict(u, self._even_characteristics()) + + def odd_theta_dict(self, u): + return self.theta_dict(u, self._odd_characteristics()) + + def even_theta_vector(self, u): + return self.theta_vector(u, self._even_characteristics()) + + def odd_theta_vector(self, u): + return self.theta_vector(u, self._odd_characteristics()) + + def theta_vector(self, u, list_of_characteristics=None): + r""" + sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6]) + sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2) + sage: RT.theta_vector(0) + (1.00171421242, 0.998608586772 + 2.06043088621e-36*I, 0.999838600402 + 6.56441798402e-36*I, 0.999838600402 + 1.11674199318e-36*I, 0.283260714578, 1.73435525065e-17 - 6.52327288837e-17*I, 0.278207060819 + 2.68012527741e-36*I, 1.7344656162e-17 - 4.25398925609e-17*I, 0.197753422929, 0.181367479122 - 3.91422304159e-34*I, 1.21087425915e-17 + 8.00949937526e-17*I, 1.21087425915e-17 + 8.00979583156e-17*I, 0.283260714578, 1.73435525065e-17 + 5.02857477884e-16*I, 1.7344656162e-17 - 5.08846937396e-16*I, 0.278207060819 - 4.93038065763e-32*I) + """ + if list_of_characteristics is None: + list_of_characteristics = list(self._all_characteristics()) + D = self.theta_dict(u, list_of_characteristics) + + L = [] + for c in list_of_characteristics: + L.append(D[c]) + return vector(L) + + def _all_characteristics(self): + r""" + sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6]) + sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2) + sage: list(RT._all_characteristics()) + [((0, 0), (0, 0)), ((0, 0), (1/2, 0)), ((0, 0), (0, 1/2)), ((0, 0), (1/2, 1/2)), + ((1/2, 0), (0, 0)), ((1/2, 0), (1/2, 0)), ((1/2, 0), (0, 1/2)), ((1/2, 0), (1/2, 1/2)), + ((0, 1/2), (0, 0)), ((0, 1/2), (1/2, 0)), ((0, 1/2), (0, 1/2)), ((0, 1/2), (1/2, 1/2)), + ((1/2, 1/2), (0, 0)), ((1/2, 1/2), (1/2, 0)), ((1/2, 1/2), (0, 1/2)), ((1/2, 1/2), (1/2, 1/2))] + """ + for u in GF(2)**(self._g): + for v in GF(2)**(self._g): + yield (tuple(u.change_ring(ZZ)/2), tuple(v.change_ring(ZZ)/2)) + + def _even_characteristics(self): + r""" + sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6]) + sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2) + sage: list(RT._even_characteristics()) + [((0, 0), (0, 0)), ((0, 0), (1/2, 0)), ((0, 0), (0, 1/2)), ((0, 0), (1/2, 1/2)), + ((1/2, 0), (0, 0)), ((1/2, 0), (0, 1/2)), ((0, 1/2), (0, 0)), ((0, 1/2), (1/2, 0)), + ((1/2, 1/2), (0, 0)), ((1/2, 1/2), (1/2, 1/2))] + """ + for u in GF(2)**(self._g): + for v in GF(2)**(self._g): + u = u.change_ring(ZZ) + v = v.change_ring(ZZ) + if u*v % 2 == 0: + yield (tuple(u/2), tuple(v/2)) + + def _odd_characteristics(self): + r""" + sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6]) + sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2) + sage: list(RT._odd_characteristics()) + [((1/2, 0), (1/2, 0)), ((1/2, 0), (1/2, 1/2)), ((0, 1/2), (0, 1/2)), + ((0, 1/2), (1/2, 1/2)), ((1/2, 1/2), (1/2, 0)), ((1/2, 1/2), (0, 1/2))] + """ + for u in GF(2)**(self._g): + for v in GF(2)**(self._g): + u = u.change_ring(ZZ) + v = v.change_ring(ZZ) + if u*v % 2 == 1: + yield (tuple(u/2), tuple(v/2)) + +def elliptic_theta(n, z, q): + r""" + Compute Jacobi elliptic theta functions. + + The Jacobi elliptic theta functions `\vartheta_n : CC \times CC + \to CC` are defined as + + - `\vartheta_1(z, q) = \sum_{n = -\infty}^{\infty} (-1)^{n - 1/2} q^{(n + 1/2)^2} e^{(2 n + 1) i z}` + + - `\vartheta_2(z, q) = \sum_{n = -\infty}^{\infty} q^{(n + 1/2)^2} e^{(2 n + 1) i z}` + + - `\vartheta_3(z, q) = \sum_{n = -\infty}^{\infty} q^{n^2} e^{2 n i z}` + + - `\vartheta_4(z, q) = \sum_{n = -\infty}^{\infty} (-1)^n q^{n^2} e^{2 n i z}` + + INPUT: + + - `n` denotes the number of the elliptic theta function, 1, 2, 3, or 4 + + - `z` denotes the argument in the complex plane + + - `q` denotes the *nome*, `q = e^{i \pi \tau}` with `\tau` in the complex upper half plane. + + OUTPUT: + + An element `\vartheta_n(z, q)` in the same complex field as the nome `q`. + + .. note: + + This notation and definition has been chosen to agree with + Mathematica's ``EllipticTheta[n, z, q]`` function. + + EXAMPLES: + + Let's check that our values of the elliptic theta functions + agree with Mathematica's:: + + sage: q = 2*CDF(e^(-pi)) + sage: elliptic_theta(1, 0, q).abs() < 1e-10 + True + sage: elliptic_theta(2, 0, q).real() + 1.09251047534 + sage: elliptic_theta(2, 0, q).imag().abs() < 1e-10 + True + sage: elliptic_theta(3, 0, q) + 1.17296726855 + sage: elliptic_theta(4, 0, q).real() + 0.827255921362 + sage: elliptic_theta(4, 0, q).imag().abs() < 1e-10 + True + + sage: mathematica.EllipticTheta(1, 0, 2*e^-pi).N(16) # optional - mathematica + 0 + sage: mathematica.EllipticTheta(2, 0, 2*e^-pi).N(16) # optional - mathematica + 1.0925104753387513 + sage: mathematica.EllipticTheta(3, 0, 2*e^-pi).N(16) # optional - mathematica + 1.1729672685486494 + sage: mathematica.EllipticTheta(4, 0, 2*e^-pi).N(16) # optional - mathematica + 0.82725592136214803 + + At a non-zero point:: + + sage: elliptic_theta(1, 1/3, q).real(), elliptic_theta(1, 1/3, q).imag().abs() < 1e-10 + (0.34799740214, True) + sage: elliptic_theta(2, 1/3, q).real(), elliptic_theta(2, 1/3, q).imag().abs() < 1e-10 + (1.02909707137, True) + sage: elliptic_theta(3, 1/3, q).real(), elliptic_theta(3, 1/3, q).imag().abs() < 1e-10 + (1.13587132251, True) + sage: elliptic_theta(4, 1/3, q).real(), elliptic_theta(4, 1/3, q).imag().abs() < 1e-10 + (0.864181180143, True) + + sage: mathematica.EllipticTheta(1, 1/3, 2*e^-pi).N(16) # optional - mathematica + 0.34799740214022461 + sage: mathematica.EllipticTheta(2, 1/3, 2*e^-pi).N(16) # optional - mathematica + 1.0290970713692956 + sage: mathematica.EllipticTheta(3, 1/3, 2*e^-pi).N(16) # optional - mathematica + 1.1358713225095661 + sage: mathematica.EllipticTheta(4, 1/3, 2*e^-pi).N(16) # optional - mathematica + 0.86418118014343551 + + Some random values, to different precisions:: + + sage: elliptic_theta(4, 1+I, CDF(e^(-pi))) + 1.13518915646 + 0.285173964442*I + sage: mathematica.EllipticTheta(4, 1+I, e^(-pi)).N(16) # optional - mathematica + 1.13518915646320072 + 0.28517396444192508*I + + sage: elliptic_theta(3, 0, -1/4*CDF(e^(-1/10*pi + 1/5))).real(), elliptic_theta(3, 0, -1/4*CDF(e^(-1/10*pi + 1/5))).imag().abs() < 1e-16 + (0.55888785571, True) + sage: mathematica.EllipticTheta(3, 0, -1/4*e^(-pi/10 + 1/5)).N(53) # optional - mathematica + 0.55888785570955688429958704552908258979432384584578182 + + sage: elliptic_theta(3, 0, 1/3*ComplexField(100)(e^(-1/10*pi))) + 1.4939685355395330059989284819 + sage: mathematica.EllipticTheta(3, 0, 1/3*e^(-pi/10)).N(53) # optional - mathematica + 1.49396853553953300599892848191670064396107631596893401 + + EXAMPLES: + + It is a theorem that + + .. math:: + + \vartheta_3( 0, e^{-\pi} ) = \frac{\pi^{1/4}}{\Gamma(3/4)} . + + :: + + sage: elliptic_theta(3, 0, CDF(e^(-pi))) + 1.08643481121 + sage: CDF(pi)^(1/4) / gamma(CDF(3/4)) + 1.08643481121 + + sage: elliptic_theta(3, 0, ComplexField(100)(e^(-pi))) + 1.0864348112133080145753161215 + sage: ComplexField(100)(pi)^(1/4) / gamma(ComplexField(100)(3/4)) + 1.0864348112133080145753161215 + + sage: elliptic_theta(3, 0, ComplexField(200)(e^(-pi))) + 1.0864348112133080145753161215102234570702057072452188859208 + sage: ComplexField(200)(pi)^(1/4) / gamma(ComplexField(200)(3/4)) + 1.0864348112133080145753161215102234570702057072452188859208 + + Which we can verify in Mathematica as well:: + + sage: ((mathematica(pi)^(1/4)) / mathematica.Gamma(3/4)).N(100) # optional - mathematica + 1.086434811213308014575316121510223457070205707245218885920790315981856732267109795960561618489679764 + sage: mathematica.EllipticTheta(3, 0, e^-pi).N(100) # optional - mathematica + 1.086434811213308014575316121510223457070205707245218885920790315981856732267109795960561618489679764 + + TESTS: + + The elliptic theta functions are quasi-doubly periodic and + satisfy various functional equations:: + + sage: II = ComplexField(100)(I) + sage: PI = ComplexField(100).pi() + sage: tau = II + sage: q = (tau * PI * II).exp() + sage: z = 2/5*tau + 1/3 + + One quasi-period is `\pi`:: + + sage: v = elliptic_theta(1, z + PI, q)/elliptic_theta(1, z, q); v.real(), v.imag().abs() < 1e-20 + (-1.0000000000000000000000000000, True) + sage: v = elliptic_theta(2, z + PI, q)/elliptic_theta(2, z, q); v.real(), v.imag().abs() < 1e-20 + (-1.0000000000000000000000000000, True) + sage: v = elliptic_theta(3, z + PI, q)/elliptic_theta(3, z, q); v.real(), v.imag().abs() < 1e-20 + (1.0000000000000000000000000000, True) + sage: v = elliptic_theta(4, z + PI, q)/elliptic_theta(4, z, q); v.real(), v.imag().abs() < 1e-20 + (1.0000000000000000000000000000, True) + + The second quasi-period is `i \pi` (in general, `\tau \pi`):: + + sage: elliptic_theta(1, z + II*PI, q)/elliptic_theta(1, z, q) + -40.47363290142101126706943987... + 31.84639025962267153445618616...*I + sage: elliptic_theta(2, z + II*PI, q)/elliptic_theta(2, z, q) + 40.47363290142101126706943987... - 31.84639025962267153445618616...*I + sage: elliptic_theta(3, z + II*PI, q)/elliptic_theta(3, z, q) + 40.47363290142101126706943987... - 31.84639025962267153445618616...*I + sage: elliptic_theta(4, z + II*PI, q)/elliptic_theta(4, z, q) + -40.47363290142101126706943987... + 31.84639025962267153445618616...*I + + The factor of automorphy in the second case is `q^{-1} e^{-2 \pi z}`:: + + sage: ~q * (-2*II*z).exp() + 40.47363290142101126706943987... - 31.84639025962267153445618616...*I + """ + z = q.parent()(z) + I = q.parent().gen() + Pi = q.parent().pi() + + if n == 1: + tau = q.log()/(I*Pi) + return -I * (I*z + Pi*I*tau/4).exp() * elliptic_theta(4, z + Pi*ZZ(1)/2*tau, q) + if n == 2: + return elliptic_theta(1, z + Pi/2, q) + if n == 3: + tau = q.log()/(I*Pi) + Omega = matrix(1, 1, [tau]) + z = ~Pi * Omega.column_space()([z]) + return RiemannTheta(Omega).value_at_point(z) + if n == 4: + return elliptic_theta(3, z - Pi/2, q) + raise ValueError, "n must be 1, 2, 3, or 4 in elliptic_theta" + +def elliptic_theta_prime(n, z, q): + r""" + Compute derivatives of Jacobi elliptic theta functions. + + The derivatives of the Jacobi elliptic theta functions + `\vartheta_n' : CC \times CC \to CC`. + + INPUT: + + - `n` denotes the number of the elliptic theta function, 1, 2, 3, or 4 + + - `z` denotes the argument in the complex plane + + - `q` denotes the *nome*, `q = e^{i \pi \tau}` with `\tau` in the + complex upper half plane. + + OUTPUT: + + An element `\vartheta_n'(z, q)` in the same complex field as the nome `q`. + + .. note: + + This notation and definition has been chosen to agree with + Mathematica's ``EllipticThetaPrime[n, z, q]`` function. + + EXAMPLES: + + Derivatives of different derivatives of theta functions:: + + sage: v = elliptic_theta_prime(1, I, CDF(e^-pi)); v.real(), v.imag().abs() < 1e-12 + (1.35566883407, True) + sage: v = elliptic_theta_prime(2, I, CDF(e^-pi)); v.imag(), v.real().abs() < 1e-12 + (-1.1228178842, True) + sage: v = elliptic_theta_prime(3, I, CDF(e^-pi)); v.imag(), v.real().abs() < 1e-12 + (-0.62768475242, True) + sage: v = elliptic_theta_prime(4, I, CDF(e^-pi)); v.imag(), v.real().abs() < 1e-12 + (0.626162043874, True) + + sage: mathematica.EllipticThetaPrime(1, I, e^(-pi)).N(16) # optional - mathematica + 1.3556688340695497 + sage: mathematica.EllipticThetaPrime(2, I, e^(-pi)).N(16) # optional - mathematica + -1.1228178842005472*I + sage: mathematica.EllipticThetaPrime(3, I, e^(-pi)).N(16) # optional - mathematica + -0.62768475242048347*I + sage: mathematica.EllipticThetaPrime(4, I, e^(-pi)).N(16) # optional - mathematica + 0.62616204387425979*I + + Derivates at different points:: + + sage: v = elliptic_theta_prime(3, 1, CDF(e^-pi)); v.real(), v.imag().abs() < 1e-12 + (-0.157156104884, True) + sage: v = elliptic_theta_prime(3, 2, CDF(e^-pi)); v.real(), v.imag().abs() < 1e-12 + (0.130790002852, True) + sage: v = elliptic_theta_prime(3, 3, CDF(e^-pi)); v.real(), v.imag().abs() < 1e-12 + (0.0483135237156, True) + sage: v = elliptic_theta_prime(3, 4, CDF(e^-pi)); v.real(), v.imag().abs() < 1e-12 + (-0.171008153468, True) + + sage: mathematica.EllipticThetaPrime(3, 1, e^(-pi)).N(16) # optional - mathematica + -0.15715610488427422 + sage: mathematica.EllipticThetaPrime(3, 2, e^(-pi)).N(16) # optional - mathematica + 0.1307900028522554 + sage: mathematica.EllipticThetaPrime(3, 3, e^(-pi)).N(16) # optional - mathematica + 0.048313523715647603 + sage: mathematica.EllipticThetaPrime(3, 4, e^(-pi)).N(16) # optional - mathematica + -0.17100815346753069 + + And at a different modulus, to higher precision:: + + sage: elliptic_theta_prime(3, 1/5 + 2*I, ComplexField(200)(e^(-pi + I))) + 2.6604481224103801421791420912309717423496926664711512178563 - + 3.8528130836653595482465286825791971210966140276918135528629*I + sage: mathematica.EllipticThetaPrime(3, 1/5 + 2*I, e^(-pi + I)).N(60) # optional - mathematica + 2.660448122410380142179142091230971742349692666471151217856297 - + 3.852813083665359548246528682579197121096614027691813552862919*I + + sage: elliptic_theta_prime(3, -9/5 + I, ComplexField(200)(1/5*e^(-2*pi + I))) + -0.0054325201726101704353440405791768933668088597580411363599553 + + 0.00053239599348395256664970257021495457788140582724507952668306*I + sage: mathematica.EllipticThetaPrime(3, -9/5 + I, 1/5*e^(-2*pi + I)).N(60) # optional - mathematica + -0.00543252017261017043534404057917689336680885975804113635995534 + + 0.000532395993483952566649702570214954577881405827245079526683083*I + """ + I = q.parent().gen() + Pi = q.parent().pi() + + if n == 1: + tau = q.log()/(I*Pi) + Omega = matrix(1, 1, [tau]) + z = Omega.base_ring()(z) + zp = Omega.base_ring()(z + Pi*ZZ(1)/2*tau) + return -I * (I*z + Pi*I*tau/4).exp() * (I*elliptic_theta(4, zp, q) + elliptic_theta_prime(4, zp, q)) + # return -I * (I*z + Pi*I*tau/4).exp() * elliptic_theta_prime(4, z + Pi*ZZ(1)/2*tau, q) + if n == 2: + return elliptic_theta_prime(1, z + Pi/2, q) + if n == 3: + tau = q.log()/(I*Pi) + Omega = matrix(1, 1, [tau]) + z = ~Pi * Omega.column_space()([z]) + return RiemannTheta(Omega).value_at_point(z, derivs=[vector([~Pi])]) + if n == 4: + return elliptic_theta_prime(3, z - Pi/2, q) + raise ValueError, "n must be 1, 2, 3, or 4 in elliptic_theta" + +def siegel_theta(Omega, z): + r""" + Calculate the Riemann-Siegel theta function associated to `\Omega` + at the point `z`. + + Let `g` be the genus (dimension) of the theta funciton. The + Riemann-Siegel theta function `\theta : CC^g \times H_g \to CC` is + defined by the infinite series + + .. math:: + + \theta(z, \Omega) = \sum_{n \in ZZ^g} \exp \left( 2 \pi i \left( + n^{t} \cdot \Omega \cdot n + n^t \cdot z \right) \right) . + + .. note: + + This notation and definition has been chosen to agree with + Mathematica's ``SiegelTheta[Omega, z]`` function. + + EXAMPLES: + + Here's a genus 2 example:: + + sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6]) + sage: siegel_theta(tau2, vector([0, 0])) + 1.00171421242 + sage: siegel_theta(tau2, vector([1/3 + 1/2*CDF(I), 0])) + 0.991161256811 - 0.0155303295308*I + sage: mathematica.SiegelTheta(tau2, vector([0, 0])).Re() # optional - mathematica + 1.0017142124239051 + sage: mathematica.SiegelTheta(tau2, vector([1/3 + 1/2*I, 0])).N(16) # optional - mathematica + 0.9911612568105869 - 0.015530329530843857*I + + And here's a more complicated genus 3 example:: + + sage: im = diagonal_matrix([1, 2, 5]) + sage: X = matrix(QQ, 3, 3, [-2, 1/2, 1, 1, 0, -2, 1, -1, 0]) + sage: im = X.transpose() * im * X; im + [ 11 -6 -6] + [ -6 21/4 1/2] + [ -6 1/2 9] + sage: im.det() + 40 + sage: re = matrix(QQ, 3, 3, [0, 1, 3, 1, 2, -1/2, 3, -1/2, -2]) + sage: tau3 = re + ComplexField(100).gen()*im; tau3.change_ring(CDF) + [ 11.0*I 1.0 - 6.0*I 3.0 - 6.0*I] + [ 1.0 - 6.0*I 2.0 + 5.25*I -0.5 + 0.5*I] + [ 3.0 - 6.0*I -0.5 + 0.5*I -2.0 + 9.0*I] + + sage: v = siegel_theta(tau3, 0); v.real(), v.imag().abs() < 1e-30 + (0.99830741581318716996854547831, True) + sage: siegel_theta(tau3, 1/10*vector([1/2 + I, 2/3*I, 1/4*I])) + 0.99708919647052881852125605664 + 0.00077789981278179739800207580479*I + sage: mathematica.SiegelTheta(tau3, vector([0, 0, 0])).Re() # optional - mathematica + 0.99830741581318716996854547831117 + + sage: v = 1/10*vector(ComplexField(100), [1/2 + I, 2/3*I, 1/4*I]) + sage: mathematica.SiegelTheta(tau3, v) # optional - mathematica + 0.99708919647052881852125605664603 + 0.00077789981278179739800207580478*I + + Now with higher precision:: + + sage: tau3p = re + ComplexField(200).gen()*im + sage: v = siegel_theta(tau3p, 0); v.real(), v.imag().abs() < 1e-60 + (0.99830741581318716996854547831116617197514122934164610471589, True) + sage: siegel_theta(tau3p, 1/10*vector([1/2 + I, 2/3*I, 1/4*I])) + 0.99708919647052881852125605664602767210203416266645540488523 + 0.00077789981278179739800207580477768179379021250559152499833916*I + + sage: def siegel_theta_m(tau, u): + ... mtau = magma(tau) + ... mu = magma(matrix(tau.base_ring(), tau.nrows(), 1, list(u))) + ... mc = magma(matrix(tau.base_ring(), 2*tau.nrows(), 1)) + ... return magma.Theta(mc, mu, mtau) + + sage: siegel_theta_m(tau3p, 1/10*vector([1/2 + I, 2/3*I, 1/4*I])) # optional - magma + 0.994949815318003698123521324524667920491325676563809855469316 + 0.00153595088862631852918813432460329503001550616121771001338577*$.1 + """ + return RiemannTheta(Omega).siegel_theta(z) + +def siegel_theta_with_char(characteristics, z, u): + r""" + EXAMPLES:: + + sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6]) + + sage: siegel_theta(tau2, vector([0, 0])) + 1.00171421242 + sage: mathematica.SiegelTheta(tau2, vector([0, 0])).Re() # optional - mathematica + 1.0017142124239051 + + sage: siegel_theta_with_char([[1/2, 0], [0, 0]], tau2, vector([0, 0])) + 0.283260714578 + sage: mathematica.SiegelTheta([[1/2, 0], [0, 0]], tau2, vector([0, 0])) # optional - mathematica + 0.2832607145783672 + 0.*I + + sage: siegel_theta_with_char([[1/2, 0], [0, 0]], 1 + 2*tau2, vector([0, 0])) + 0.0278618215705 + 0.0278618215705*I + sage: mathematica.SiegelTheta([[1/2, 0], [0, 0]], 1 + 2*tau2, vector([0, 0])) # optional - mathematica + 0.02786182157051039 + 0.027861821570510384*I + + Odd theta functions are zero at the origin: + + sage: im = diagonal_matrix([1, 2, 5]) + sage: X = matrix(QQ, 3, 3, [-2, 1/2, 1, 1, 0, -2, 1, -1, 0]) + sage: im = X.transpose() * im * X; im + [ 11 -6 -6] + [ -6 21/4 1/2] + [ -6 1/2 9] + sage: im.det() + 40 + sage: re = matrix(QQ, 3, 3, [0, 1, 3, 1, 2, -1/2, 3, -1/2, -2]) + sage: tau3 = re + ComplexField(100).gen()*im; tau3.change_ring(CDF) + [ 11.0*I 1.0 - 6.0*I 3.0 - 6.0*I] + [ 1.0 - 6.0*I 2.0 + 5.25*I -0.5 + 0.5*I] + [ 3.0 - 6.0*I -0.5 + 0.5*I -2.0 + 9.0*I] + + sage: v = siegel_theta_with_char([ [1/2, 0, 0], [1/2, 0, 0] ], tau3, vector([0, 0, 0])); v.abs() < 1e-30 + True + sage: mathematica.SiegelTheta([ [1/2, 0, 0], [1/2, 0, 0] ], tau3, vector([0, 0, 0])).Abs() < 1e-30 # optional - mathematica + True + """ + return RiemannTheta(z).siegel_theta_with_char(characteristics, u) + +def shimura_phi_with_char(characteristics, z, u): + return RiemannTheta(z).shimura_phi_with_char(characteristics, u) + +def shimura_phi(z, u): + return RiemannTheta(z).shimura_phi(u) diff -r 341d08e69736 -r d6e4600f2353 sage/functions/riemann_theta_tests.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/riemann_theta_tests.py Thu Jun 04 22:47:59 2009 -0700 @@ -0,0 +1,702 @@ +r""" +Functional tests of Riemann theta functions. + +See ``sage.functions.riemann_theta_tests_{sage, maple, mathematica}`` +for docstrings that actually exercise these tests. They have been +split out so that parallel doctesting is faster. + +AUTHORS: + +- Nick Alexander (2009-03) +""" + +from sage.functions.all import * +from sage.functions.riemann_theta import invert, RiemannTheta + +from sage.rings.all import ZZ, QQ, ComplexField, CC, CDF, RealField, RR, RDF +from sage.modules.all import vector +from sage.matrix.all import matrix, identity_matrix, MatrixSpace, diagonal_matrix +from sage.interfaces.all import maple, mathematica, magma + +from sage.misc.sage_eval import sage_eval +from sage.misc.misc import verbose +from sage.misc.cachefunc import cached_function, cached_method +from sage.misc.prandom import choice + +def maple_to_CDF(m): + t = m.evalf(16) + v = sage_eval(repr(t), locals={'I':CDF.gen()}) + return v + +def msiegel_theta(Omega, z): + return mathematica.SiegelTheta(Omega, z).N(16) + +def melliptic_theta(n, z, q): + return mathematica.EllipticTheta(n, z, q).N(16) + +def melliptic_theta_prime(n, z, q): + return mathematica.EllipticThetaPrime(n, z, q).N(16) + +def rough_compare(m, s): + if s.abs() > 1e30: + # too hard to check + return True + # we check if the values or close, and failing that, if the + # first several digits agree. Removing the periods + # compensates for scientific notation + a = (m.Re() - s.real()).Abs() < 1e-7 + b = repr(m.Re()).replace('.', '')[:8] == repr(s.real()).replace('.', '')[:8] + c = (m.Im() - s.imag()).Abs() < 1e-7 + d = repr(m.Im()).replace('.', '')[:8] == repr(s.imag()).replace('.', '')[:8] + return (a or b) and (c or d) + +def elliptic_thetas(z, q): + ms = [] + for v in [ melliptic_theta(i, z, q) for i in range(1, 5) ]: + ms.append(v) + ss = [] + for v in [ elliptic_theta(i, z, CDF(q)) for i in range(1, 5) ]: + ss.append(v) + for m, s in zip(ms, ss): + assert rough_compare(m, s) + return zip(ms, ss) + +def elliptic_theta_primes(z, q): + ms = [] + for v in [ melliptic_theta_prime(i, z, q) for i in range(1, 5) ]: + ms.append(v) + ss = [] + for v in [ elliptic_theta_prime(i, z, CDF(q)) for i in range(1, 5) ]: + ss.append(v) + for m, s in zip(ms, ss): + assert rough_compare(m, s) + return zip(ms, ss) + +# def assert_sage_and_mathematica_agree_on_elliptic_theta(): +# I = CDF.gen() +# PI = CDF.pi() +# print "Asserting that sage and mathematica agree on elliptic_theta.." +# elliptic_thetas(1/2*I + 1.222, CDF(e**(-PI))) +# elliptic_thetas(10*I + 10, CDF(e**(-2*PI))) +# elliptic_thetas(1/3*I + 5, CDF(e**(-3*PI))) +# print "Asserting that sage and mathematica agree on elliptic_theta.. DONE" + +# def assert_sage_and_mathematica_agree_on_elliptic_theta_prime(): +# I = CDF.gen() +# PI = CDF.pi() +# print "Asserting that sage and mathematica agree on elliptic_theta_prime.." +# elliptic_theta_primes(1/2*I + 1.222, CDF(e**(-PI))) +# elliptic_theta_primes(10*I + 10, CDF(e**(-2*PI))) +# elliptic_theta_primes(1/3*I + 5, CDF(e**(-3*PI))) +# print "Asserting that sage and mathematica agree on elliptic_theta_prime.. DONE" + +def good_random_elliptic_theta_element_iterator(R=None): + if R is None: + R = CDF + II = R.gen() + PI = R.pi() + V = (QQ**2) + while True: + z = V.random_element() * vector([II, 1])/4 + q = ZZ.random_element(1, 10)/12 + II*ZZ.random_element(1, 10)/12 + q += II + + if z == 0 or z.abs() > 10 or q.abs() > 10: + continue + + q2 = (II*PI*q).exp() + yield z, q2 + +def assert_sage_and_mathematica_agree_on_elliptic_theta(n=5): + print "Asserting that sage and mathematica agree on elliptic_theta.." + it = good_random_elliptic_theta_element_iterator() + for i in range(n): + z, q = it.next() + elliptic_thetas(z, q) + print "Asserting that sage and mathematica agree on elliptic_theta.. DONE" + +def assert_sage_and_mathematica_agree_on_elliptic_theta_prime(n=5): + print "Asserting that sage and mathematica agree on elliptic_theta_prime.." + it = good_random_elliptic_theta_element_iterator() + for i in range(n): + z, q = it.next() + elliptic_theta_primes(z, q) + print "Asserting that sage and mathematica agree on elliptic_theta_prime.. DONE" + +def assert_sage_periods_for_elliptic_theta(n=5): + R = CDF + PI = R.pi() + II = R.gen() + print "Asserting that sage's elliptic_theta has correct periods.." + it = good_random_elliptic_theta_element_iterator(R=R) + for i in range(n): + z, q = it.next() + + v = elliptic_theta(1, z + PI, q)/elliptic_theta(1, z, q) + assert (v.real() + 1).abs() < 1e-10 + assert v.imag().abs() < 1e-10 + v = elliptic_theta(2, z + PI, q)/elliptic_theta(2, z, q) + assert (v.real() + 1).abs() < 1e-10 + assert v.imag().abs() < 1e-10 + v = elliptic_theta(3, z + PI, q)/elliptic_theta(3, z, q) + assert (v.real() - 1).abs() < 1e-10 + assert v.imag().abs() < 1e-10 + v = elliptic_theta(4, z + PI, q)/elliptic_theta(4, z, q) + assert (v.real() - 1).abs() < 1e-10 + assert v.imag().abs() < 1e-10 + + tau = q.log()/(II*PI) + auto = 1/q * (-2*II*z).exp() + v = elliptic_theta(1, z + tau*PI, q)/elliptic_theta(1, z, q) / -auto + assert (v - 1).abs() < 1e-10 + v = elliptic_theta(2, z + tau*PI, q)/elliptic_theta(2, z, q) / auto + assert (v - 1).abs() < 1e-10 + v = elliptic_theta(3, z + tau*PI, q)/elliptic_theta(3, z, q) / auto + assert (v - 1).abs() < 1e-10 + v = elliptic_theta(4, z + tau*PI, q)/elliptic_theta(4, z, q) / -auto + assert (v - 1).abs() < 1e-10 + + print "Asserting that sage's elliptic_theta has correct periods.. DONE" + +def assert_sage_and_mathematica_agree(n=5): + assert_sage_and_mathematica_agree_on_elliptic_theta(n=n) + assert_sage_and_mathematica_agree_on_elliptic_theta_prime(n=n) + assert_sage_periods_for_elliptic_theta(n=n) + +# # prettify the answer a little +# def cap(x): +# # indistinguishable from 0 at this precision +# R = x.parent() +# print x.abs().log(2), R(len(S)).log(2) + -(self._ring.prec() - 1) +# if x.abs().log(2) < R(len(S)).log(2) + -(self._ring.prec() - 1): +# return R(0) +# return x +# u, v = cap(u), cap(v) + +def good_random_siegel_theta_element_iterator(g, R=None): + if R is None: + R = CDF + II = R.gen() + PI = R.pi() + V = (QQ**g) + M = MatrixSpace(QQ, g, g) + while True: + D = diagonal_matrix([ ZZ.random_element(3, 20)/4 for _ in range(g) ]) + A = M.random_element() + while A.det() == 0: + A = M.random_element() + im = A * D * A.transpose() + + # the running time of these algorithms is dominated by the + # smallest eigenvalue of the imaginary part; here we artificially keep + # them large + mn = min(im.change_ring(CDF).eigenvalues()) + im = im / R(mn) + if not all(e > 1 for e in im.change_ring(CDF).eigenvalues()): + continue + + A = M.random_element() + re = A + A.transpose() + z = re + II*im + + p, q = V.random_element(), V.random_element() + u = z*p + q + yield z, u + +def assert_sage_and_mathematica_agree_on_siegel_theta(g, n=5): + print "Asserting that sage and mathematica agree on siegel_theta in genus %s.." % g + it = good_random_siegel_theta_element_iterator(g, R=CDF) + + cnt = 0 + mma_failures = 0 + while cnt < n: + z, u = it.next() + # the following strange procedure checks that values "around" + # the point particular point agree. What can happen is that a + # component (the imaginary component in my tests) is zero, + # which is returned only very approximately by mathematica. + try: + ms = [ mathematica.SiegelTheta(z, ZZ(i)/100*u) for i in range(98, 103) ] + except TypeError: + # mathematica fails to evaluate (with problems in some Range[] function) *alot* + mma_failures += 1 + if mma_failures > 5*n: + raise RuntimeError, "Mathematica has failed too often!" + continue + + ss = [ siegel_theta(z, ZZ(i)/100*u) for i in range(98, 103) ] + bad = [ rough_compare(m, s) for m, s in zip(ms, ss) ].count(False) + assert bad <= 1 + cnt += 1 + print "Asserting that sage and mathematica agree on siegel_theta in genus %s.. DONE" % g + +def abshift(u, z, a, b, r=ZZ(0), s=ZZ(0)): + R = z.base_ring() + II = R.gen() + PI = R.pi() + u = z.column_space()(u) + a = z.column_space()(a) + b = z.column_space()(b) + r = z.column_space()(r) + s = z.column_space()(s) + return (2*PI*II * (-ZZ(1)/2*a*z*a - a*u + r*b - s*a) ).exp() + +def assert_sage_periods_for_siegel_theta(g, n=10): + print "Asserting that sage's siegel_theta has correct periods in genus %s.." % g + it = good_random_siegel_theta_element_iterator(g, R=CDF) + cnt = 0 + while cnt < n: + z, u = it.next() + # the following strange procedure checks that values "around" + # the point particular point agree. What can happen is that a + # component (the imaginary component in my tests) is zero, + # which is returned only very approximately by mathematica. + orig = siegel_theta(z, u) + if orig.abs() < 1e-1 or orig.abs() > 1e10: # comparing such things is too error prone + continue + + try: + for a in (ZZ**g).basis(): + for b in (ZZ**g).basis(): + scale = abshift(u, z, a, b) + if scale.abs() < 1e-2 or scale.abs() > 1e20: # comparing such things is too error prone + raise ValueError + comp = siegel_theta(z, u + z*a + b) * ~scale + assert (orig - comp).abs() < 1e-5 or (orig/comp - 1) < 1e-5, \ + "I expected %s but got %s; difference is %s, ratio is %s" % (orig, comp, orig - comp, orig/comp) + cnt += 1 + except ValueError: + continue + print "Asserting that sage's siegel_theta has correct periods in genus %s.. DONE" % g + +def assert_sage_periods_for_siegel_theta_with_char(g, n=10): + print "Asserting that sage's siegel_theta_with_char has correct periods in genus %s.." % g + it = good_random_siegel_theta_element_iterator(g, R=ComplexField(100)) + cnt = 0 + skipped = 0 + + V = (ZZ**g) + characteristics = [] + characteristics += [ (V.random_element() / 2, V.random_element() / 2) for _ in range(3) ] + characteristics += [ (V.random_element() / 3, V.random_element() / 3) for _ in range(3) ] + + while cnt < n: + z, u = it.next() + # the following strange procedure checks that values "around" + # the point particular point agree. What can happen is that a + # component (the imaginary component in my tests) is zero, + # which is returned only very approximately by mathematica. + + r, s = choice(characteristics) + orig = siegel_theta_with_char((r, s), z, u) + if orig.abs() < 1e-3 or orig.abs() > 1e30: # comparing such things is too error prone + skipped += 1 + continue + + try: + for a in (ZZ**g).basis(): + for b in (ZZ**g).basis(): + scale = abshift(u, z, a, b, r, s) + if scale.abs() < 1e-3 or scale.abs() > 1e30: # comparing such things is too error prone + skipped += 1 + raise ValueError + comp = siegel_theta_with_char((r, s), z, u + z*a + b) * ~scale + assert (orig - comp).abs() < 1e-5 or (orig/comp - 1) < 1e-5, \ + "I expected %s but got %s; difference is %s, ratio is %s" % (orig, comp, orig - comp, orig/comp) + cnt += 1 + except ValueError: + continue + print "Asserting that sage's siegel_theta_with_char has correct periods in genus %s.. DONE (skipped %s)" % (g, skipped) + + +def num_deriv(f, u, v, h=ZZ(1)/100): + num = -f(u+2*h*v) + 8*f(u+h*v) - 8*f(u-h*v) + f(u-2*h*v) + den = 12*h + return num/den + +def num_second_deriv(f, u, v1, v2, h=ZZ(1)/100): + num = f(u - h*v1) - 2*f(u) + f(u + h*v2) + den = h**2 + return num/den + +def assert_sage_first_derivatives_for_siegel_theta(g, n=5): + print "Asserting that sage calculates siegel_theta first derivatives correctly in genus %s" % g + R = ComplexField(200) + TPII = 2 * R.gen() * R.pi() + it = good_random_siegel_theta_element_iterator(g, R=R) + cnt = 0 + while cnt < n: + z, u = it.next() +# # the following strange procedure checks that values "around" +# # the point particular point agree. What can happen is that a +# # component (the imaginary component in my tests) is zero, +# # which is returned only very approximately by mathematica. +# orig = siegel_theta(z, u) +# if orig.abs() < 1e-2 or orig.abs() > 1e10: # comparing such things is too error prone +# continue + cnt += 1 + + RT = RiemannTheta(z) + + vs1 = [] + vs2 = [] + vs3 = [] + for v in (QQ**g).basis(): + d = num_deriv(lambda t: RT.value_at_point(t), u, v, ZZ(1)/500) + vs1.append(d) + e = RT.value_at_point(u, derivs=[v]) + vs2.append(e) + m = maple.RiemannTheta(u.list(), z, [v.list()]) + vs3.append(maple_to_CDF(m)) + + grad_appx = vector(vs1) + grad_indiv = vector(vs2) + grad = RT.gradient_at_point(u) + grad_maple = vector(vs3) + + # we *really* can't expect much accuracy + e1 = (grad_appx - grad).norm() / grad_appx.norm() + e2 = (grad_appx - grad_indiv).norm() / grad_appx.norm() + e3 = (grad_appx - grad_maple).norm() / grad_appx.norm() + e4 = all(gr.norm() < 1e-5 for gr in [ grad_appx, grad_indiv, grad, grad_maple ]) + + if not (e4 or (e1 < 1e-5 and e2 < 1e-5 and e3 < 1e-5)): + s = "\n" + s += "%s \t- approximation\n" % grad_appx + s += "%s \t- maple\n" % grad_maple + s += "%s \t- individual (ratio %s)\n" % (grad_indiv, grad_indiv.norm() / grad_appx.norm()) + s += "%s \t- gradient (ratio %s)" % (grad, grad.norm() / grad_appx.norm()) + raise ValueError, s + + print "Asserting that sage calculates siegel_theta first derivatives correctly in genus %s.. DONE" % g + +def assert_sage_second_derivatives_for_siegel_theta(g, n=5): + print "Asserting that sage calculates siegel_theta second derivatives correctly in genus %s" % g + R = ComplexField(200) + TPII = 2 * R.gen() * R.pi() + it = good_random_siegel_theta_element_iterator(g, R=R) + cnt = 0 + while cnt < n: + z, u = it.next() + # the following strange procedure checks that values "around" + # the point particular point agree. What can happen is that a + # component (the imaginary component in my tests) is zero, + # which is returned only very approximately by mathematica. + orig = siegel_theta(z, u) + if orig.abs() < 1e-1 or orig.abs() > 1e6: # comparing such things is too error prone + continue + cnt += 1 + + RT = RiemannTheta(z) + + vs1 = [] + vs2 = [] + vs3 = [] + for v2 in (QQ**g).basis(): + for v1 in (QQ**g).basis(): + d = num_deriv(lambda t: RT.value_at_point(t, derivs=[v1]), u, v2, ZZ(1)/500) + # d = num_second_deriv(lambda t: RT.value_at_point(t), u, v1, v2, ZZ(1)/500) + vs1.append(d) + e = RT.value_at_point(u, derivs=[v1, v2]) + vs2.append(e) + m = maple.RiemannTheta(u.list(), z, [v1.list(), v2.list()]) + vs3.append(maple_to_CDF(m)) + + hess_appx = matrix(g, g, vs1) + hess_indiv = matrix(g, g, vs2) + hess = RT.hessian_at_point(u) + hess_maple = matrix(g, g, vs3) + + # we *really* can't expect much accuracy + e1 = (hess_appx - hess).norm() / hess_appx.norm() + e2 = (hess_appx - hess_indiv).norm() / hess_appx.norm() + e3 = (hess_appx - hess_maple).norm() / hess_appx.norm() + e4 = all(gr.norm() < 1e-5 for gr in [ hess_appx, hess_indiv, hess, hess_maple ]) + + if not (e4 or (e1 < 1e-5 and e2 < 1e-5 and e3 < 1e-5)): + s = "\n" + s += "%s \t- approximation\n" % hess_appx + s += "%s \t- maple\n" % hess_maple + s += ("%s \t- individual (ratio %s) (scalar %s)\n" + % (hess_indiv, hess_appx.norm() / hess_indiv.norm(), CDF(hess_indiv[0, 0] / hess_appx[0, 0]) )) + s += ("%s \t- hessian (ratio %s) (scalar %s)" + % (hess, hess_appx.norm() / hess.norm(), CDF(hess[0, 0] / hess_appx[0, 0]) )) + raise ValueError, s + + print "Asserting that sage calculates siegel_theta second derivatives correctly in genus %s.. DONE" % g + + +def assert_sage_third_derivatives_for_siegel_theta(g, n=5): + print "Asserting that sage calculates siegel_theta third derivatives correctly in genus %s" % g + R = ComplexField(200) + TPII = 2 * R.gen() * R.pi() + it = good_random_siegel_theta_element_iterator(g, R=R) + cnt = 0 + while cnt < n: + z, u = it.next() + # the following strange procedure checks that values "around" + # the point particular point agree. What can happen is that a + # component (the imaginary component in my tests) is zero, + # which is returned only very approximately by mathematica. + orig = siegel_theta(z, u) + if orig.abs() < 1e-1 or orig.abs() > 1e6: # comparing such things is too error prone + continue + cnt += 1 + + RT = RiemannTheta(z) + + vs1 = [] + vs2 = [] + vs3 = [] + D = {} + for v3_ in range(g): + for v2_ in range(g): + for v1_ in range(g): + kk = tuple(list(sorted([ v1_, v2_, v3_]))) + if D.has_key(kk): + continue + D[kk] = True + v1 = (QQ**g).basis()[v1_] + v2 = (QQ**g).basis()[v2_] + v3 = (QQ**g).basis()[v3_] + # print "kk=", kk + + m = maple.RiemannTheta(u.list(), z, [v1.list(), v2.list(), v3.list()]) + # m = maple.RiemannTheta(maple('[x, y]'), z).fdiff(maple(ss), maple('{x=%s, y=%s}' % tuple(u))) + vs2.append(maple_to_CDF(m)) + + e = RT.value_at_point(u, derivs=[v1, v2, v3]) + vs1.append(e) + E = {0:'x', 1:'y'} + ss = "'[" + ",".join([ E[v1_], E[v2_], E[v3_] ]) + "]'" + + vs1 = vector(vs1) + vs2 = vector(vs2) + e1 = (vs1 - vs2).norm() + e2 = (vs1 - vs2).norm() / vs2.norm() + if not (e1 < 1e-5 or e2 < 1e-5): + raise ValueError, "\n%s - calculated\n%s - maple\n%s - scalar\n%s - difference" % (vs1, vs2, vs1[0] / vs2[0], vs1-vs2) + + print "Asserting that sage calculates siegel_theta third derivatives correctly in genus %s.. DONE" % g + +def msiegel_theta_with_char(L, z, u): + return mathematica.SiegelTheta(L, z, u).N(16) + +def assert_automorphy_type1(n=5): + g = 2 # XXX + print "Asserting that siegel_theta is automorphic in genus %s.." % g + R = ComplexField(100) + TPII = 2 * R.gen() * R.pi() + it = good_random_siegel_theta_element_iterator(g, R=R) + cnt = 0 + while cnt < n: + z, u = it.next() + verbose("z = \n%s" % z, level=2) + verbose("u = \n%s" % u, level=2) + + for A in [ matrix(2, 2, [0, 1, -1, 0]), matrix(2, 2, [1, 1, 0, 1]) ]: + z_ = A * z * A.transpose() + u_ = A*u + v1 = shimura_phi(z, u) + v2 = shimura_phi(z_, u_) + + verbose("v1 = %s" % v1, level=1) + verbose("v2 = %s" % v2, level=1) + + e1 = v1.abs() < 1e-10 and v2.abs() < 1e-10 + e2 = (v1 - v2).abs() / v1.abs() < 1e-1 + assert e1 or e2 or (v1.abs().is_NaN() or v2.abs().is_NaN()) + + cnt += 1 + print "Asserting that siegel_theta is automorphic in genus %s.. DONE" % g + +def assert_automorphy_type2(n=5): + g = 2 # XXX + print "Asserting that siegel_theta is automorphic in genus %s.." % g + R = ComplexField(200) + TPII = 2 * R.gen() * R.pi() + it = good_random_siegel_theta_element_iterator(g, R=R) + cnt = 0 + while cnt < n: + z, u = it.next() + verbose("z = \n%s" % z, level=2) + verbose("u = \n%s" % u, level=2) + + M = MatrixSpace(ZZ, 2, 2).random_element() + B = M + M.transpose() + B = identity_matrix(2) + assert B.is_symmetric() + + z_ = z + B + u_ = u + v1 = shimura_phi(z, u) + v2 = shimura_phi_with_char([[0, 0], [ZZ(1)/2*B[0, 0], ZZ(1)/2*B[1, 1]]], z_, u_) + verbose("v1 = %s" % v1, level=1) + verbose("v2 = %s" % v2, level=1) + + e1 = v1.abs() < 1e-10 and v2.abs() < 1e-10 + e2 = (v1 / v2 - 1).abs() < 1e-10 + assert e1 or e2 or (v1.abs().is_NaN() or v2.abs().is_NaN()) + + cnt += 1 + print "Asserting that siegel_theta is automorphic in genus %s.. DONE" % g + +def assert_automorphy_type3(n=5): + g = 2 # XXX + print "Asserting that siegel_theta is automorphic in genus %s.." % g + R = ComplexField(250) + II = R.gen() + TPII = 2 * R.gen() * R.pi() + it = good_random_siegel_theta_element_iterator(g, R=R) + cnt = 0 + while cnt < n: + z, u = it.next() + # u = z.column_space()(0) + verbose("z = \n%s" % z, level=2) + verbose("u = \n%s" % u, level=2) + + z_ = -invert(z) + u_ = invert(z)*u + + if True: + v1 = shimura_phi(z, u) + v2 = shimura_phi(z_, u_) + else: + v1 = siegel_theta(z, u) + v2 = siegel_theta(z_, u_) + v1 = (TPII / 2 * u * invert(z) * u).exp() * v1 + + v1 = (-II*z).det().sqrt() * v1 + + verbose("v1 = %s" % v1, level=1) + verbose("v2 = %s" % v2, level=1) + rat = (v1 / v2) + verbose("ratio = %s" % CDF(rat), level=1) + verbose("ratio.abs() = %s" % CDF(rat.abs())) + + e1 = ((v1 / v2).abs() - 1).abs() + e2 = v1.abs() < 1e-10 and v2.abs() < 1e-10 + assert e1 or e2 or (v1.abs().is_NaN() or v2.abs().is_NaN()) + + cnt += 1 + print "Asserting that siegel_theta is automorphic in genus %s.. DONE" % g + +def log_shimura_phi(z, u): + R = z.base_ring() + II = R.gen() + PI = R.pi() + u = z.column_space()(u) + Z = z - z.apply_map(lambda t: t.conjugate()) + Z = invert(Z) + scale = (II * PI * (u * Z * u)) # .exp() + # print "scale", scale + return scale + siegel_theta(z, u).log() + +def log_shimura_phi_twisted(lattice, u): + return log_shimura_phi(lattice._tau, invert(lattice._w2) * u) + +class log_shimura_phi_twisted_symbolic: + def __init__(self, rs, lattice): + self._lattice = lattice + self._z = self._lattice._tau + self._R = self._z.base_ring() + II = self._R.gen() + PI = self._R.pi() + self._CS = self._z.column_space() + self._r, self._s = map(self._CS, rs) + + u0, u1 = var('u0, u1', ns=1) + u = vector([u0, u1]) + v = u * invert(self._lattice._w2).transpose() + + self._Z = self._z - self._z.apply_map(lambda t: t.conjugate()) + self._Z = invert(self._Z) + self._scale = (II * PI * (v * self._Z * v)) # .exp() + + # scale_ = R['u0, u1'](scale)(u0=u_[0], u1=u_[1]) + self._func = sage.symbolic.function.function('theta') + self._theta = self._func(*v) + self._total = self._theta.log() + self._scale + self._RT = sage.functions.riemann_theta.RiemannTheta(self._z) + + # cache with memory one + self._last_u = None + self._last_subs_dict = None + + def polynomial_for_derivatives(self, derivs): + assert len(derivs) >= 2 + + u0, u1 = var('u0, u1', ns=1) + u = vector([u0, u1]) + v = u * invert(self._lattice._w2).transpose() + func = self._func + theta = self._theta + + total_symb = self._total + for deriv in derivs: + total_symb = total_symb.derivative(u[deriv]) # .derivative(u0) # XXX + V0, V1 = var('V0, V1', ns=1) + XX0 = func(u0, u1).derivative(u0).subs(u0=V0, u1=V1).subs(V0=v[0], V1=v[1]) + XX1 = func(u0, u1).derivative(u1).subs(u0=V0, u1=V1).subs(V0=v[0], V1=v[1]) + total_symb = total_symb.subs(XX0 == var('X0', ns=1)) + total_symb = total_symb.subs(XX1 == var('X1', ns=1)) + + XX00 = func(u0, u1).derivative(u0).derivative(u0).subs(u0=V0, u1=V1).subs(V0=v[0], V1=v[1]) + XX01 = func(u0, u1).derivative(u0).derivative(u1).subs(u0=V0, u1=V1).subs(V0=v[0], V1=v[1]) + XX11 = func(u0, u1).derivative(u1).derivative(u1).subs(u0=V0, u1=V1).subs(V0=v[0], V1=v[1]) + total_symb = total_symb.subs(XX00 == var('X00', ns=1)) + total_symb = total_symb.subs(XX01 == var('X01', ns=1)) + total_symb = total_symb.subs(XX11 == var('X11', ns=1)) + + total_symb = total_symb.subs(theta.log() == var('logX', ns=1)) + total_symb = total_symb.subs(theta == var('X', ns=1)) + total_symb = total_symb.subs(u[0] == var('u0', ns=1)) + total_symb = total_symb.subs(u[1] == var('u1', ns=1)) + + S = self._R['u0, u1, logX, X, X0, X1, X00, X01, X11'].fraction_field() + total_poly = (total_symb)._polynomial_(S) + return total_poly + + def _subs_dict_for(self, u): + u_ = self._CS(u) + if self._last_u is not None and (self._last_u - u_).norm() < 1e-20: + return self._last_subs_dict + v_ = self._CS(invert(self._lattice._w2) * u_) + + v__ = v_ + self._z*self._r + self._s # XXX this is *not enough* if we're not taking second or higher derivatives! + X = self._RT.value_at_point(v__) + logX = X.log() + X0, X1 = self._RT.gradient_at_point(v__) + X00, X01, _, X11 = self._RT.hessian_at_point(v__).list() + + subs_dict = { 'u0':u_[0], 'u1':u_[1], 'logX':logX, 'X':X, 'X0':X0, 'X1':X1, 'X00':X00, 'X01':X01, 'X11':X11 } + self._last_u = u_ + self._last_subs_dict = subs_dict + return subs_dict + + def __call__(self, u, derivs=[0, 0]): + return self._R(self.polynomial_for_derivatives(tuple(derivs)).subs(**self._subs_dict_for(u))) + val_ = self._R( val ) + return val_ + +def automorphy_basis(): + Gs = [] + for A in [ matrix(2, 2, [0, 1, -1, 0]), matrix(2, 2, [1, 1, 0, 1]) ]: + G = block_diagonal_matrix([A, ~A.transpose()]) + Gs.append(G) + + # the form [1, b, 0, 1], b symmetric + for B in [ matrix(2, 2, [1, 0, 0, 0]), matrix(2, 2, [0, 0, 0, 1]), matrix(2, 2, [0, 1, 1, 0]) ]: + assert B.is_symmetric() + G = identity_matrix(4, 4) + G.set_block(0, 2, B) + # assert G in GammaTheta(2) + Gs.append(G) + + G = matrix(4, 4, [0, 0, 1, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, -1, 0, 0]) + Gs.append(G) + + assert all(G.det() == 1 for G in Gs) + assert all(is_symplectic(G) for G in Gs) + + return Gs diff -r 341d08e69736 -r d6e4600f2353 sage/functions/riemann_theta_tests_maple.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/riemann_theta_tests_maple.py Thu Jun 04 22:47:59 2009 -0700 @@ -0,0 +1,37 @@ +r""" +TESTS: + + Random functional tests for elliptic theta functions:: + + sage: from sage.functions.riemann_theta_tests import * + +# First derivatives:: + +# sage: assert_sage_first_derivatives_for_siegel_theta(2, n=2) # optional - maple, long time +# Asserting that sage calculates siegel_theta first derivatives correctly in genus 2 +# Asserting that sage calculates siegel_theta first derivatives correctly in genus 2.. DONE + +# sage: assert_sage_first_derivatives_for_siegel_theta(3, n=2) # optional - maple, long time +# Asserting that sage calculates siegel_theta first derivatives correctly in genus 3 +# Asserting that sage calculates siegel_theta first derivatives correctly in genus 3.. DONE + +# Second derivatives:: + +# sage: assert_sage_second_derivatives_for_siegel_theta(2, n=2) # optional - maple, long time +# Asserting that sage calculates siegel_theta second derivatives correctly in genus 2 +# Asserting that sage calculates siegel_theta second derivatives correctly in genus 2.. DONE + +# sage: assert_sage_second_derivatives_for_siegel_theta(3, n=2) # optional - maple, long time +# Asserting that sage calculates siegel_theta second derivatives correctly in genus 3 +# Asserting that sage calculates siegel_theta second derivatives correctly in genus 3.. DONE + +# Third derivatives:: + + sage: assert_sage_third_derivatives_for_siegel_theta(2, n=1) # optional - maple, long time + Asserting that sage calculates siegel_theta third derivatives correctly in genus 2 + Asserting that sage calculates siegel_theta third derivatives correctly in genus 2.. DONE + +AUTHORS: + +- Nick Alexander (2009-03) +""" diff -r 341d08e69736 -r d6e4600f2353 sage/functions/riemann_theta_tests_mathematica.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/riemann_theta_tests_mathematica.py Thu Jun 04 22:47:59 2009 -0700 @@ -0,0 +1,29 @@ +r""" +TESTS: + + Random functional tests for elliptic theta functions:: + + sage: from sage.functions.riemann_theta_tests import * + + sage: assert_sage_and_mathematica_agree_on_elliptic_theta(n=3) # optional - mathematica, long time + Asserting that sage and mathematica agree on elliptic_theta.. + Asserting that sage and mathematica agree on elliptic_theta.. DONE + + sage: assert_sage_and_mathematica_agree_on_elliptic_theta_prime(n=3) # optional - mathematica, long time + Asserting that sage and mathematica agree on elliptic_theta_prime.. + Asserting that sage and mathematica agree on elliptic_theta_prime.. DONE + + Random functional tests for siegel theta functions:: + + sage: assert_sage_and_mathematica_agree_on_siegel_theta(2, n=3) # optional - mathematica, long time + Asserting that sage and mathematica agree on siegel_theta in genus 2.. + Asserting that sage and mathematica agree on siegel_theta in genus 2.. DONE + + sage: assert_sage_and_mathematica_agree_on_siegel_theta(3, n=3) # optional - mathematica, long time + Asserting that sage and mathematica agree on siegel_theta in genus 3.. + Asserting that sage and mathematica agree on siegel_theta in genus 3.. DONE + +AUTHORS: + +- Nick Alexander (2009-03) +""" diff -r 341d08e69736 -r d6e4600f2353 sage/functions/riemann_theta_tests_sage.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/riemann_theta_tests_sage.py Thu Jun 04 22:47:59 2009 -0700 @@ -0,0 +1,51 @@ +r""" +TESTS: + + Random functional tests for elliptic theta functions:: + + sage: from sage.functions.riemann_theta_tests import * + + Periods:: + + sage: assert_sage_periods_for_elliptic_theta(n=10) # long time + Asserting that sage's elliptic_theta has correct periods.. + Asserting that sage's elliptic_theta has correct periods.. DONE + + Random functional tests for siegel theta functions: + + Periods:: + + sage: assert_sage_periods_for_siegel_theta(2, n=10) # long time + Asserting that sage's siegel_theta has correct periods in genus 2.. + Asserting that sage's siegel_theta has correct periods in genus 2.. DONE + + sage: assert_sage_periods_for_siegel_theta_with_char(2, n=5) # long time + Asserting that sage's siegel_theta_with_char has correct periods in genus 2.. + Asserting that sage's siegel_theta_with_char has correct periods in genus 2.. DONE (skipped ...) + + sage: assert_sage_periods_for_siegel_theta(3, n=10) # long time + Asserting that sage's siegel_theta has correct periods in genus 3 + Asserting that sage's siegel_theta has correct periods in genus 3.. DONE + + sage: assert_sage_periods_for_siegel_theta_with_char(3, n=4) # long time + Asserting that sage's siegel_theta_with_char has correct periods in genus 3.. + Asserting that sage's siegel_theta_with_char has correct periods in genus 3.. DONE (skipped ...) + + Automorphy:: + + sage: assert_automorphy_type1(n=10) + Asserting that siegel_theta is automorphic in genus 2.. + Asserting that siegel_theta is automorphic in genus 2.. DONE + + sage: assert_automorphy_type2(n=10) + Asserting that siegel_theta is automorphic in genus 2.. + Asserting that siegel_theta is automorphic in genus 2.. DONE + + sage: assert_automorphy_type3(n=10) + Asserting that siegel_theta is automorphic in genus 2.. + Asserting that siegel_theta is automorphic in genus 2.. DONE + +AUTHORS: + +- Nick Alexander (2009-03) +"""