<p>这个答案基于之前发布的@romeric代码。我更正了代码并简化了它,并添加了<code>cdivision</code>编译器指令。在</p>
<pre><code>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def permfunc_modified_2(np.ndarray [double, ndim =2, mode='c'] M):
cdef:
int n = M.shape[0], s=1, i, j
int *f = <int*>malloc(n*sizeof(int))
double *d = <double*>malloc(n*sizeof(double))
double *v = <double*>malloc(n*sizeof(double))
double p = 1, prod
for i in range(n):
v[i] = 0.
for j in range(n):
v[i] += M[j,i]
p *= v[i]
f[i] = i
d[i] = 1
j = 0
while (j < n-1):
prod = 1.
for i in range(n):
v[i] -= 2.*d[j]*M[j, i]
prod *= v[i]
d[j] = -d[j]
s = -s
p += s*prod
f[0] = 0
f[j] = f[j+1]
f[j+1] = j+1
j = f[0]
free(d)
free(f)
free(v)
return p/pow(2.,(n-1))
</code></pre>
<p>@romeric的原始代码没有初始化<code>v</code>,因此有时会得到不同的结果。另外,我分别组合了<code>while</code>之前的两个循环和{<cd3>}内部的两个循环。在</p>
<p>最后,比较</p>
^{pr2}$