Перманент матрицы

[maxal]

wrest в сообщении #1600055 писал(а):

Скажите пож-ста, 40 секунд Sage и 8 секунд PARI/GP -- это один и тот же комп? В коде Sage как я понимаю включена мемоизация (хотя там мне кажется и хитов нету), так что рекурсия против прямого цикла вроде не проблема. Или разные компы?

Комп один и тот же. Просто Sage проигрывает в производительности PARI/GP, ничего удивительного.

[Null]

Geen в сообщении #1600042 писал(а):

На самом деле меньше

Длинна чисел $O(n\log(n))$

[vicvolf]

maxal в сообщении #1600061 писал(а):

Комп один и тот же. Просто Sage проигрывает в производительности PARI/GP, ничего удивительного.

Конечно время счета одного и того же алгоритма зависит от программы и производительности компа. Сложность алгоритма зависит от количества длинных и коротких операций или от средней длины числа. В данном случае наилучший показатель:

-- 06.07.2023, 08:09 --

Null в сообщении #1600064 писал(а):

$O(n\log(n))$

[kthxbye]

Ну что ж, наверное стоит добавить результаты в OEIS. Предлагаю следующие варианты программ на PARI/GP:

A204262

upto(n)=my(v1); v1=vector(n+1, i, i--; i!*x^i); for(i=1, n, for(j=i, n, my(A=intformal((j-i)^2*v1[j])); v1[j+1] = A + subst(v1[j+1] - A, x, i))); v1

A204264

a(n)=if(n==0, 1, my(v1); v1=vector(n+1, i, i--; i!*x^i); for(i=1, n, for(j=i, n, my(A=intformal((j-i)^2*v1[j])); v1[j+1] = A + subst(v1[j+1] - A, x, n - i + 1))); v1[n+1])

Добавил бы сам, но у меня все 3 слота правок заняты. Желательно также добавить ссылку на эту тему. Если кто-нибудь откроет правку к A204262, то я смогу добавить туда свою формулу (если успею до аппрува). Результат участника Null тоже желательно оформить формулой (при этом я не уверен как это правильно сделать).

[wrest]

kthxbye
Сделайте так, чтобы код не зависел от того, определён ранее x или нет.
Иначе может быть так:

? x=6;
? upto(3)
  ***   at top-level: upto(3)
  ***                 ^-------
  ***   in function upto: ...1[j+1]=A+subst(v1[j+1]-A,x,i)));v1
  ***                                                 ^---------
  ***   incorrect type in evaluator [variable name expected] (t_INT).
  ***   Break loop: type 'break' to go back to GP prompt
break>

И это в лучшем случае. В худшем - может неверно считать.

Например, вместо upto(n)=my(v1); сделайте upto(n)=my(v1, x='x);

Да, и по правилам хорошего тона, вместо A=intformal((j-i)^2*v1[j]) лучше писать A=intformal((j-i)^2*v1[j],x)

P.S. Но в целом, хотел вам сказать что я весьма впечатлён. Хороший код, ничего лишнего. Прогресс налицо!
P.P.S. Меня смущает intformal() но я не знаю, ускорит ли тут что-то "прямой" код.

[wrest]

Ещё процентов на 20 побыстрее:

{upto(n)=
my(v,x='x,A);
v1=vector(n+1,i,i--;i!*x^i);
for(i=1,n,
  for(j=i,n,
    if(i==j,A=0,A=intformal((j-i)^2*v[j],x));
    if(j>(2*i),v[j+1]=A,v1[j+1]=A+subst(v[j+1]-A,x,i))
  )
 );
return(v)}

[Null]

wrest в сообщении #1600404 писал(а):

if(j>(2*i),v[j+1]=A,v1[j+1]=A+subst(v[j+1]-A,x,i))

Да. Матрица с большим квадратом из 0 имеет нулевой перманент.

[wrest]

Нельзя ли снизить сложность до $O(n^2)$?

[kthxbye]

Задал вопрос на MathOverflow и получил ответ от Теренса Тао. Оказывается, моя формула непосредственно связана с формулой участника Null.

[flamehowk]

Geen в сообщении #1600026 писал(а):

Алгоритм Null на питоне:

Используется синтаксис Python

import datetime
t = datetime.datetime.now();
dt = datetime.timedelta(seconds=1);
sys.set_int_max_str_digits(10000);
N = 1000;
N += 1;
res = [0]*N;
tt = [0]*N;
Fl = [[0]*N];
Fl[0][-1] = 1;
for n in range(1,N):
        f = 0;
        for l in range(n):
                F = Fl[l];
                k = (n-l)*(n-l);
                r = 0;  ff = 0;
                for i in range(N-n,N):
                        F[i] *= k;
                        F[i-1] = ( _ := F[i]//(N-i) );
                        r = r*l+_;
                        ff = ff*(l+1)+_;
                F[-1] = f-r*l;
                f = ff*(l+1)+F[-1];
        ff = [0]*N;     ff[-1] = f;
        Fl.append(ff);
        res[n] = f;
        tt[n] = (datetime.datetime.now()-t)/dt;
print(len(str(res[-1])));
print(tt[::100]);
 

Результат:

4941
[0, 0.077507, 0.773941, 3.224009, 9.611014, 22.661707, 45.289982, 81.614997, 136.412417, 214.808085, 321.546689]

Null в сообщении #1600012 писал(а):

Сложность в любом случае $O(n^3\times(n))$(на длинную арифметику).

Получилось $\approx 3.8$

Добрый день.
Попробовал Ваш алгоритм. Он работает, но я никак не могу понять - где он держит саму матрицу. Соответственно - ее нельзя поменять на произвольную.
Отсюда возникает два вопроса:
1) Может ли этот алгоритм быть применен для вычисления перманента любых матриц или только той самой, которая задана в топике?
2.1) Если - да, не могли бы Вы подсказать - как его переписать так, чтобы в него можно было вставить любую матрицу (я так понимаю - только квадратную, но и это уже будет хорошо)?
2.2) Если - нет, не могли бы Вы подсказать где можно найти другой достаточно быстрый алгоритм вычисления перманента? Мне не нужно высчитывать большие числа, но перманентов нужно вычислить много (счет как минимум на сотни тысяч) и чем быстрее я буду их считать - тем лучше.

Заранее благодарю.

[Null]

flamehowk в сообщении #1655484 писал(а):

1) Может ли этот алгоритм быть применен для вычисления перманента любых матриц или только той самой, которая задана в топике?

Этот алгоритм работает только для матриц вида:
$ \begin{array}{cccc} a & b& c & d\\ b & b &c &d\\ c & c &c &d\\d&d&d&d\end{array}$
Про общие алгоритмы расчета перманентов ничего не знаю.

[flamehowk]

Благодарю Вас.
Я так и предположил, потому что похоже этот код генерирует матрицу по ходу вычислений.
Что ж - пойду разбираться с формулами Райзера.

[drzewo]

kthxbye в сообщении #1599718 писал(а):

ничения: алгоритм должен вычислять вышеуказанное количество членов максимум за минуту.

а разве это не от железа зависит?

[flamehowk]

Если вдруг кому понадобится универсальный алгоритм для расчета перманентов матриц - вот нашел код сразу на нескольких языках. Вариант на Питоне есть и для Райзера и для Глинна.
https://codegolf.stackexchange.com/questions/97060/calculate-the-permanent-as-quickly-as-possible

[wrest]

flamehowk
По Райзеру, как тут ранее упоминалось, в PARI/GP считает встроенная (а значит - сравнительно быстрая) функция

Цитата:

matpermanent(x)
Permanent of the square matrix x using Ryser's formula in Gray code order.

? n = 20; m = matrix(n,n,i,j, i!=j);
? matpermanent(m)
%2 = 895014631192902121
? n! * sum(i=0,n, (-1)^i/i!)
%3 = 895014631192902121
This function runs in time $O(2^n n)$ for a matrix of size $n$ and is not implemented for $n$ large.