ENH: Use low level access in complex normal

Switch to low level interface to improve performance of generated
complex normal

Author: bashtage

randomstate/randomstate.pyx

@@ -1782,9 +1782,12 @@ cdef class RandomState:
 
         >>> s = np.random.complex_normal(size=1000)
         """
-        cdef np.ndarray ogamma, orelation, oloc
+        # TODO: What about loc??
+        cdef np.ndarray ogamma, orelation, oloc, randoms
+        cdef double *randoms_data
         cdef double fgamma_r, fgamma_i, frelation_r, frelation_i, frho, f_v_real , f_v_imag, \
-            floc_r, floc_i, f_real, f_imag
+            floc_r, floc_i, f_real, f_imag, i_r_scale, r_scale, i_scale
+        cdef np.npy_intp i, j, n,
 
         oloc = <np.ndarray>np.PyArray_FROM_OTF(loc, np.NPY_COMPLEX128, np.NPY_ALIGNED)
         ogamma = <np.ndarray>np.PyArray_FROM_OTF(gamma, np.NPY_COMPLEX128, np.NPY_ALIGNED)
@@ -1813,30 +1816,32 @@ cdef class RandomState:
                 raise ValueError('Im(relation) ** 2 > Re(gamma ** 2 - relation** 2)')
 
             if size is None:
-                f_real, f_imag = self.standard_normal(size=2, method=method)
-                
+                random_gauss_zig_double_fill(&self.rng_state, 1, &f_real)
+                random_gauss_zig_double_fill(&self.rng_state, 1, &f_imag)
                 f_imag *= sqrt(1 - f_rho * f_rho)
                 f_imag += f_rho * f_real
                 f_real *= sqrt(0.5 * f_v_real)
                 f_imag *= sqrt(0.5 * f_v_imag)
 
                 return PyComplex_FromDoubles(f_real, f_imag)
 
-            if np.PyArray_IsAnyScalar(size):
-                size = (size,)
-            else:
-                size = tuple(size)
-
-            norms = self.standard_normal(size=size + (2,), method=method)
-            real = norms[...,0]
-            imag = norms[...,1]
+            randoms = <np.ndarray>np.empty(size, np.complex128)
+            randoms_data = <double *>np.PyArray_DATA(randoms)
+            n = np.PyArray_SIZE(randoms)
 
-            imag *= sqrt(1 - f_rho * f_rho)
-            imag += f_rho * real
-            real *= sqrt(0.5 * f_v_real)
-            imag *= sqrt(0.5 * f_v_imag)
+            i_r_scale = sqrt(1 - f_rho * f_rho)
+            r_scale = sqrt(0.5 * f_v_real)
+            i_scale = sqrt(0.5 * f_v_imag)
+            j = 0
+            with self.lock, nogil:
+                for i in range(n):
+                    random_gauss_zig_double_fill(&self.rng_state, 1, &f_real)
+                    random_gauss_zig_double_fill(&self.rng_state, 1, &f_imag)
+                    randoms_data[j+1] = i_scale * (f_rho * f_real + i_r_scale * f_imag)
+                    randoms_data[j] = r_scale * f_real
+                    j += 2
 
-            return floc_r + real + (floc_i + imag) * (0+1.0j)
+            return randoms
 
         gpc = ogamma + orelation
         gmc = ogamma - orelation

setup.py

@@ -240,7 +240,7 @@ def cythonize(e, *args, **kwargs):
         with open(output_file_name, 'w') as output_file:
             output_file.write(template.substitute())
 
-ext_modules = cythonize(extensions, force=not DEVELOP)
+ext_modules = cythonize(extensions, force=not DEVELOP, annotate=True)
 
 classifiers = ['Development Status :: 5 - Production/Stable',
                'Environment :: Console',