diff --git a/healpy/sphtfunc.py b/healpy/sphtfunc.py
index 6a4b515d0..4bb49d299 100644
--- a/healpy/sphtfunc.py
+++ b/healpy/sphtfunc.py
@@ -347,7 +347,7 @@ def alm2map(
         return np.array(output)
 
 
-def synalm(cls, lmax=None, mmax=None, new=False, verbose=True):
+def synalm(cls, lmax=None, mmax=None, new=False, rng=None, verbose=True):
     """Generate a set of alm given cl.
     The cl are given as a float array. Corresponding alm are generated.
     If lmax is None, it is assumed lmax=cl.size-1
@@ -369,6 +369,10 @@ def synalm(cls, lmax=None, mmax=None, new=False, verbose=True):
       (e.g. TT, EE, BB, TE, EB, TB or TT, EE, BB, TE if 4 cl as input).
       If False, use the old ordering, ie by row
       (e.g. TT, TE, TB, EE, EB, BB or TT, TE, EE, BB if 4 cl as input).
+      
+    rng: numpy random number generator state. Useful when generating multiple
+       maps in parallel. Typically, `rng=np.random.RandomState(seed)`
+       
 
     Returns
     -------
@@ -408,8 +412,14 @@ def synalm(cls, lmax=None, mmax=None, new=False, verbose=True):
         cls_list = [np.asarray(cls, dtype=np.float64)]
         szalm = Alm.getsize(lmax, mmax)
         alm = np.zeros(szalm, "D")
-        alm.real = np.random.standard_normal(szalm)
-        alm.imag = np.random.standard_normal(szalm)
+        
+        if rng is None:
+            alm.real = np.random.standard_normal(szalm)
+            alm.imag = np.random.standard_normal(szalm)
+        else:
+            alm.real = rng.standard_normal(szalm)
+            alm.imag = rng.standard_normal(szalm)
+            
         alms_list = [alm]
         sphtlib._synalm(cls_list, alms_list, lmax, mmax)
         return alm
@@ -442,8 +452,14 @@ def synalm(cls, lmax=None, mmax=None, new=False, verbose=True):
     alms_list = []
     for i in six.moves.xrange(Nspec):
         alm = np.zeros(szalm, "D")
-        alm.real = np.random.standard_normal(szalm)
-        alm.imag = np.random.standard_normal(szalm)
+        
+        if rng is None:
+            alm.real = np.random.standard_normal(szalm)
+            alm.imag = np.random.standard_normal(szalm)
+        else:
+            alm.real = rng.standard_normal(szalm)   
+            alm.imag = rng.standard_normal(szalm)
+        
         alms_list.append(alm)
     if new:  # new input order: input given by diagonal, should be given by row
         cls_list = new_to_old_spectra_order(cls_list)
@@ -467,6 +483,7 @@ def synfast(
     fwhm=0.0,
     sigma=None,
     new=False,
+    rng=None,
     verbose=True,
 ):
     """Create a map(s) from cl(s).
@@ -497,6 +514,7 @@ def synfast(
     sigma : float, scalar, optional
       The sigma of the Gaussian used to smooth the map (applied on alm)
       [in radians]
+    rng: Numpy random number generator state, passed to synalm.
 
     Returns
     -------
@@ -524,7 +542,7 @@ def synfast(
     cls_lmax = cb.len_array_or_arrays(cls) - 1
     if lmax is None or lmax < 0:
         lmax = min(cls_lmax, 3 * nside - 1)
-    alms = synalm(cls, lmax=lmax, mmax=mmax, new=new, verbose=verbose)
+    alms = synalm(cls, lmax=lmax, mmax=mmax, new=new, rng=rng, verbose=verbose)
     maps = alm2map(
         alms,
         nside,