Математика - Рассчитайте Перманент Как Можно Быстрее

Задача состоит в том, чтобы написать максимально быстрый код для вычисления перманента матрицы.

Перманент

 
 
 
 
 
 
 
 cat /proc/cpuinfo/|grep flags 
 
 
 
 
 
 
 

-by- n матрица n = ( n 1021509632 ) определяется как

Здесь [[1, -1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, 1, 1, -1], [1, -1, 1, 1, 1, 1, 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1, -1], [-1, -1, 1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, 1, -1, -1, -1, 1, -1, -1], [-1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, 1, -1], [-1, 1, 1, 1, -1, 1, 1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, 1], [1, -1, 1, 1, -1, -1, 1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, -1], [1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1], [1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1], [1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, 1, 1, -1], [-1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, 1, -1, 1, 1], [-1, -1, -1, -1, -1, -1, -1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1], [1, 1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, 1, 1, 1], [-1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, 1], [1, 1, -1, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, 1, -1, 1], [1, 1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1], [1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, -1, -1, -1], [-1, 1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1], [1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 1, 1, 1, -1, 1, 1, -1, 1, -1, 1], [1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1], [-1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, -1]] represents the set of all permutations of 192 .

В качестве примера (из вики):

В этом вопросе матрицы все квадратные и будут иметь только значения [[ 1 -1 1 -1 -1 -1 -1 -1] [-1 -1 1 1 -1 1 1 -1] [ 1 -1 -1 -1 -1 1 1 1] [-1 -1 -1 1 -1 1 1 1] [ 1 -1 -1 1 1 1 1 -1] [-1 1 -1 1 -1 1 1 -1] [ 1 -1 1 -1 1 -1 1 -1] [-1 -1 1 -1 1 1 1 1]] and 0 в них.

Примеры

Вход:

[[-1 -1 -1 -1] [-1 1 -1 -1] [ 1 -1 -1 -1] [ 1 -1 1 -1]]

Постоянный:

-4

Вход:

[[ 1 -1 -1 1] [-1 -1 -1 1] [-1 1 -1 1] [ 1 -1 -1 1]]

Постоянный:

1

Вход:

-1

Постоянный:

[1, n]

Вход:

S_n

Постоянный:

i,j

Задача

Вам следует написать код, который, учитывая a by A матрица, выводит ее постоянную.

Поскольку мне нужно будет протестировать ваш код, было бы полезно, если бы вы могли предоставить мне простой способ ввести матрицу в качестве входных данных для вашего кода, например, путем чтения из стандартного входа.

Имейте в виду, что перманент может быть большим (матрица из единиц — крайний случай).

Счета и ничьи

Я протестирую ваш код на случайных +-1 матрицах возрастающего размера и остановлю работу, когда ваш код в первый раз займет на моем компьютере более 1 минуты. Матрицы оценок будут одинаковыми для всех заявок, чтобы обеспечить справедливость.

Если два человека набирают одинаковое количество очков, то победителем становится тот, кто быстрее всех при этом значении n . If those are within 1 second of each other then it is the one posted first.

Языки и библиотеки

Вы можете использовать любой доступный язык и библиотеки, которые вам нравятся, но не использовать ранее существовавшую функцию для вычисления постоянного значения. Там, где это возможно, было бы хорошо иметь возможность запускать ваш код, поэтому, пожалуйста, включите полное объяснение того, как запустить/скомпилировать ваш код в Linux, если это вообще возможно.`

Эталонные реализации

уже есть вопрос по код-гольфу вопрос с большим количеством кода на разных языках для вычисления перманента для маленьких матриц. Математика и Клен также оба имеют постоянные реализации, если вы можете получить к ним доступ.

Моя машина Тайминги будут работать на моей 64-битной машине. Это стандартная установка Ubuntu с 8 ГБ ОЗУ, восьмиядерным процессором AMD FX-8350 и Radeon HD 4250. Это также означает, что мне нужно иметь возможность запускать ваш код.

Низкоуровневая информация о моей машине

n gives

флаги: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov

pat pse36 clflush mmx fxsr sse sse2 ht системный вызов nx mmxext fxsr_opt

• pdpe1gb rdtscp lm константа_tsc Rep_good nopl нонстоп_tsc extd_apicid aperfmperf pni pclmulqdq монитор ssse3 fma cx16 sse4_1 sse4_2 popcnt
• aes xsave avx f16c lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs xop Skinit wdt lwp fma4 tce nodeid_msr tbm topoext perfctr_core perfctr_nb cpb hw_pstate vmmcall
• bmi1 arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flashbyasid декодирование, пауза, фильтр, pfthreshold
• Я задам тесно связанный дополнительный многоязычный вопрос, который не страдает от большой проблемы с Int, чтобы любители Scala, Nim, Julia, Rust, Bash тоже могли продемонстрировать свои языки. Таблица лидеров n = 33 (45 секунд. 64 секунды для n = 34). Тон Хоспел
• в С++ с г++ 5.4.0. n = 32 (32 секунды).
• Деннис в С с gcc 5.4.0 с использованием флагов gcc Ton Hospel.
• n = 31 (54 секунды). Кристиан Сиверс в Хаскелл

n = 31 (60 секунд).первородный

в

gcc C++ n ≈ 36 (57 секунд в моей системе)

Использует формулу Глинна с кодом Грея для обновлений, если все суммы столбцов четные, в противном случае используется метод Райзера. Резьбовые и векторизованные. Оптимизирован для AVX, поэтому не ждите многого от старых процессоров. Не беспокойтесь о

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 (f '[[1 -1 1 -1 -1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 1 1 -1]

[1 -1 1 1 1 1 1 -1 1 -1 -1 1 1 1 -1 -1 1 1 1 -1]

[-1 -1 1 1 1 -1 -1 -1 -1 1 -1 1 1 1 -1 -1 -1 1 -1 -1]

[-1 -1 -1 1 1 -1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 1 -1]

[-1 1 1 1 -1 1 1 1 -1 -1 -1 1 -1 1 -1 1 1 1 1 1]

[1 -1 1 1 -1 -1 1 -1 1 1 1 1 -1 1 1 -1 1 -1 -1 -1]

[1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 1 -1 -1 -1 1 1 1 -1]

[1 -1 -1 1 -1 1 1 -1 1 1 1 -1 1 -1 1 1 1 -1 1 1]

[1 -1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 -1 1 1 -1 1 1 -1]

[-1 -1 1 -1 1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 1 1]

[-1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1]

[1 1 -1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 1 1 1 1 1 1]

[-1 1 1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 1 -1 1]

[1 1 -1 -1 -1 1 -1 1 -1 -1 -1 -1 1 -1 1 1 -1 1 -1 1]

[1 1 1 1 1 -1 -1 -1 1 1 1 -1 1 -1 1 1 1 -1 1 1]

[1 -1 -1 1 -1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 -1 -1 -1]

[-1 1 1 1 -1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 1 1]

[1 1 -1 -1 1 1 -1 1 1 -1 1 1 1 -1 1 1 -1 1 -1 1]

[1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 -1 -1 1]

[-1 1 1 1 -1 -1 -1 -1 -1 -1 -1 1 1 -1 1 1 -1 1 -1 -1]])
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

for a matrix with only +1's even if your system is fast enough since the signed 128 bit accumulator will overflow. For random matrices you will probably not hit the overflow. For -4 192 внутренние множители начнут переполняться для всех (f '[[ 1 -1 -1 1] [-1 -1 -1 1] [-1 1 -1 1] [ 1 -1 -1 1]]) (f '[[ 1 -1 1 -1 -1 -1 -1 -1] [-1 -1 1 1 -1 1 1 -1] [ 1 -1 -1 -1 -1 1 1 1] [-1 -1 -1 1 -1 1 1 1] [ 1 -1 -1 1 1 1 1 -1] [-1 1 -1 1 -1 1 1 -1] [ 1 -1 1 -1 1 -1 1 -1] [-1 -1 1 -1 1 1 1 1]]) matrix. So only use this program for (define (f ll) (for/sum ((p (permutations (range (length ll))))) (for/product ((l ll)(c p)) (list-ref l c)))) .

Просто укажите элементы матрицы в STDIN, разделенные любыми пробелами.

(for/sum((p(permutations(range(length l)))))(for/product((k l)(c p))(list-ref k c)))

$ time ./matrix-permanent-c 8 30 8395059644858368 real 0m8.582s user 1m8.656s sys 0m0.000s :

1 -1 1 -1 -1 1 1 1 -1 -1 -1 -1 1 1 1 1 -1 1 1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 1 1 1 -1 -1 1 1 1 -1 -1 -1 1 1 1 -1 -1 -1 -1 1 -1 1 1 1 -1 -1 -1 1 -1 -1 -1 -1 -1 1 1 -1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 1 -1 -1 1 1 1 -1 1 1 1 -1 -1 -1 1 -1 1 -1 1 1 1 1 1 1 -1 1 1 -1 -1 1 -1 1 1 1 1 -1 1 1 -1 1 -1 -1 -1 1 -1 -1 1 -1 -1 -1 1 -1 1 1 1 1 -1 -1 -1 1 1 1 -1 1 -1 -1 1 -1 1 1 -1 1 1 1 -1 1 -1 1 1 1 -1 1 1 1 -1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 -1 1 1 -1 1 1 -1 -1 -1 1 -1 1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 1 1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 -1 -1 -1 1 1 -1 -1 1 -1 1 1 -1 1 1 1 1 1 1 -1 1 1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 1 -1 1 1 1 -1 -1 -1 1 -1 1 -1 -1 -1 -1 1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 -1 1 1 1 -1 1 -1 1 1 1 -1 1 1 1 -1 -1 1 -1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 -1 -1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1 1 -1 -1 1 1 -1 -1 1 1 1 1 -1 -1 1 1 -1 1 1 -1 1 1 1 -1 1 1 -1 1 -1 1 1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 1 1 1 -1 -1 -1 -1 -1 -1 -1 1 1 -1 1 1 -1 1 -1 -1 ||answer||

С99, n ≈ 33 (35 секунд)

n

Ввод в настоящее время немного громоздок; он принимается со строками в качестве аргументов командной строки, где каждая запись представлена ​​своим знаком, т. е. + указывает на 1 и - указывает на -1.

Тестовый запуск

unresolved external symbol _fork ||answer||

Питон 2, n ≈ 28

fork

Использует Формула Глинна с Код Грея для обновлений. Добегает до matrix-permanent-c in a minute on my machine. One can surely do better implementing this in a faster language and with better data structures. This doesn't use that the matrix is ±1-valued.

Реализация формулы Райзера очень похожа: суммирование производится по всем векторам коэффициентов 0/1, а не по векторам ±1. Это занимает примерно вдвое больше времени, чем формула Глинна, поскольку складывает все 2^n таких векторов, тогда как Глинн делит это пополам, используя симметрию только для тех, которые начинаются с pypy /path/to/pypy-src/rpython/bin/rpython matrix-permanent.py .

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 ||answer||

Хаскелл, n=31 (54 с)

С большим неоценимым вкладом @Angs: используйте from random import * from time import time from itertools import* def perm(a): # naive method, recurses over submatrices, slow if len(a) == 1: return a[0][0] elif len(a) == 2: return a[0][0]*a[1][1]+a[1][0]*a[0][1] else: tsum = 0 for i in xrange(len(a)): transposed = [zip(*a)[j] for j in xrange(len(a)) if j != i] tsum += a[0][i] * perm(zip(*transposed)[1:]) return tsum def perm_ryser(a): # Ryser's formula, using matrix entries maxn = len(a) n_list = range(1,maxn+1) s_list = chain.from_iterable(combinations(n_list,i) for i in range(maxn+1)) total = 0 for st in s_list: stotal = (-1)**len(st) for i in xrange(maxn): stotal *= sum(a[i][j-1] for j in st) total += stotal return total*((-1)**maxn) def genmatrix(d): mat = [] for x in xrange(d): row = [] for y in xrange(d): row.append([-1,1][randrange(0,2)]) mat.append(row) return mat def main(): for i in xrange(1,24): k = genmatrix(i) print 'Matrix: (%dx%d)'%(i,i) print '\n'.join('['+', '.join(`j`.rjust(2) for j in a)+']' for a in k) print 'Permanent:', t = time() p = perm_ryser(k) print p,'(took',time()-t,'seconds)' if __name__ == '__main__': main() , use short circuit products, look at odd n.

s_list

Мои первые попытки параллелизма в Haskell. В истории изменений можно увидеть множество шагов по оптимизации. Удивительно, но в основном это были очень небольшие изменения. Код основан на формуле из раздела «Формула Баласубраманиана-Бакса/Франклина-Глинна» статьи Википедии на вычисление постоянного.

{1,...,n} computes the permanent. It is called via S который преобразует матрицу всегда допустимым способом, но особенно полезным для матриц, которые мы здесь получаем.

Скомпилировать с n=23 . To run with parallelisation, give it runtime parameters like this: ListConvolve . Ввод осуществляется со стандартного ввода с вложенными списками, разделенными запятыми, в скобках, например Timing as in the last example (newlines allowed everywhere).

Руст + эксприм

Эта простая реализация Ryser с кодом Грея занимает около 65 90 секунд для запуска n=31 на моем ноутбуке. Я предполагаю, что ваша машина доберется туда гораздо раньше, чем за 60-е годы. Я использую extprim 1.1.1 для p[m_] := Last[Fold[Take[ListConvolve[##, {1, -1}, 0], 2^Length[m]]&, Table[If[IntegerQ[Log2[k]], m[[j, Log2[k] + 1]], 0], {j, n}, {k, 0, 2^Length[m] - 1}]]] .

Я никогда не использовал Rust и понятия не имею, что делаю. Никаких опций компилятора, кроме каких-либо use std::env; use std::thread; use std::sync::Arc; use std::sync::mpsc; extern crate extprim; use extprim::i128::i128; static THREADS : i64 = 8; // keep this a power of 2. fn main() { // Read command line args for the matrix, specified like // "++- --- -+-" for [[1, 1, -1], [-1, -1, -1], [-1, 1, -1]]. let mut args = env::args(); args.next(); let mat : Arc<Vec<Vec<i64>>> = Arc::new(args.map( |ss| ss.trim().bytes().map( |cc| if cc == b'+' {1} else {-1}).collect() ).collect()); // Figure how many iterations each thread has to do. let size = 2i64.pow(mat.len() as u32); let slice_size = size / THREADS; // Assumes divisibility. let mut accumulator : i128; if slice_size >= 4 { // permanent() requires 4 divides slice_size let (tx, rx) = mpsc::channel(); // Launch threads. for ii in 0..THREADS { let mat = mat.clone(); let tx = tx.clone(); thread::spawn(move || tx.send(permanent(&mat, ii * slice_size, (ii+1) * slice_size)) ); } // Accumulate results. accumulator = extprim::i128::ZERO; for _ in 0..THREADS { accumulator += rx.recv().unwrap(); } } else { // Small matrix, don't bother threading. accumulator = permanent(&mat, 0, size); } println!("{}", accumulator); } fn permanent(mat: &Vec<Vec<i64>>, start: i64, end: i64) -> i128 { let size = mat.len(); let sentinel = std::i64::MAX / size as i64; let mut bits : Vec<bool> = Vec::with_capacity(size); let mut sums : Vec<i64> = Vec::with_capacity(size); // Initialize gray code bits. let gray_number = start ^ (start / 2); for row in 0..size { bits.push((gray_number >> row) % 2 == 1); sums.push(0); } // Initialize column sums for row in 0..size { if bits[row] { for column in 0..size { sums[column] += mat[row][column]; } } } // Do first two iterations with initial sums let mut total = product(&sums, sentinel); for column in 0..size { sums[column] += mat[0][column]; } bits[0] = true; total -= product(&sums, sentinel); // Do rest of iterations updating gray code bits incrementally let mut gray_bit : usize; let mut idx = start + 2; while idx < end { gray_bit = idx.trailing_zeros() as usize; if bits[gray_bit] { for column in 0..size { sums[column] -= mat[gray_bit][column]; } bits[gray_bit] = false; } else { for column in 0..size { sums[column] += mat[gray_bit][column]; } bits[gray_bit] = true; } total += product(&sums, sentinel); if bits[0] { for column in 0..size { sums[column] -= mat[0][column]; } bits[0] = false; } else { for column in 0..size { sums[column] += mat[0][column]; } bits[0] = true; } total -= product(&sums, sentinel); idx += 2; } return if size % 2 == 0 {total} else {-total}; } #[inline] fn product(sums : &Vec<i64>, sentinel : i64) -> i128 { let mut ret : Option<i128> = None; let mut tally = sums[0]; for ii in 1..sums.len() { if tally.abs() >= sentinel { ret = Some(ret.map_or(i128::new(tally), |n| n * i128::new(tally))); tally = sums[ii]; } else { tally *= sums[ii]; } } if ret.is_none() { return i128::new(tally); } return ret.unwrap() * i128::new(tally); } does. Comments/suggestions/optimizations are appreciated.

Призыв идентичен программе Денниса.

cargo build --release ||answer||

Математика, n ≈ 20

i128

Используя [[1,2],[3,4]] command, a 20x20 matrix requires about 48 seconds on my system. This is not exactly as efficient as the other since it relies on the fact that the permanent can be found as the coefficient of the product of polymomials from each row of the matrix. Efficient polynomial multiplication is performed by creating the coefficient lists and performing convolution using ./<name> +RTS -N . Для этого требуется около О(2^н н²) время при условии, что свертка выполняется с использованием быстрого преобразования Фурье или аналогичного метода, который требует О(н бревно н) время.

Python 2, n = 22 [Справочник]

Это «эталонная» реализация, которой я вчера поделился с Лембиком, она не доходит до ghc -O2 -threaded -fllvm -feager-blackholing -o <name> <name>.hs by a few seconds on his machine, on my machine it does it in about 52 seconds. To achieve these speeds you need to run this through PyPy.

Первая функция вычисляет перманент аналогично тому, как можно вычислить определитель, просматривая каждую подматрицу, пока не останется 2x2, к которому можно применить основное правило. Это невероятно медленно.

Вторая функция реализует функцию Райзера (второе уравнение, указанное в Википедии). Набор pt is essentially the powerset of the numbers p (переменная import Control.Parallel.Strategies import qualified Data.Vector.Unboxed as V import Data.Int type Row = V.Vector Int8 x :: Row -> [Row] -> Integer -> Int -> Integer x p (v:vs) m c = let c' = c - 1 r = if c>0 then parTuple2 rseq rseq else r0 (a,b) = ( x p vs m c' , x (V.zipWith(-) p v) vs (-m) c' ) `using` r in a+b x p _ m _ = prod m p prod :: Integer -> Row -> Integer prod m p = if 0 `V.elem` p then 0 else V.foldl' (\a b->a*fromIntegral b) m p p, pt :: [Row] -> Integer p (v:vs) = x (foldl (V.zipWith (+)) v vs) (map (V.map (2*)) vs) 1 11 `div` 2^(length vs) p [] = 1 -- handle 0x0 matrices too :-) pt (v:vs) | even (length vs) = p ((V.map (2*) v) : vs ) `div` 2 pt mat = p mat main = getContents >>= print . pt . map V.fromList . read in the code).

Vector ||answer||

RPython 5.4.1, n ≈ 32 (37 секунд)

from operator import mul def fast_ryser_perm(M): n=len(M) row_comb = [0] * n total = 0 old_grey = 0 sign = +1 binary_power_dict = {2**i:i for i in range(n)} num_loops = 2**n for bin_index in range(1, num_loops) + [0]: total += sign * reduce(mul, row_comb) new_grey = bin_index^(bin_index/2) grey_diff = old_grey ^ new_grey grey_diff_index = binary_power_dict[grey_diff] new_vector = M[grey_diff_index] direction = cmp(old_grey, new_grey) for i in range(n): row_comb[i] += new_vector[i] * direction sign = -sign old_grey = new_grey return total * (-1)**n

Чтобы скомпилировать, скачать самый последний исходный код PyPy и выполните следующее:

+1

Полученный исполняемый файл будет называться n=23 or similiar in the current working directory.

Начиная с PyPy 5.0, примитивы потоков RPython стали намного менее примитивными, чем раньше. Вновь созданным потокам требуется GIL, который более или менее бесполезен для параллельных вычислений. я использовал from operator import mul def fast_glynn_perm(M): row_comb = [sum(c) for c in zip(*M)] n=len(M) total = 0 old_grey = 0 sign = +1 binary_power_dict = {2**i:i for i in range(n)} num_loops = 2**(n-1) for bin_index in xrange(1, num_loops + 1): total += sign * reduce(mul, row_comb) new_grey = bin_index^(bin_index/2) grey_diff = old_grey ^ new_grey grey_diff_index = binary_power_dict[grey_diff] new_vector = M[grey_diff_index] direction = 2 * cmp(old_grey,new_grey) for i in range(n): row_comb[i] += new_vector[i] * direction sign = -sign old_grey = new_grey return total/num_loops instead, so it may not work as expected on Windows, although I haven't tested fails to compile ( $ gcc -Wall -std=c99 -march=native -Ofast -fopenmp -fwrapv -o permanent permanent.c $ ./permanent +--+ ---+ -+-+ +--+ -4 $ ./permanent ---- -+-- +--- +-+- 0 $ ./permanent +-+----- --++-++- +----+++ ---+-+++ +--++++- -+-+-++- +-+-+-+- --+-++++ 192 $ ./permanent +-+--+++----++++-++- +-+++++-+--+++--+++- --+++----+-+++---+-- ---++-++++++------+- -+++-+++---+-+-+++++ +-++--+-++++-++-+--- +--+---+-++++---+++- +--+-++-+++-+-+++-++ +-----+++-----++-++- --+-+-++-+-++++++-++ -------+----++++---- ++---++--+-++-++++++ -++-----++++-----+-+ ++---+-+----+-++-+-+ +++++---+++-+-+++-++ +--+----+--++-+----- -+++-++--+++--++--++ ++--++-++-+++-++-+-+ +++---+--++---+----+ -+++-------++-++-+-- 1021509632 $ time ./permanent +++++++++++++++++++++++++++++++{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,} # 31 8222838654177922817725562880000000 real 0m8.365s user 1m6.504s sys 0m0.000s $ time ./permanent ++++++++++++++++++++++++++++++++{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,} # 32 263130836933693530167218012160000000 real 0m17.013s user 2m15.226s sys 0m0.001s $ time ./permanent +++++++++++++++++++++++++++++++++{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,} # 33 8683317618811886495518194401280000000 real 0m34.592s user 4m35.354s sys 0m0.001s ).

Исполняемый файл принимает до двух параметров командной строки. Первый — количество потоков, второй необязательный параметр — #include <stdint.h> #include <stdio.h> #define CHUNK_SIZE 12 #define NUM_THREADS 8 #define popcnt __builtin_popcountll #define BILLION (1000 * 1000 * 1000) #define UPDATE_ROW_PPROD() \ update_row_pprod(row_pprod, row, rows, row_sums, mask, mask_popcnt) typedef __int128 int128_t; static inline int64_t update_row_pprod ( int64_t* row_pprod, int64_t row, int64_t* rows, int64_t* row_sums, int64_t mask, int64_t mask_popcnt ) { int64_t temp = 2 * popcnt(rows[row] & mask) - mask_popcnt; row_pprod[0] *= temp; temp -= 1; row_pprod[1] *= temp; temp -= row_sums[row]; row_pprod[2] *= temp; temp += 1; row_pprod[3] *= temp; return row + 1; } int main(int argc, char* argv[]) { int64_t size = argc - 1, rows[argc - 1]; int64_t row_sums[argc - 1]; int128_t permanent = 0, sign = size & 1 ? -1 : 1; if (argc == 2) { printf("%d\n", argv[1][0] == '-' ? -1 : 1); return 0; } for (int64_t row = 0; row < size; row++) { char positive = argv[row + 1][0] == '+' ? '-' : '+'; sign *= ',' - positive; rows[row] = row_sums[row] = 0; for (char* p = &argv[row + 1][1]; *p; p++) { rows[row] <<= 1; rows[row] |= *p == positive; row_sums[row] += *p == positive; } row_sums[row] = 2 * row_sums[row] - size; } #pragma omp parallel for reduction(+:permanent) num_threads(NUM_THREADS) for (int64_t mask = 1; mask < 1LL << (size - 1); mask += 2) { int64_t mask_popcnt = popcnt(mask); int64_t row = 0; int128_t row_prod = 1 - 2 * (mask_popcnt & 1); int128_t row_prod_high = -row_prod; int128_t row_prod_inv = row_prod; int128_t row_prod_inv_high = -row_prod; for (int64_t chunk = 0; chunk < size / CHUNK_SIZE; chunk++) { int64_t row_pprod[4] = {1, 1, 1, 1}; for (int64_t i = 0; i < CHUNK_SIZE; i++) row = UPDATE_ROW_PPROD(); row_prod *= row_pprod[0], row_prod_high *= row_pprod[1]; row_prod_inv *= row_pprod[3], row_prod_inv_high *= row_pprod[2]; } int64_t row_pprod[4] = {1, 1, 1, 1}; while (row < size) row = UPDATE_ROW_PPROD(); row_prod *= row_pprod[0], row_prod_high *= row_pprod[1]; row_prod_inv *= row_pprod[3], row_prod_inv_high *= row_pprod[2]; permanent += row_prod + row_prod_high + row_prod_inv + row_prod_inv_high; } permanent *= sign; if (permanent < 0) printf("-"), permanent *= -1; int32_t output[5], print = 0; output[0] = permanent % BILLION, permanent /= BILLION; output[1] = permanent % BILLION, permanent /= BILLION; output[2] = permanent % BILLION, permanent /= BILLION; output[3] = permanent % BILLION, permanent /= BILLION; output[4] = permanent % BILLION; if (output[4]) printf("%u", output[4]), print = 1; if (print) printf("%09u", output[3]); else if (output[3]) printf("%u", output[3]), print = 1; if (print) printf("%09u", output[2]); else if (output[2]) printf("%u", output[2]), print = 1; if (print) printf("%09u", output[1]); else if (output[1]) printf("%u", output[1]), print = 1; if (print) printf("%09u\n", output[0]); else printf("%u\n", output[0]); } . 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:

/* Compile using something like: g++ -Wall -O3 -march=native -fstrict-aliasing -std=c++11 -pthread -s permanent.cpp -o permanent */ #include <iostream> #include <iomanip> #include <cstdlib> #include <cstdint> #include <climits> #include <array> #include <vector> #include <thread> #include <future> #include <ctgmath> #include <immintrin.h> using namespace std; bool const DEBUG = false; int const CACHE = 64; using Index = int_fast32_t; Index glynn; // Number of elements in our vectors Index const POW = 3; Index const ELEMS = 1 << POW; // Over how many floats we distribute each row Index const WIDTH = 9; // Number of bits in the fraction part of a floating point number int const FLOAT_MANTISSA = 23; // Type to use for the first add/multiply phase using Sum = float; using SumN = __restrict__ Sum __attribute__((vector_size(ELEMS*sizeof(Sum)))); // Type to convert to between the first and second phase using ProdN = __restrict__ int32_t __attribute__((vector_size(ELEMS*sizeof(int32_t)))); // Type to use for the third and last multiply phase. // Also used for the final accumulator using Value = __int128; using UValue = unsigned __int128; // Wrap Value so C++ doesn't really see it and we can put it in vectors etc. // Needed since C++ doesn't fully support __int128 struct Number { Number& operator+=(Number const& right) { value += right.value; return *this; } // Output the value void print(ostream& os, bool dbl = false) const; friend ostream& operator<<(ostream& os, Number const& number) { number.print(os); return os; } Value value; }; using ms = chrono::milliseconds; auto nr_threads = thread::hardware_concurrency(); vector<Sum> input; // Allocate cache aligned datastructures template<typename T> T* alloc(size_t n) { T* mem = static_cast<T*>(aligned_alloc(CACHE, sizeof(T) * n)); if (mem == nullptr) throw(bad_alloc()); return mem; } // Work assigned to thread k of nr_threads threads Number permanent_part(Index n, Index k, SumN** more) { uint64_t loops = (UINT64_C(1) << n) / nr_threads; if (glynn) loops /= 2; Index l = loops < ELEMS ? loops : ELEMS; loops /= l; auto from = loops * k; auto to = loops * (k+1); if (DEBUG) cout << "From=" << from << "\n"; uint64_t old_gray = from ^ from/2; uint64_t bit = 1; bool bits = (to-from) & 1; Index nn = (n+WIDTH-1)/WIDTH; Index ww = nn * WIDTH; auto column = alloc<SumN>(ww); for (Index i=0; i<n; ++i) for (Index j=0; j<ELEMS; ++j) column[i][j] = 0; for (Index i=n; i<ww; ++i) for (Index j=0; j<ELEMS; ++j) column[i][j] = 1; Index b; if (glynn) { b = n > POW+1 ? n - POW - 1: 0; auto c = n-1-b; for (Index k=0; k<l; k++) { Index gray = k ^ k/2; for (Index j=0; j< c; ++j) if (gray & 1 << j) for (Index i=0; i<n; ++i) column[i][k] -= input[(b+j)*n+i]; else for (Index i=0; i<n; ++i) column[i][k] += input[(b+j)*n+i]; } for (Index i=0; i<n; ++i) for (Index k=0; k<l; k++) column[i][k] += input[n*(n-1)+i]; for (Index k=1; k<l; k+=2) column[0][k] = -column[0][k]; for (Index i=0; i<b; ++i, bit <<= 1) { if (old_gray & bit) { bits = bits ^ 1; for (Index j=0; j<ww; ++j) column[j] -= more[i][j]; } else { for (Index j=0; j<ww; ++j) column[j] += more[i][j]; } } for (Index i=0; i<n; ++i) for (Index k=0; k<l; k++) column[i][k] /= 2; } else { b = n > POW ? n - POW : 0; auto c = n-b; for (Index k=0; k<l; k++) { Index gray = k ^ k/2; for (Index j=0; j<c; ++j) if (gray & 1 << j) for (Index i=0; i<n; ++i) column[i][k] -= input[(b+j)*n+i]; } for (Index k=1; k<l; k+=2) column[0][k] = -column[0][k]; for (Index i=0; i<b; ++i, bit <<= 1) { if (old_gray & bit) { bits = bits ^ 1; for (Index j=0; j<ww; ++j) column[j] -= more[i][j]; } } } if (DEBUG) { for (Index i=0; i<ww; ++i) { cout << "Column[" << i << "]="; for (Index j=0; j<ELEMS; ++j) cout << " " << column[i][j]; cout << "\n"; } } --more; old_gray = (from ^ from/2) | UINT64_C(1) << b; Value total = 0; SumN accu[WIDTH]; for (auto p=from; p<to; ++p) { uint64_t new_gray = p ^ p/2; uint64_t bit = old_gray ^ new_gray; Index i = __builtin_ffsl(bit); auto diff = more[i]; auto c = column; if (new_gray > old_gray) { // Phase 1 add/multiply. // Uses floats until just before loss of precision for (Index i=0; i<WIDTH; ++i) accu[i] = *c++ -= *diff++; for (Index j=1; j < nn; ++j) for (Index i=0; i<WIDTH; ++i) accu[i] *= *c++ -= *diff++; } else { // Phase 1 add/multiply. // Uses floats until just before loss of precision for (Index i=0; i<WIDTH; ++i) accu[i] = *c++ += *diff++; for (Index j=1; j < nn; ++j) for (Index i=0; i<WIDTH; ++i) accu[i] *= *c++ += *diff++; } if (DEBUG) { cout << "p=" << p << "\n"; for (Index i=0; i<ww; ++i) { cout << "Column[" << i << "]="; for (Index j=0; j<ELEMS; ++j) cout << " " << column[i][j]; cout << "\n"; } } // Convert floats to int32_t ProdN prod32[WIDTH] __attribute__((aligned (32))); for (Index i=0; i<WIDTH; ++i) // Unfortunately gcc doesn't recognize the static_cast<int32_t> // as a vector pattern, so force it with an intrinsic #ifdef __AVX__ //prod32[i] = static_cast<ProdN>(accu[i]); reinterpret_cast<__m256i&>(prod32[i]) = _mm256_cvttps_epi32(accu[i]); #else // __AVX__ for (Index j=0; j<ELEMS; ++j) prod32[i][j] = static_cast<int32_t>(accu[i][j]); #endif // __AVX__ // Phase 2 multiply. Uses int64_t until just before overflow int64_t prod64[3][ELEMS]; for (Index i=0; i<3; ++i) { for (Index j=0; j<ELEMS; ++j) prod64[i][j] = static_cast<int64_t>(prod32[i][j]) * prod32[i+3][j] * prod32[i+6][j]; } // Phase 3 multiply. Collect into __int128. For large matrices this will // actually overflow but that's ok as long as all 128 low bits are // correct. Terms will cancel and the final sum can fit into 128 bits // (This will start to fail at n=35 for the all 1 matrix) // Strictly speaking this needs the -fwrapv gcc option for (Index j=0; j<ELEMS; ++j) { auto value = static_cast<Value>(prod64[0][j]) * prod64[1][j] * prod64[2][j]; if (DEBUG) cout << "value[" << j << "]=" << static_cast<double>(value) << "\n"; total += value; } total = -total; old_gray = new_gray; } return bits ? Number{-total} : Number{total}; } // Prepare datastructures, Assign work to threads Number permanent(Index n) { Index nn = (n+WIDTH-1)/WIDTH; Index ww = nn*WIDTH; Index rows = n > (POW+glynn) ? n-POW-glynn : 0; auto data = alloc<SumN>(ww*(rows+1)); auto pointers = alloc<SumN *>(rows+1); auto more = &pointers[0]; for (Index i=0; i<rows; ++i) more[i] = &data[ww*i]; more[rows] = &data[ww*rows]; for (Index j=0; j<ww; ++j) for (Index i=0; i<ELEMS; ++i) more[rows][j][i] = 0; Index loops = n >= POW+glynn ? ELEMS : 1 << (n-glynn); auto a = &input[0]; for (Index r=0; r<rows; ++r) { for (Index j=0; j<n; ++j) { for (Index i=0; i<loops; ++i) more[r][j][i] = j == 0 && i %2 ? -*a : *a; for (Index i=loops; i<ELEMS; ++i) more[r][j][i] = 0; ++a; } for (Index j=n; j<ww; ++j) for (Index i=0; i<ELEMS; ++i) more[r][j][i] = 0; } if (DEBUG) for (Index r=0; r<=rows; ++r) for (Index j=0; j<ww; ++j) { cout << "more[" << r << "][" << j << "]="; for (Index i=0; i<ELEMS; ++i) cout << " " << more[r][j][i]; cout << "\n"; } // Send work to threads... vector<future<Number>> results; for (auto i=1U; i < nr_threads; ++i) results.emplace_back(async(DEBUG ? launch::deferred: launch::async, permanent_part, n, i, more)); // And collect results auto r = permanent_part(n, 0, more); for (auto& result: results) r += result.get(); free(data); free(pointers); // For glynn we should double the result, but we will only do this during // the final print. This allows n=34 for an all 1 matrix to work // if (glynn) r *= 2; return r; } // Print 128 bit number void Number::print(ostream& os, bool dbl) const { const UValue BILLION = 1000000000; UValue val; if (value < 0) { os << "-"; val = -value; } else val = value; if (dbl) val *= 2; uint32_t output[5]; for (int i=0; i<5; ++i) { output[i] = val % BILLION; val /= BILLION; } bool print = false; for (int i=4; i>=0; --i) { if (print) { os << setfill('0') << setw(9) << output[i]; } else if (output[i] || i == 0) { print = true; os << output[i]; } } } // Read matrix, check for sanity void my_main() { Sum a; while (cin >> a) input.push_back(a); size_t n = sqrt(input.size()); if (input.size() != n*n) throw(logic_error("Read " + to_string(input.size()) + " elements which does not make a square matrix")); vector<double> columns_pos(n, 0); vector<double> columns_neg(n, 0); Sum *p = &input[0]; for (size_t i=0; i<n; ++i) for (size_t j=0; j<n; ++j, ++p) { if (*p >= 0) columns_pos[j] += *p; else columns_neg[j] -= *p; } std::array<double,WIDTH> prod; prod.fill(1); int32_t odd = 0; for (size_t j=0; j<n; ++j) { prod[j%WIDTH] *= max(columns_pos[j], columns_neg[j]); auto sum = static_cast<int32_t>(columns_pos[j] - columns_neg[j]); odd |= sum; } glynn = (odd & 1) ^ 1; for (Index i=0; i<WIDTH; ++i) // A float has an implicit 1. in front of the fraction so it can // represent 1 bit more than the mantissa size. And 1 << (mantissa+1) // itself is in fact representable if (prod[i] && log2(prod[i]) > FLOAT_MANTISSA+1) throw(range_error("Values in matrix are too large. A subproduct reaches " + to_string(prod[i]) + " which doesn't fit in a float without loss of precision")); for (Index i=0; i<3; ++i) { auto prod3 = prod[i] * prod[i+3] * prod[i+6]; if (log2(prod3) >= CHAR_BIT*sizeof(int64_t)-1) throw(range_error("Values in matrix are too large. A subproduct reaches " + to_string(prod3) + " which doesn't fit in an int64")); } nr_threads = pow(2, ceil(log2(static_cast<float>(nr_threads)))); uint64_t loops = UINT64_C(1) << n; if (glynn) loops /= 2; if (nr_threads * ELEMS > loops) nr_threads = max(loops / ELEMS, UINT64_C(1)); // if (DEBUG) nr_threads = 1; cout << n << " x " << n << " matrix, method " << (glynn ? "Glynn" : "Ryser") << ", " << nr_threads << " threads" << endl; // Go for the actual calculation auto start = chrono::steady_clock::now(); auto perm = permanent(n); auto end = chrono::steady_clock::now(); auto elapsed = chrono::duration_cast<ms>(end-start).count(); cout << "Permanent="; perm.print(cout, glynn); cout << " (" << elapsed / 1000. << " s)" << endl; } // Wrapper to print any exceptions int main() { try { my_main(); } catch(exception& e) { cerr << "Error: " << e.what() << endl; exit(EXIT_FAILURE); } exit(EXIT_SUCCESS); }

---

Пример использования

permanent.cpp

---

Метод

Я использовал Формула Баласубраманиана-Бакса/Франклина-Глинна, со сложностью выполнения О(2^нп). Однако вместо повторения δ в порядке кода Грея я вместо этого заменил умножение векторных строк одной операцией xor (отображение (1, -1) → (0, 1)). Векторную сумму также можно найти за одну операцию, взяв n минус удвоенное значение popcount.

Ракетка 84 байта

Следующая простая функция работает для матриц меньшего размера, но зависает на моей машине для матриц большего размера:

permanent 1 2 3 4 ^D

Не в гольфе:

n<=36

Код можно легко изменить для неодинакового количества строк и столбцов.

Тестирование:

1/-1

Выход:

n>=37

Как я уже упоминал выше, при тестировании он зависает следующим образом:

n>=35