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

[Null]

$f_{4,0}=24x^4$
$f_{4,1}=6x^3+18x^2$
$f_{4,2}=6x^2+40x+16$
$f_{4,3}=19x+133$
$f_{4,4}=209$
Похоже работает.

[kthxbye]

Кое-как слепил программку на PARI/GP:

f1(n,l,x)=if(l==0,n!*x^n,f2(n,l,x)+f1(n,l-1,l)-f2(n,l,l))
f2(n,l,m)=if(n>0,my(A=intformal((n-l)^2*f1(n-1,l,x))); sum(j=1,n-l,m^j*polcoeff(A,j)))
f3(n,l,x)=if(l==0,n!*x^n,f1(n,l,x)+f3(n,l-1,l)-f1(n,l,l))
b(n)=f1(n,n,n)

Считает медленно, но значения согласуются. Завтра напишу реализацию с использованием векторов и сравню скорость со своим алгоритмом.

[maxal]

Вот реализация на Sage
Для примера вычисляет несколько первых членов и значение для $n=100$. В онлайне $n=500$ осилить не может, но локально оно тем же кодом вычисляется за 40 секунд.

-- Tue Jul 04, 2023 13:11:43 --

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

Считает медленно, но значения согласуются.

Это из-за рекурсии. Либо разверните ее в цикл, либо используйте мемоизацию.

-- Tue Jul 04, 2023 13:22:21 --

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

maxal, нужно чтобы кто-нибудь проверил хотя-бы.

Проверка прошла успешно.

[wrest]

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

Вот реализация на Sage

Вам, как зачинателю темы «интерактивный курс: введение в программирование на PARI/GP» самое то написать на pari :)

[maxal]

wrest, именно это kthxbye пообещал сделать.

[kthxbye]

Вот реализация с использованием массивов:

f1(n,l,m,A)=sum(j=0,n-l,m^j*polcoeff(A,j))
f2(n,l,m,A)=my(B=(n-l)^2*intformal(A)); sum(j=1,n-l,m^j*polcoeff(B,j))
b(n)=my(v1, v2); v1=vector(n+1,i,i--; i!*x^i); v2=v1; for(i=1,n,for(j=i,n,v2[j+1]=f2(j,i,x,v2[j])+f1(j,i-1,i,v1[j+1])-f2(j,i,i,v2[j])); v1=v2); v1

К сожалению она проигрывает моему алгоритму, поскольку считает гораздо медленнее.

А можно ли на PARI/GP переменной многочлена внутри элемента массива присваивать значение? Причем так, чтобы рядом был (ну т.е. суммировался) многочлен внутри элемента массива, где она остается переменной. Или даже просто не внутри элемента, а скажем внутри какой-нибудь буквы A, которой мы присвоили элемент массива.

[Null]

Не знаю точных алгоритмов PARI/GP но сразу 2 вопроса:
1.У вас подсчет многочлена в точке не оптимальный(Надо по схеме Горнера за $O(n)$).
2. При интегрировании будут нецелые коэффициенты(Это замедляет?), надо в начале умножить на $(n-l)^2$ а потом интегрировать.
3. $f2(j,i,x,v2[j])$ -оно что строит сам многочлен по его коэффициентам после интегрирования? $O(n^2)$?
Наверное так быстрее? Ни когда не прогал на PARI/GP. Как сделать for по убыванию?

f1(n,l,m,A)={s=0;for(j=0,n-l,s=s*m+polcoeff(A,n-l-j));s}
f2(n,l,m,A)=my(B=intformal((n-l)^2*A)); s=0;for(j=0,n-l,s=s*m+polcoeff(B,n-l-j));B-s
b(n)=my(v1, v2); v1=vector(n+1,i,i--; i!*x^i); v2=v1; for(i=1,n,for(j=i,n,v2[j+1]=f1(j,i-1,i,v1[j+1])+f2(j,i,i,v2[j])); v1=v2); v1
b(200)

[kthxbye]

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

1.У вас подсчет многочлена в точке не оптимальный(Надо по схеме Горнера за $O(n)$).

Я не математик, а только экспериментатор. Можно то же самое, но более простыми словами? Вы имеете ввиду что надо использовать s?

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

2. При интегрировании будут нецелые коэффициенты(Это замедляет?), надо в начале умножить на $(n-l)^2$ а потом интегрировать.

Я вынес $(n-l)^2$ за интеграл и проверил. Скорость та же.

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

3. f2(j,i,x,v2[j]) -оно что строит сам многочлен по его коэффициентам после интегрирования? $O(n^2)$?

Да, тут я чуть-чуть не доглядел. Ваша f3(n,l,A)=intformal((n-l)^2*A) здесь очень кстати.

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

Наверное так быстрее? Ни когда не прогал на PARI/GP. Как сделать for по убыванию?

Да, использование f3 (и если я правильно понял, то еще и схемы Горнера) ускоряет процесс. Но все равно $500$ первых значений за минуту на моем ноуте прога не считает. Кажется for по убыванию нет, но можно просто вычислять вот так:

for(j=0,n,v[n-j+1]=A+B)

А зачем вы поправили код? Он же вычислял все правильно:

f1(n,l,m,A)={s=0;for(j=0,n-l,s=s*m+polcoeff(A,n-l-j));s}
f2(n,l,m,A)=my(B=intformal((n-l)^2*A)); s=0;for(j=0,n-l,s=s*m+polcoeff(B,n-l-j));s
f3(n,l,A)=intformal((n-l)^2*A)
b(n)=my(v1, v2); v1=vector(n+1,i,i--; i!*x^i); v2=v1; for(i=1,n,for(j=i,n,v2[j+1]=f3(j,i,v2[j])+f1(j,i-1,i,v1[j+1])-f2(j,i,i,v2[j])); v1=v2); v1

А вот ваша конечная версия уже дает ошибочные значения.

P.S. Ой, не заметил что вы поменяли f2. Тогда все считается правильно, да.

[wrest]

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

А можно ли на PARI/GP переменной многочлена внутри элемента массива присваивать значение? Причем так, чтобы рядом был (ну т.е. суммировался) многочлен внутри элемента массива, где она остается переменной. Или даже просто не внутри элемента, а скажем внутри какой-нибудь буквы A, которой мы присвоили элемент массива.

Допустим у нас в векторе есть коэффициенты полинома $3x^2+2x+10$ и мы хотим вычислить его значение при $x=3$, равное $43$. Это можно сделать так:
? v=[3,2,10];x=3;eval(Pol(v,'x))
%1 = 43
?

-- 05.07.2023, 13:37 --

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

Кажется for по убыванию нет,

Есть forstep().
? forstep(i=10,1,-1,print1(i," "))
10 9 8 7 6 5 4 3 2 1
?

[kthxbye]

wrest, благодарю! Оказывается можно даже без вектора с коэффициентами:

b(n)=my(v1, v2); v1=vector(n+1,i,i--; i!*x^i); v2=v1; for(i=1,n,for(j=i,n,v2[j+1]=my(A=intformal((j-i)^2*v2[j]), f(x)=eval(A), g(x)=eval(v1[j+1])); v2[j+1]=A+g(i)-f(i)); v1=v2); v1

Считает, правда, крайне долго.

[Null]

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

Считает, правда, крайне долго.

Дольше чем было до этого. Оптимизировали.

[wrest]

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

Считает, правда, крайне долго.

И крайне неправильно

[Null]

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

И крайне неправильно

Проверяю на сайте https://pari.math.u-bordeaux.fr/gp.html - начало правильное как A204262

[wrest]

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

Проверяю на сайте https://pari.math.u-bordeaux.fr/gp.html
- начало правильное как A204262

Повезло просто. ТС использует необъявленный x который неизвестно что. У вас он оказался неинициализированным. А мог бы и не оказаться, как у меня. Например попробуйте
x=5;
b(n)=my(v1, v2); v1=vector(n+1,i,i--; i!*x^i); v2=v1; for(i=1,n,for(j=i,n,v2[j+1]=my(A=intformal((j-i)^2*v2[j]), f(x)=eval(A), g(x)=eval(v1[j+1])); v2[j+1]=A+g(i)-f(i)); v1=v2); v1

Впрочем, замена v1=vector(n+1,i,i--; i!*x^i); на v1=vector(n+1,i,(i-1)!); должна помочь с лишним иксом.

Дальше там вообще "непереводимый фольклор", удивительно что оно вообще работает.

[kthxbye]

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

Дальше там вообще "непереводимый фольклор", удивительно что оно вообще работает.

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

ТС использует необъявленный x который неизвестно что.

Он нужен для вычисления интегралов.