diff -r c1f419157c16 -r f6e95f0166b4 sage/functions/riemann_theta.py --- a/sage/functions/riemann_theta.py Sat Apr 11 15:07:08 2009 -0700 +++ b/sage/functions/riemann_theta.py Tue Apr 14 09:56:01 2009 -0700 @@ -238,7 +238,9 @@ return self._Omega.columns() def upper_gamma_1(a, b): - return a.gamma_inc(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): @@ -279,6 +281,10 @@ verbose("Cholesky decomposition is:", level=2) verbose(self._T, level=2) + # cached + 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() @@ -453,36 +459,60 @@ d = ZZ(1) / min( self._T[i, i] for i in range(0, g) ) detT = self._T.det() - # 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 - RRR = ComplexField((ZZ(3)/2 * self._ring.prec()).ceil()) - # 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()) + 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) - from sage.numerical.optimize import find_root - try: - rt = find_root(func, left, self._ring.prec()*left) - except RuntimeError: - rt = -left + left = (RRR(g + 2*num_derivs + (RRR(g)**2 + 8*num_derivs).sqrt()).sqrt() + r)/2 + left = RDF(left.real()) - R = max(left, rt) - assert func(R) < 2**(-20), "Could not find an accurate bound for the radius R" - R += ZZ(1)/100 # plus a slight epsilon, to deal with errors in root finding to low precision - return RDF(R) + 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""" @@ -571,7 +601,10 @@ 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 + # 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) @@ -607,7 +640,10 @@ 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 + # 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) diff -r c1f419157c16 -r f6e95f0166b4 sage/functions/riemann_theta_tests.py --- a/sage/functions/riemann_theta_tests.py Sat Apr 11 15:07:08 2009 -0700 +++ b/sage/functions/riemann_theta_tests.py Tue Apr 14 09:56:01 2009 -0700 @@ -238,26 +238,34 @@ return (2*PI*II * (-ZZ(1)/2*a*z*a - a*u) ).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 + print "Asserting that sage's siegel_theta has correct periods in genus %s.." % g it = good_random_siegel_theta_element_iterator(g, R=CDF) - for i in range(n): + 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 - 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() > 1e10: # comparing such things is too error prone - continue - 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) + 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 + 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, \ + "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 num_deriv(f, u, v, h=ZZ(1)/100): @@ -477,7 +485,7 @@ def assert_automorphy_type1(n=5): g = 2 # XXX - print "Asserting that siegel_theta is automorphic in genus %s.. DONE" % g + 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) @@ -505,7 +513,7 @@ def assert_automorphy_type2(n=5): g = 2 # XXX - print "Asserting that siegel_theta is automorphic in genus %s.. DONE" % g + 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) @@ -536,7 +544,7 @@ def assert_automorphy_type3(n=5): g = 2 # XXX - print "Asserting that siegel_theta is automorphic in genus %s" % g + print "Asserting that siegel_theta is automorphic in genus %s.." % g R = ComplexField(250) II = R.gen() TPII = 2 * R.gen() * R.pi() @@ -573,3 +581,192 @@ 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 + + @cached_method + 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 + + + +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 c1f419157c16 -r f6e95f0166b4 sage/functions/riemann_theta_tests_sage.py --- a/sage/functions/riemann_theta_tests_sage.py Sat Apr 11 15:07:08 2009 -0700 +++ b/sage/functions/riemann_theta_tests_sage.py Tue Apr 14 09:56:01 2009 -0700 @@ -26,15 +26,15 @@ Automorphy:: sage: assert_automorphy_type1(n=10) - Asserting that siegel_theta is automorphic in genus 2.. DONE + 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.. DONE + 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.. DONE + Asserting that siegel_theta is automorphic in genus 2.. Asserting that siegel_theta is automorphic in genus 2.. DONE AUTHORS: