@@ -17,11 +17,12 @@ from libc.stdint cimport (uint8_t, uint16_t, uint32_t, uint64_t,
1717 int8_t , int16_t , int32_t , int64_t , intptr_t )
1818from libc .stdlib cimport malloc , free
1919from libc .math cimport sqrt
20- from cpython cimport Py_INCREF , PyComplex_FromDoubles
20+ from cpython cimport Py_INCREF , PyComplex_FromDoubles , PyComplex_RealAsDouble , \
21+ PyComplex_ImagAsDouble , PyInt_AsLong , PyFloat_AsDouble
22+
2123
2224import randomstate
2325from binomial cimport binomial_t
24- from cython_overrides cimport PyFloat_AsDouble , PyInt_AsLong , PyComplex_RealAsDouble , PyComplex_ImagAsDouble
2526from randomstate .entropy import random_entropy
2627
2728np .import_array ()
@@ -1810,7 +1811,7 @@ cdef class RandomState:
18101811 raise ValueError ("method must be either 'bm' or 'zig'" )
18111812 cdef np .ndarray ogamma , orelation , oloc , randoms , v_real , v_imag , rho
18121813 cdef double * randoms_data
1813- cdef double fgamma_r , fgamma_i , frelation_r , frelation_i , frho , f_v_real , f_v_imag , \
1814+ cdef double fgamma_r , fgamma_i , frelation_r , frelation_i , frho , fvar_r , fvar_i , \
18141815 floc_r , floc_i , f_real , f_imag , i_r_scale , r_scale , i_scale , f_rho
18151816 cdef np .npy_intp i , j , n
18161817 cdef np .broadcast it
@@ -1825,19 +1826,19 @@ cdef class RandomState:
18251826 fgamma_r = PyComplex_RealAsDouble (gamma )
18261827 fgamma_i = PyComplex_ImagAsDouble (gamma )
18271828 frelation_r = PyComplex_RealAsDouble (relation )
1828- frelation_i = PyComplex_ImagAsDouble (relation )
1829+ frelation_i = 0.5 * PyComplex_ImagAsDouble (relation )
18291830
1830- f_v_real = fgamma_r + frelation_r
1831- f_v_imag = fgamma_r - frelation_r
1831+ fvar_r = 0.5 * ( fgamma_r + frelation_r )
1832+ fvar_i = 0.5 * ( fgamma_r - frelation_r )
18321833 if fgamma_i != 0 :
18331834 raise ValueError ('Im(gamma) != 0' )
1834- if f_v_imag < 0 :
1835+ if fvar_i < 0 :
18351836 raise ValueError ('Re(gamma - relation) < 0' )
1836- if f_v_real < 0 :
1837+ if fvar_r < 0 :
18371838 raise ValueError ('Re(gamma + relation) < 0' )
18381839 f_rho = 0.0
1839- if f_v_imag > 0 and f_v_real > 0 :
1840- f_rho = frelation_i / sqrt (f_v_imag * f_v_real )
1840+ if fvar_i > 0 and fvar_r > 0 :
1841+ f_rho = frelation_i / sqrt (fvar_i * fvar_r )
18411842 if f_rho > 1.0 or f_rho < - 1.0 :
18421843 raise ValueError ('Im(relation) ** 2 > Re(gamma ** 2 - relation** 2)' )
18431844
@@ -1849,16 +1850,16 @@ cdef class RandomState:
18491850 random_gauss_fill (& self .rng_state , 1 , & f_real )
18501851 random_gauss_fill (& self .rng_state , 1 , & f_imag )
18511852
1852- compute_complex (& f_real , & f_imag , floc_r , floc_i , f_v_real , f_v_imag , f_rho )
1853+ compute_complex (& f_real , & f_imag , floc_r , floc_i , fvar_r , fvar_i , f_rho )
18531854 return PyComplex_FromDoubles (f_real , f_imag )
18541855
18551856 randoms = < np .ndarray > np .empty (size , np .complex128 )
18561857 randoms_data = < double * > np .PyArray_DATA (randoms )
18571858 n = np .PyArray_SIZE (randoms )
18581859
18591860 i_r_scale = sqrt (1 - f_rho * f_rho )
1860- r_scale = sqrt (0.5 * f_v_real )
1861- i_scale = sqrt (0.5 * f_v_imag )
1861+ r_scale = sqrt (fvar_r )
1862+ i_scale = sqrt (fvar_i )
18621863 j = 0
18631864 with self .lock , nogil :
18641865 if method == u'zig' :
@@ -1916,10 +1917,10 @@ cdef class RandomState:
19161917 for i in range (n ):
19171918 floc_r = (< double * > np .PyArray_MultiIter_DATA (it , 1 ))[0 ]
19181919 floc_i = (< double * > np .PyArray_MultiIter_DATA (it , 1 ))[1 ]
1919- f_v_real = (< double * > np .PyArray_MultiIter_DATA (it , 2 ))[0 ]
1920- f_v_imag = (< double * > np .PyArray_MultiIter_DATA (it , 3 ))[0 ]
1920+ fvar_r = (< double * > np .PyArray_MultiIter_DATA (it , 2 ))[0 ]
1921+ fvar_i = (< double * > np .PyArray_MultiIter_DATA (it , 3 ))[0 ]
19211922 f_rho = (< double * > np .PyArray_MultiIter_DATA (it , 4 ))[0 ]
1922- compute_complex (& randoms_data [j ], & randoms_data [j + 1 ], floc_r , floc_i , f_v_real , f_v_imag , f_rho )
1923+ compute_complex (& randoms_data [j ], & randoms_data [j + 1 ], floc_r , floc_i , fvar_r , fvar_i , f_rho )
19231924 j += 2
19241925 np .PyArray_MultiIter_NEXT (it )
19251926
0 commit comments