# RPython 5.4.1, n ≈ 32 (37 seconds)
<!-- language-all: lang-python -->
from rpython.rlib.rtime import time
from rpython.rlib.rarithmetic import r_int, r_uint
from rpython.rlib.rrandom import Random
from rpython.rlib.rposix import pipe, close, read, write, fork, waitpid
from rpython.rlib.rbigint import rbigint
from math import log, ceil
from struct import pack
bitsize = len(pack('l', 1)) * 8 - 1
bitcounts = bytearray([0])
for i in range(16):
b = bytearray([j+1 for j in bitcounts])
bitcounts += b
def bitcount(n):
bits = 0
while n:
bits += bitcounts[n & 65535]
n >>= 16
return bits
def main(argv):
if len(argv) < 2:
write(2, 'Usage: %s NUM_THREADS [N]'%argv[0])
return 1
threads = int(argv[1])
if len(argv) > 2:
n = int(argv[2])
rnd = Random(r_uint(time()*1000))
m = []
for i in range(n):
row = []
for j in range(n):
row.append(1 - r_int(rnd.genrand32() & 2))
m.append(row)
else:
m = []
strm = ""
while True:
buf = read(0, 4096)
if len(buf) == 0:
break
strm += buf
rows = strm.split("\n")
for row in rows:
r = []
for val in row.split(' '):
r.append(int(val))
m.append(r)
n = len(m)
a = []
for row in m:
val = 0
for v in row:
val = (val << 1) | -(v >> 1)
a.append(val)
batches = int(ceil(n * log(n) / (bitsize * log(2))))
pids = []
handles = []
total = rbigint.fromint(0)
for i in range(threads):
r, w = pipe()
pid = fork()
if pid:
close(w)
pids.append(pid)
handles.append(r)
else:
close(r)
total = run(n, a, i, threads, batches)
write(w, total.str())
close(w)
return 0
for pid in pids:
waitpid(pid, 0)
for handle in handles:
strval = read(handle, 256)
total = total.add(rbigint.fromdecimalstr(strval))
close(handle)
print total.rshift(n-1).str()
return 0
def run(n, a, mynum, threads, batches):
start = (1 << n-1) * mynum / threads
end = (1 << n-1) * (mynum+1) / threads
dtotal = rbigint.fromint(0)
for delta in range(start, end):
pdelta = rbigint.fromint(1 - ((bitcount(delta) & 1) << 1))
for i in range(batches):
pbatch = 1
for j in range(i, n, batches):
pbatch *= n - (bitcount(delta ^ a[j]) << 1)
pdelta = pdelta.int_mul(pbatch)
dtotal = dtotal.add(pdelta)
return dtotal
def target(*args):
return main
To compile, [download](http://pypy.org/download.html) the most recent PyPy source, and execute the following:
pypy /path/to/pypy-src/rpython/bin/rpython matrix-permanent.py
The resulting executable will be named `matrix-permanent-c` or similiar in the current working directory.
As of PyPy 5.0, RPython's threading primitives are a lot less primitive than they used to be. Newly spawned threads require the GIL, which is more or less useless for parallel computations. I've used `fork` instead, so it may not work as expected on Windows, <s>although I haven't tested</s> fails to compile (`unresolved external symbol _fork`).
The executable accepts up to two command line parameters. The first is the number of threads, the second optional parameter is `n`. If it is provided, a random matrix will be generated, otherwise it will be read from stdin. Each row must be newline separated (without a trailing newline), and each value space separated. The third example input would be given as:
1 -1 1 -1 -1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 1 1 -1
1 -1 1 1 1 1 1 -1 1 -1 -1 1 1 1 -1 -1 1 1 1 -1
-1 -1 1 1 1 -1 -1 -1 -1 1 -1 1 1 1 -1 -1 -1 1 -1 -1
-1 -1 -1 1 1 -1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 1 -1
-1 1 1 1 -1 1 1 1 -1 -1 -1 1 -1 1 -1 1 1 1 1 1
1 -1 1 1 -1 -1 1 -1 1 1 1 1 -1 1 1 -1 1 -1 -1 -1
1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 1 -1 -1 -1 1 1 1 -1
1 -1 -1 1 -1 1 1 -1 1 1 1 -1 1 -1 1 1 1 -1 1 1
1 -1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 -1 1 1 -1 1 1 -1
-1 -1 1 -1 1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 1 1
-1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1
1 1 -1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 1 1 1 1 1 1
-1 1 1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 1 -1 1
1 1 -1 -1 -1 1 -1 1 -1 -1 -1 -1 1 -1 1 1 -1 1 -1 1
1 1 1 1 1 -1 -1 -1 1 1 1 -1 1 -1 1 1 1 -1 1 1
1 -1 -1 1 -1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 -1 -1 -1
-1 1 1 1 -1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 1 1
1 1 -1 -1 1 1 -1 1 1 -1 1 1 1 -1 1 1 -1 1 -1 1
1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 -1 -1 1
-1 1 1 1 -1 -1 -1 -1 -1 -1 -1 1 1 -1 1 1 -1 1 -1 -1
---
**Sample Usage**
$ time ./matrix-permanent-c 8 30
8395059644858368
real 0m8.582s
user 1m8.656s
sys 0m0.000s
---
**Method**
I've used the [Balasubramanian-Bax/Franklin-Glynn formula](https://en.wikipedia.org/wiki/Computing_the_permanent#Balasubramanian-Bax.2FFranklin-Glynn_formula), with a runtime complexity of _O(2<sup>n</sup>n)_. However, instead of iterating the _δ_ in grey code order, I've instead replaced vector-row multiplication with a single xor operation (mapping (1, -1) → (0, 1)). The vector sum can likewise be found in a single operation, by taking n minus twice the popcount.