diff -r b4439a628ed2 -r a5e0c4d1c443 sage/functions/riemann_theta.py
--- a/sage/functions/riemann_theta.py	Fri Apr 24 20:58:34 2009 -0700
+++ b/sage/functions/riemann_theta.py	Tue Apr 28 22:28:22 2009 -0700
@@ -257,6 +257,7 @@
         verbose(self._T, level=2)
 
         # cached
+        self._no_derivs_radius = None
         self._gradient_radius = None
         self._hessian_radius = None
 
@@ -535,7 +536,12 @@
         shift = self._Xin * x
         intshift = shift.apply_map(lambda t: ZZ(t.round()))
 
-        R = self.radius(derivs)
+        if not derivs:
+            if self._no_derivs_radius is None:
+                self._no_derivs_radius = self.radius(derivs)
+            R = self._no_derivs_radius
+        else:
+            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)
@@ -774,6 +780,40 @@
                 m[i, j] = m[j, i]
         return m
 
+    def siegel_theta(self, u):
+        return self.value_at_point(u)
+
+    def siegel_theta_with_char(self, characteristics, u):
+        r, s = characteristics
+        CS = self._Omega.column_space()
+        r = CS(r)
+        s = CS(s)
+        u = CS(u)
+
+        R = self._ring
+        II = R.gen()
+        PI = R.pi()
+        z = (self._Omega/(-2*PI*II))
+        scale = (2*PI*II * ( ZZ(1)/2*r*z*r + r*(u + s) )).exp()
+        u_ = u + z*r + s
+        return scale * self.siegel_theta(u_)
+
+    def shimura_phi_with_char(self, characteristics, u):
+        R = self._ring
+        II = R.gen()
+        PI = R.pi()
+        z = (self._Omega/(-2*PI*II))
+        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 * self.siegel_theta_with_char(characteristics, u)
+
+    def shimura_phi(self, z, u):
+        CS = self._Omega.column_space()
+        characteristics = (CS(0), CS(0))
+        return self.shimura_phi_with_char(characteristics, u)
+
 def elliptic_theta(n, z, q):
     r"""
     Compute Jacobi elliptic theta functions.
@@ -1129,7 +1169,7 @@
     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)
+    return RiemannTheta(Omega).siegel_theta(z)
 
 def siegel_theta_with_char(characteristics, z, u):
     r"""                                                                                                                                
@@ -1173,28 +1213,10 @@
     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_)
+    return RiemannTheta(z).siegel_theta_with_char(characteristics, 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)
+    return RiemannTheta(z).shimura_phi_with_char(characteristics, u)
 
 def shimura_phi(z, u):
-    return shimura_phi_with_char([0, 0], z, u)
+    return RiemannTheta(z).shimura_phi(u)