diff -r f6e95f0166b4 sage/functions/riemann_theta.py
--- a/sage/functions/riemann_theta.py	Tue Apr 14 09:56:01 2009 -0700
+++ b/sage/functions/riemann_theta.py	Tue Apr 14 11:53:38 2009 -0700
@@ -450,6 +450,7 @@
         num_derivs = len(derivs)
         norm_derivs = [ deriv.norm() for deriv in derivs ]
         assert all( norm_deriv > 1e-10 for norm_deriv in norm_derivs ), "Directions need to be non-zero vectors"
+        prod_norm_derivs = prod( norm_derivs )
 
         U = self._T._pari_().qflll()._sage_() # find shortest vector in lattice
         v = (U * self._T).column(0)
@@ -463,7 +464,7 @@
             # 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 = -(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
@@ -578,13 +579,128 @@
         return exp_part, osc_part
 
     def value_at_point(self, z, derivs=[]):
+        r"""
+        EXAMPLES:
+            sage: from sage.functions.riemann_theta import RiemannTheta
+            sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I])
+            sage: E = RiemannTheta(1/10 * Omega)
+            sage: E.value_at_point(0)
+            4.41764896374 + 0.409494157593*I
+            sage: E.value_at_point([1/3, 2/3 + 1/2*I])
+            1582.31272298 + 996.188182401*I
+        """
         exp_part, osc_part = self.exp_and_osc_at_point(z, derivs=derivs)
         scale = 1
         if derivs:
             scale = z.base_ring().gen()**(len(derivs)-1) # experimentally derived fudge factor
         return scale * exp_part.exp() * osc_part
 
+    def derivatives_at_point(self, z, list_of_derivs=[]):
+        r"""
+        XXX
+        sage: from sage.functions.riemann_theta import RiemannTheta
+        sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I])
+        sage: E = RiemannTheta(1/10 * Omega)
+        sage: u = vector(CDF, [1/3, 2/3 + 1/2*I])
+        sage: X, X0, X1, X00, X01, X11 = E.derivatives_at_point(u, [ [(1, 0)], [(0, 1)], [(1, 0), (1, 0)], [(1, 0), (0, 1)], [(0, 1), (0, 1)] ])
+        sage: X
+        1582.31272298 + 996.188182401*I
+        sage: E.value_at_point(u)
+        1582.31272298 + 996.188182401*I
+        sage: vector([X0, X1])
+        (6869.0440534 - 11042.144279*I, 61131.0804751 - 49270.0715471*I)
+        sage: E.gradient_at_point(u)
+        (6869.0440534 - 11042.144279*I, 61131.0804751 - 49270.0715471*I)
+        sage: matrix(2, 2, [X00, X01, X01, X11])
+        [ 34666.7793688 - 139103.389251*I -481117.989651 + 81904.8136253*I]
+        [-481117.989651 + 81904.8136253*I  -1460873.9145 - 2861525.00205*I]
+        sage: E.hessian_at_point(u)
+        [ 34666.7793688 - 139103.389251*I -481117.989651 + 81904.8136253*I]
+        [-481117.989651 + 81904.8136253*I  -1460873.9145 - 2861525.00205*I]
+
+        sage: X, X1, X00, X0, X11, X01 = E.derivatives_at_point(u, [ [(0, 1)], [(1, 0), (1, 0)], [(1, 0)], [(0, 1), (0, 1)], [(1, 0), (0, 1)] ])
+        sage: X
+        1582.31272298 + 996.188182401*I
+        sage: vector([X0, X1])
+        (6869.0440534 - 11042.144279*I, 61131.0804751 - 49270.0715471*I)
+        sage: matrix(2, 2, [X00, X01, X01, X11])
+        [ 34666.7793688 - 139103.389251*I -481117.989651 + 81904.8136253*I]
+        [-481117.989651 + 81904.8136253*I  -1460873.9145 - 2861525.00205*I]
+        """
+        Pi = self._ring.pi()
+        I = self._ring.gen()
+        TPII = 2*Pi*I
+
+        # scale derivatives appropriately
+        CS = self._Omega.column_space()
+        for i in range(len(list_of_derivs)):
+            derivs = list_of_derivs[i]
+            derivs = [ TPII * CS(deriv) for deriv in derivs ]
+            list_of_derivs[i] = derivs
+        # scale input point appropriately
+        z = TPII * self._Omega.column_space()(z)
+
+        x = z.apply_map(lambda t: t.real())
+        y = z.apply_map(lambda t: t.imag())
+        shift = self._Xin * x
+        intshift = shift.apply_map(lambda t: ZZ(t.round()))
+
+        # the radius should be the largest product of norms *with the most terms in the product*
+        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
+
+        R = self.radius(max_elem)
+
+        verbose("Radius of S_R is %s" % R, level=2)
+        S = self.integer_points(z, R)
+        verbose("S_R consists of %s points" % len(S), level=2)
+        verbose("S_R = %s" % S, level=3)
+
+        uvs = finite_sum_with_many_derivatives(shift, intshift, x, y, self._X, self._Y, S, list_of_derivs, x.base_ring())
+        vals = []
+        for i in range(len(uvs)): # the first item is the value of the function at the point
+            u, v = uvs[i]
+
+            scale = 1
+            if i > 0:
+                derivs = list_of_derivs[i-1] # again, first item is value of function
+                scale = z.base_ring().gen()**(len(derivs)-1) # experimentally derived fudge factor
+            osc_part = u + I*v
+            exp_part = ZZ(1)/2 * (x * shift)
+            verbose("scale, exponential part, oscillating part = (%s, %s, %s)" % (scale, exp_part, osc_part), level=2)
+            vals.append(scale * exp_part.exp() * osc_part)
+        return vals
+
     def gradient_at_point(self, z):
+        r"""
+        EXAMPLES:
+            sage: from sage.functions.riemann_theta import RiemannTheta
+            sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I])
+            sage: E = RiemannTheta(1/10 * Omega)
+            sage: E.gradient_at_point([1/3, 2/3 + 1/2*I])
+            (6869.0440534 - 11042.144279*I, 61131.0804751 - 49270.0715471*I)
+            sage: E.gradient_at_point([1/3*I, 2/3*I + 1/2])
+            (706078145.06 - 649539106.278*I, 2333263039.98 - 412406488.368*I)
+        """
         Pi = self._ring.pi()
         I = self._ring.gen()
         TPII = 2*Pi*I
@@ -620,6 +736,18 @@
         return vector(self._ring, vals)
 
     def hessian_at_point(self, z):
+        r"""
+        EXAMPLES:
+            sage: from sage.functions.riemann_theta import RiemannTheta
+            sage: Omega = matrix(CDF, 2, 2, [1 + I, 2, 2, I])
+            sage: E = RiemannTheta(1/10 * Omega)
+            sage: E.hessian_at_point([1/3, 2/3 + 1/2*I])
+            [ 34666.7793688 - 139103.389251*I -481117.989651 + 81904.8136253*I]
+            [-481117.989651 + 81904.8136253*I  -1460873.9145 - 2861525.00205*I]
+            sage: E.hessian_at_point([1/3*I, 2/3*I + 1/2])
+            [-16895727870.4 - 11148594973.2*I -21658283152.2 - 41419586987.9*I]
+            [-21658283152.2 - 41419586987.9*I -3006810401.8 - 130403155418.0*I]
+        """
         Pi = self._ring.pi()
         I = self._ring.gen()
         TPII = 2*Pi*I