@@ -17,11 +17,12 @@ from libc.stdint cimport (uint8_t, uint16_t, uint32_t, uint64_t,
17
17
int8_t , int16_t , int32_t , int64_t , intptr_t )
18
18
from libc .stdlib cimport malloc , free
19
19
from 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
+
21
23
22
24
import randomstate
23
25
from binomial cimport binomial_t
24
- from cython_overrides cimport PyFloat_AsDouble , PyInt_AsLong , PyComplex_RealAsDouble , PyComplex_ImagAsDouble
25
26
from randomstate .entropy import random_entropy
26
27
27
28
np .import_array ()
@@ -1810,7 +1811,7 @@ cdef class RandomState:
1810
1811
raise ValueError ("method must be either 'bm' or 'zig'" )
1811
1812
cdef np .ndarray ogamma , orelation , oloc , randoms , v_real , v_imag , rho
1812
1813
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 , \
1814
1815
floc_r , floc_i , f_real , f_imag , i_r_scale , r_scale , i_scale , f_rho
1815
1816
cdef np .npy_intp i , j , n
1816
1817
cdef np .broadcast it
@@ -1825,19 +1826,19 @@ cdef class RandomState:
1825
1826
fgamma_r = PyComplex_RealAsDouble (gamma )
1826
1827
fgamma_i = PyComplex_ImagAsDouble (gamma )
1827
1828
frelation_r = PyComplex_RealAsDouble (relation )
1828
- frelation_i = PyComplex_ImagAsDouble (relation )
1829
+ frelation_i = 0.5 * PyComplex_ImagAsDouble (relation )
1829
1830
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 )
1832
1833
if fgamma_i != 0 :
1833
1834
raise ValueError ('Im(gamma) != 0' )
1834
- if f_v_imag < 0 :
1835
+ if fvar_i < 0 :
1835
1836
raise ValueError ('Re(gamma - relation) < 0' )
1836
- if f_v_real < 0 :
1837
+ if fvar_r < 0 :
1837
1838
raise ValueError ('Re(gamma + relation) < 0' )
1838
1839
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 )
1841
1842
if f_rho > 1.0 or f_rho < - 1.0 :
1842
1843
raise ValueError ('Im(relation) ** 2 > Re(gamma ** 2 - relation** 2)' )
1843
1844
@@ -1849,16 +1850,16 @@ cdef class RandomState:
1849
1850
random_gauss_fill (& self .rng_state , 1 , & f_real )
1850
1851
random_gauss_fill (& self .rng_state , 1 , & f_imag )
1851
1852
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 )
1853
1854
return PyComplex_FromDoubles (f_real , f_imag )
1854
1855
1855
1856
randoms = < np .ndarray > np .empty (size , np .complex128 )
1856
1857
randoms_data = < double * > np .PyArray_DATA (randoms )
1857
1858
n = np .PyArray_SIZE (randoms )
1858
1859
1859
1860
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 )
1862
1863
j = 0
1863
1864
with self .lock , nogil :
1864
1865
if method == u'zig' :
@@ -1916,10 +1917,10 @@ cdef class RandomState:
1916
1917
for i in range (n ):
1917
1918
floc_r = (< double * > np .PyArray_MultiIter_DATA (it , 1 ))[0 ]
1918
1919
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 ]
1921
1922
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 )
1923
1924
j += 2
1924
1925
np .PyArray_MultiIter_NEXT (it )
1925
1926
0 commit comments