diff -r 6f3e19a9467b -r 03b9d8113657 sage/functions/all.py --- a/sage/functions/all.py Tue Apr 14 11:53:38 2009 -0700 +++ b/sage/functions/all.py Thu Apr 23 15:19:33 2009 -0700 @@ -41,4 +41,8 @@ from spike_function import spike_function -from riemann_theta import elliptic_theta, elliptic_theta_prime, siegel_theta +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 diff -r 6f3e19a9467b -r 03b9d8113657 sage/functions/complex_lattice.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/complex_lattice.py Thu Apr 23 15:19:33 2009 -0700 @@ -0,0 +1,108 @@ +# -*- sage -*- + +import sage.matrix.matrix +from sage.rings.all import ZZ, QQ, ComplexField +from riemann_theta import invert + +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 6f3e19a9467b -r 03b9d8113657 sage/functions/klein_p.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/klein_p.py Thu Apr 23 15:19:33 2009 -0700 @@ -0,0 +1,355 @@ +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 + +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] + + + + + + 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: siegel_theta_with_char(c, T._lattice._tau, 0) + sage: T._rational_function_for_derivative(0, 0) + + sage: [ T([ 1/(10*n^4), 0 ]) for n in range(1, 10, 2) ] + + # sage: u = 0*u + + sage: w2i = sage.functions.riemann_theta.invert(T._lattice._w2) + sage: v0 = siegel_theta_with_char(c, T._lattice._tau, w2i*u) + sage: vs0 = [ siegel_theta_with_char(c, T._lattice._tau, w2i*(u + Omega*v))/v0 for v in (QQ^4).basis() ]; vs0 + + sage: T(u, theta=True) + sage: T(u) * siegel_theta_with_char(c, T._lattice._tau, w2i*u)^2 + + sage: v1 = T(u, theta=True); v1 + sage: T._last_subs_dict['X'] + 4.7152730063487808517364990257033652942008921406590784528275e-6 - 3.1990704034116803052447891933386212484433501980263454265724e-6*I + sage: siegel_theta(T._lattice._tau, w2i*u + T._z*T._r + T._s) + 4.7152730063487808517364990257033652942008921406590784528275e-6 - 3.1990704034116803052447891933386212484433501980263454265724e-6*I + + sage: vs1 = [ T((u + Omega*v), theta=True)/v1 for v in (QQ^4).basis() ]; vs1 + + sage: 2 + """ + def __init__(self, rs, lattice): + self._lattice = ComplexLattice(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 = 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._total = -1 * (self._theta.log() + self._scale) + 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 _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 + + @cached_method + def _rational_function_for_derivative(self, *derivs): + r""" + Calculate the rational function giving $p_{derivs}$ in terms of theta + derivatives. + """ + assert len(derivs) >= 2 + + u0, u1 = myvar('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]) + + # 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 = 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(theta.log() == myvar('logX', ns=1)) # important that log is substitutated + total_symb = total_symb.subs(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] + # one last bit of hijinks, to avoid computing GCDs over inexact rings + total_symb = (X**len(derivs) * total_symb).expand() + total_poly = (total_symb)._polynomial_(S) / S.gen(0)**len(derivs) + 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=False): + 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)) + if theta: + poly = poly.numerator() + 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=False): + 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 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 + +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): + d = act_on_char(vector(flatten(c)), G) + Omega_ = T._lattice._Omega * G.transpose() + d1, d2, d3, d4 = d + T2 = KleinPFunction([(d1, d2), (d3, d4)], Omega_) + 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 6f3e19a9467b -r 03b9d8113657 sage/functions/riemann_theta.py --- a/sage/functions/riemann_theta.py Tue Apr 14 11:53:38 2009 -0700 +++ b/sage/functions/riemann_theta.py Thu Apr 23 15:19:33 2009 -0700 @@ -207,35 +207,10 @@ else: return finite_sum_without_derivatives(shift, intshift, x, y, X, Y, S, domain) -class ComplexLattice: - def __init__(self, tau): - self._tau = tau - self._C = tau.column_space() - self._g = tau.nrows() - assert tau.ncols() == self._g - self._Omega = self._tau.augment(identity_matrix(self._g)) - - 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 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 periods(self): - return self._Omega.columns() +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) @@ -285,11 +260,6 @@ self._gradient_radius = None self._hessian_radius = None - 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)) - self._Xin = invert(self._X) self._Tin = invert(self._T) @@ -626,6 +596,10 @@ 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() @@ -646,27 +620,30 @@ 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* - 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] + 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 + # 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) @@ -1153,3 +1130,71 @@ 0.99494981531800369821577755004768563987129755626530324502775 + 0.0015359508886263184979475482892828605395586367866932565119700*I """ return RiemannTheta(Omega).value_at_point(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])) # optional - mathematica + 1.0017142124239051 + 0.*I + + 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])) # optional - mathematica + 2.9891194456216093*^-19 + 3.1440919555182356*^-19*I + """ + r, s = characteristics + CS = z.column_space() + r = CS(r) + s = CS(s) + u = CS(u) + u_ = u + z*r + s + + R = z.base_ring() + II = R.gen() + PI = R.pi() + scale = (2*PI*II * (+ZZ(1)/2*r*z*r + r*(u + s) )).exp() + return scale * siegel_theta(z, u_) + +def shimura_phi_with_char(characteristics, 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() + return scale * siegel_theta_with_char(characteristics, z, u) + +def shimura_phi(z, u): + return shimura_phi_with_char([0, 0], z, u) diff -r 6f3e19a9467b -r 03b9d8113657 sage/functions/riemann_theta_tests.py --- a/sage/functions/riemann_theta_tests.py Tue Apr 14 11:53:38 2009 -0700 +++ b/sage/functions/riemann_theta_tests.py Thu Apr 23 15:19:33 2009 -0700 @@ -10,8 +10,8 @@ - Nick Alexander (2009-03) """ -from sage.functions.all import elliptic_theta, elliptic_theta_prime, siegel_theta -from sage.functions.riemann_theta import RiemannTheta +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 @@ -228,14 +228,16 @@ cnt += 1 print "Asserting that sage and mathematica agree on siegel_theta in genus %s.. DONE" % g -def abshift(u, z, a, b): +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) - return (2*PI*II * (-ZZ(1)/2*a*z*a - a*u) ).exp() + 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 @@ -243,14 +245,12 @@ cnt = 0 while cnt < n: z, u = it.next() - print "." # 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 - print "XXX", cnt continue try: @@ -258,7 +258,6 @@ 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 - print "SCALE TOO LARGE", scale raise ValueError comp = siegel_theta(z, u + z*a + b) * ~scale assert (orig - comp).abs() < 1e-5 or (orig/comp - 1) < 1e-5, \ @@ -268,6 +267,46 @@ 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 @@ -438,48 +477,6 @@ print "Asserting that sage calculates siegel_theta third derivatives correctly in genus %s.. DONE" % g -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 shimura_phi_with_char(L, 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_with_char(L, z, u) - -def shimura_phi(z, u): - return shimura_phi_with_char([0, 0], z, u) - -def abshift(u, z, a, b): - 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) - return (2*PI*II * (-ZZ(1)/2*a*z*a - a*u) ).exp() - -def siegel_theta_with_char(L, z, u): - r, s = L - CS = z.column_space() - r = CS(r) - s = CS(s) - u = CS(u) - u_ = u + z*r + s - - R = z.base_ring() - II = R.gen() - PI = R.pi() - scale = (2*PI*II * (+ZZ(1)/2*r*z*r + r*(u + s) )).exp() - return scale * siegel_theta(z, u_) - def msiegel_theta_with_char(L, z, u): return mathematica.SiegelTheta(L, z, u).N(16) @@ -624,7 +621,6 @@ self._last_u = None self._last_subs_dict = None - @cached_method def polynomial_for_derivatives(self, derivs): assert len(derivs) >= 2 @@ -702,71 +698,3 @@ assert all(is_symplectic(G) for G in Gs) return Gs - - - -def assert_automorphy_type1_log(n=5): - r""" - sage: from sage.functions.riemann_theta_tests import assert_automorphy_type1_log - sage: assert_automorphy_type1_log(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 - - c = Chr([ZZ(0), ZZ(1)/2, 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) - - w1, w2 = z, 0*z + 1 - L = ComplexLattice(w1, w2) - - try: - v1 = log_shimura_phi_twisted_symbolic(c.characteristics(), L)(u) - failed = False - - if v1.abs() > 1e30: - raise RuntimeError - -# for A in [ matrix(2, 2, [0, 1, -1, 0]), matrix(2, 2, [1, 1, 0, 1]) ]: -# assert A.det() == 1 -# G = block_diagonal_matrix([A, ~A.transpose()]) - for G in automorphy_basis(): - d = c.acted_on_by(G) # .reduced() - # print "c=", c - # print "d=", d - - Omega_ = L._Omega * G.transpose() - w1_, w2_ = Omega_.submatrix(0, 0, 2, 2), Omega_.submatrix(0, 2, 2, 2) - L_ = ComplexLattice(w1_, w2_) - v2 = log_shimura_phi_twisted_symbolic(d.characteristics(), L_)(u) - # v2 = log_shimura_phi_twisted(L_, u) - # v2 = log_shimura_phi(z_, u_) - - print cap(v2) - # 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() < 1e-6 - - if not (e1 or e2 or (v1.abs().is_NaN() or v2.abs().is_NaN())): - print "\nv1=%s\nv2=%s\ne2=%s" % (v1, v2, e2) - failed = True - - # assert e1 or e2 or (v1.abs().is_NaN() or v2.abs().is_NaN()), "\nv1=%s\nv2=%s\ne2=%s" % (v1, v2, e2) - print cap(v1) - print - assert not failed - cnt += 1 - except RuntimeError: - continue - print "Asserting that siegel_theta is automorphic in genus %s.. DONE" % g - -assert_automorphy_type1_log() diff -r 6f3e19a9467b -r 03b9d8113657 sage/functions/riemann_theta_tests_sage.py --- a/sage/functions/riemann_theta_tests_sage.py Tue Apr 14 11:53:38 2009 -0700 +++ b/sage/functions/riemann_theta_tests_sage.py Thu Apr 23 15:19:33 2009 -0700 @@ -19,10 +19,18 @@ 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 + 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 + Automorphy:: sage: assert_automorphy_type1(n=10)