# HG changeset patch # User Nick Alexander # Date 1239241790 25200 # Node ID 81449f6b701126f8a6a56a34441547171e346804 # Parent 8f949a1a2b192f6375b9a86fb891bca3411c9866 [mq]: riemann_theta.patch diff -r 8f949a1a2b19 -r 81449f6b7011 sage/functions/all.py --- a/sage/functions/all.py Fri Feb 08 17:51:16 2008 -0800 +++ b/sage/functions/all.py Wed Apr 08 18:49:50 2009 -0700 @@ -40,3 +40,5 @@ i = I # alias from spike_function import spike_function + +from riemann_theta import elliptic_theta, elliptic_theta_prime, siegel_theta diff -r 8f949a1a2b19 -r 81449f6b7011 sage/functions/riemann_theta.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/riemann_theta.py Wed Apr 08 18:49:50 2009 -0700 @@ -0,0 +1,991 @@ +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 +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 + +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) + +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 upper_gamma_1(a, b): + return a.gamma_inc(b) + +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) + + 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) + + 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" + + 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() + + # 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()) + + from sage.numerical.optimize import find_root + try: + rt = find_root(func, left, self._ring.prec()*left) + except RuntimeError: + rt = -left + + 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) + + 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())) + + 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=[]): + 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 gradient_at_point(self, z): + 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 + 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): + 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 + 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 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])) # optional - mathematica + 1.0017142124239051 + 0.*I + 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.222*I])) + 0.99494981531800369821577755005 + 0.0015359508886263184979475482893*I + sage: mathematica.SiegelTheta(tau3, vector([0, 0, 0])) # optional - mathematica + 0.9983074158131862 + 1...-18*I + sage: mathematica.SiegelTheta(tau3, 1/10*vector([1/2 + I, 2/3*I, 1.222*I])) # optional - mathematica + 0.9949498153180014 + 0.001535950888626739*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.222*I])) + 0.99494981531800369821577755004768563987129755626530324502775 + 0.0015359508886263184979475482892828605395586367866932565119700*I + """ + return RiemannTheta(Omega).value_at_point(z) diff -r 8f949a1a2b19 -r 81449f6b7011 sage/functions/riemann_theta_tests.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/functions/riemann_theta_tests.py Wed Apr 08 18:49:50 2009 -0700 @@ -0,0 +1,575 @@ +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 elliptic_theta, elliptic_theta_prime, siegel_theta +from sage.functions.riemann_theta import 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 + +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 = 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 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) + for i in range(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 + + 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) + 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): + 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 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) + +def assert_automorphy_type1(n=5): + g = 2 # XXX + print "Asserting that siegel_theta is automorphic in genus %s.. DONE" % 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.. DONE" % 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 diff -r 8f949a1a2b19 -r 81449f6b7011 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 Wed Apr 08 18:49:50 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 8f949a1a2b19 -r 81449f6b7011 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 Wed Apr 08 18:49:50 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 8f949a1a2b19 -r 81449f6b7011 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 Wed Apr 08 18:49:50 2009 -0700 @@ -0,0 +1,43 @@ +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(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 + + 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.. 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.. 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.. DONE + +AUTHORS: + +- Nick Alexander (2009-03) +"""