diff -r a5e0c4d1c443 -r 37aa26cf30a6 sage/functions/riemann_theta.py
--- a/sage/functions/riemann_theta.py	Tue Apr 28 22:28:22 2009 -0700
+++ b/sage/functions/riemann_theta.py	Tue May 05 11:14:25 2009 -0700
@@ -42,11 +42,16 @@
 """
 
 from sage.ext.fast_callable import fast_callable
-from sage.rings.all import RDF, CDF, RealField, ComplexField, ZZ, QQ
+from sage.rings.all import RDF, CDF, RealField, ComplexField, ZZ, QQ, GF
 from sage.misc.misc import verbose
 from sage.misc.misc_c import prod # sum is built in
 from sage.modules.free_module_element import vector
 from sage.matrix.constructor import matrix, identity_matrix
+from sage.parallel.decorate import parallel
+
+@parallel()
+def _global_theta_for_parallel(characteristics, RT, u):
+    return RT.siegel_theta_with_char(characteristics, u)
 
 def finite_sum_without_derivatives(shift, intshift, x, y, X, Y, S, domain):
     r"""
@@ -814,6 +819,105 @@
         characteristics = (CS(0), CS(0))
         return self.shimura_phi_with_char(characteristics, u)
 
+    def theta_dict(self, u, list_of_characteristics=None):
+        r"""
+            sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6])
+            sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2)
+            sage: D = RT.theta_dict(0, [ ((0, 0), (0, 0)), ((0, 1/2), (0, 1/2)) ])
+            sage: D[((0, 0), (0, 0))]
+            1.00171421242
+            sage: D[((0, 1/2), (0, 1/2))].abs() < 1e-15
+            True
+
+            sage: RT.theta_dict([1/2, 1/5 + CDF(I)*0.234], [ ((0, 1/2), (0, 1/2)), ((0, 0), (0, 0)) ])
+            {((0, 0), (0, 0)): 0.998788360598 + 0.00120506039896*I, ((0, 1/2), (0, 1/2)): -0.136735147549 - 0.117846448587*I}
+
+            sage: RT.theta_vector([1/2, 1/5 + CDF(I)*0.234], [ ((0, 0), (0, 0)), ((0, 1/2), (0, 1/2)) ])
+            (0.998788360598 + 0.00120506039896*I, -0.136735147549 - 0.117846448587*I)
+            sage: RT.theta_vector([1/2, 1/5 + CDF(I)*0.234], [ ((0, 1/2), (0, 1/2)), ((0, 0), (0, 0)) ])
+            (-0.136735147549 - 0.117846448587*I, 0.998788360598 + 0.00120506039896*I)
+        """
+        if list_of_characteristics is None:
+            list_of_characteristics = list(self._all_characteristics())
+        L = list(_global_theta_for_parallel([ (c, self, u) for c in list_of_characteristics ]))
+        D = {}
+        for (args, kwds), val in L:
+            D[args[0]] = val
+        return D
+
+    def even_theta_dict(self, u):
+        return self.theta_dict(u, self._even_characteristics())
+
+    def odd_theta_dict(self, u):
+        return self.theta_dict(u, self._odd_characteristics())
+
+    def even_theta_vector(self, u):
+        return self.theta_vector(u, self._even_characteristics())
+
+    def odd_theta_vector(self, u):
+        return self.theta_vector(u, self._odd_characteristics())
+
+    def theta_vector(self, u, list_of_characteristics=None):
+        r"""
+            sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6])
+            sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2)
+            sage: RT.theta_vector(0)
+            (1.00171421242, 0.998608586772 + 2.06043088621e-36*I, 0.999838600402 + 6.56441798402e-36*I, 0.999838600402 + 1.11674199318e-36*I, 0.283260714578, 1.73435525065e-17 - 6.52327288837e-17*I, 0.278207060819 + 2.68012527741e-36*I, 1.7344656162e-17 - 4.25398925609e-17*I, 0.197753422929, 0.181367479122 - 3.91422304159e-34*I, 1.21087425915e-17 + 8.00949937526e-17*I, 1.21087425915e-17 + 8.00979583156e-17*I, 0.283260714578, 1.73435525065e-17 + 5.02857477884e-16*I, 1.7344656162e-17 - 5.08846937396e-16*I, 0.278207060819 - 4.93038065763e-32*I)
+        """
+        if list_of_characteristics is None:
+            list_of_characteristics = list(self._all_characteristics())
+        D = self.theta_dict(u, list_of_characteristics)
+
+        L = []
+        for c in list_of_characteristics:
+            L.append(D[c])
+        return vector(L)
+
+    def _all_characteristics(self):
+        r"""
+            sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6])
+            sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2)
+            sage: list(RT._all_characteristics())
+            [((0, 0), (0, 0)), ((0, 0), (1/2, 0)), ((0, 0), (0, 1/2)), ((0, 0), (1/2, 1/2)),
+             ((1/2, 0), (0, 0)), ((1/2, 0), (1/2, 0)), ((1/2, 0), (0, 1/2)), ((1/2, 0), (1/2, 1/2)),
+             ((0, 1/2), (0, 0)), ((0, 1/2), (1/2, 0)), ((0, 1/2), (0, 1/2)), ((0, 1/2), (1/2, 1/2)),
+             ((1/2, 1/2), (0, 0)), ((1/2, 1/2), (1/2, 0)), ((1/2, 1/2), (0, 1/2)), ((1/2, 1/2), (1/2, 1/2))]
+        """
+        for u in GF(2)**(self._g):
+            for v in GF(2)**(self._g):
+                yield (tuple(u.change_ring(ZZ)/2), tuple(v.change_ring(ZZ)/2))
+
+    def _even_characteristics(self):
+        r"""
+            sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6])
+            sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2)
+            sage: list(RT._even_characteristics())
+            [((0, 0), (0, 0)), ((0, 0), (1/2, 0)), ((0, 0), (0, 1/2)), ((0, 0), (1/2, 1/2)),
+             ((1/2, 0), (0, 0)), ((1/2, 0), (0, 1/2)), ((0, 1/2), (0, 0)), ((0, 1/2), (1/2, 0)),
+             ((1/2, 1/2), (0, 0)), ((1/2, 1/2), (1/2, 1/2))]
+        """
+        for u in GF(2)**(self._g):
+            for v in GF(2)**(self._g):
+                u = u.change_ring(ZZ)
+                v = v.change_ring(ZZ)
+                if u*v % 2 == 0:
+                    yield (tuple(u/2), tuple(v/2))
+
+    def _odd_characteristics(self):
+        r"""
+            sage: tau2 = CDF(I)*1/2*matrix(CDF, 2, 2, [5, 3, 3, 6])
+            sage: RT = sage.functions.riemann_theta.RiemannTheta(tau2)
+            sage: list(RT._odd_characteristics())
+            [((1/2, 0), (1/2, 0)), ((1/2, 0), (1/2, 1/2)), ((0, 1/2), (0, 1/2)),
+             ((0, 1/2), (1/2, 1/2)), ((1/2, 1/2), (1/2, 0)), ((1/2, 1/2), (0, 1/2))]
+        """
+        for u in GF(2)**(self._g):
+            for v in GF(2)**(self._g):
+                u = u.change_ring(ZZ)
+                v = v.change_ring(ZZ)
+                if u*v % 2 == 1:
+                    yield (tuple(u/2), tuple(v/2))
+
 def elliptic_theta(n, z, q):
     r"""
     Compute Jacobi elliptic theta functions.