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

[Null]

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

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

Разве это не свободная переменная?

Цитата:

2.3.10 Polynomials (t_POL). Type the polynomial in a natural way, not forgetting to put a “∗”
between a coefficient and a formal variable;
? 1 + 2*x + 3*x^2
%1 = 3*x^2 + 2*x + 1
This assumes that x is still a ”free variable”.

[wrest]

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

Разве это не свободная переменная?

Свободная, пока вы ей не присвоите что-нибудь (об этом и говорит "This assumes that x is still a ”free variable”."). Соответственно, ваша программа должна за этим следить.
Например, "освободить" через kill(x) или x='x

[maxal]

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

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

eval лучше не использовать для вычисления значений многочленов (оно не использует схему Горнера). Подставлять значения $x$ в $f(x)$ можно так: subst(f,x,3) или если нет уверенности в имени переменной, то subst(f,variable(f),3).

[kthxbye]

maxal, спасибо за информацию!

Тогда вот это наверное лучший алгоритм на PARI/GP, работающий по методу участника Null:

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])); v2[j+1]=A+subst(v1[j+1],x,i)-subst(A,x,i)); v1=v2); v1

На моем ноуте он считает $400$ членов за $43$ секунды. Примерно за это же время мой алгоритм вычисляет $500$ членов.

[maxal]

Вот мой код на PARI/GP. На моем компе вычисляет 500 членов за 8 секунд.

{ a204262terms(N) = my(a,f,g,x='x);
  a = vector(N);
  a[1] = 1;
  f = [x,1];
  for(n=2,N,
    g = List([n!*x^n]);
    for(l=1,n-1,
        listput(g, intformal((n-l)^2*f[l+1],x));
        g[#g] += subst(g[#g-1]-g[#g],x,l);
    );
    listput(g, a[n] = subst(g[#g],x,n) );
    f = g;
  );
  a;
}

[Null]

Сложность в любом случае $O(n^3\times(n))$(на длинную арифметику).
kthxbye, какой у вас алгоритм? Или подождем еще решателей?

[Geen]

Алгоритм 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$

[wrest]

Ну в общем у меня считает все b до 500 за 14 секунд

? b(n)=my(v1, v2,x='x,f11,f22,f24); v1=vector(n+1,i,(i-1)!); v2=v1; for(i=1,n,for(j=i,n,x=i;f11=eval(v1[j+1]);x='x;if(i==j,f22=0,f22=(i-j)^2*intformal(v2[j]));x=i;f24=eval(f22);x='x;v2[j+1]=f22+f11-f24); v1=v2); return(v1);
? b(500);
? ##
  ***   last result computed in 13,982 ms.
?

-- 05.07.2023, 20:08 --

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

Вот мой код на PARI/GP. На моем компе вычисляет 500 членов за 8 секунд.

На мой планшет похоже, где-то 9 секунд ваш код у меня бежит.

-- 05.07.2023, 20:10 --

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

eval лучше не использовать для вычисления значений многочленов (оно не использует схему Горнера). Подставлять значения $x$ в $f(x)$ можно так: subst(f,x,3) или если нет уверенности в имени переменной, то subst(f,variable(f),3).

Да, я мучился с subst по всякому, но так и не победил...

-- 05.07.2023, 20:16 --

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

На моем ноуте он считает $400$ членов за $43$ секунды.

Хм... у вас небыстрый ноут, однако. Ваш код до 400 у меня проходит на планшете за 7 секунд...

[Null]

(A204264)

Можно аналогично A204264. Только отдельный элемент считает.

{ a204264(N) = my(a,f,g,x='x);
  a = vector(N);
  a[1] = 1;
  f = [x,N];
  for(n=2,N,
    g = List([n!*x^n]);
    for(l=1,n-1,
        listput(g, intformal((n-l)^2*f[l+1],x));
        g[#g] += subst(g[#g-1]-g[#g],x,N+1-l);
    );
    listput(g, a[n] = subst(g[#g],x,N+1-n) );
    f = g;
  );
  a[N];
}

[wrest]

Null
Что-то не так. Ошибку ищите
? print(a204264(5))
78745
?

P.S. А, сорри, у вас другая последовательность, не та что в стартовом посте у ТС. Тогда ошибки нет.

-- 05.07.2023, 20:59 --

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

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

Кстати, а сколько цифр в числе $b(500)$?

Судя по всему $2172$.

2167 вышло.

[kthxbye]

Стоило назвать алгоритм лучшим, как сразу же появились альтернативы (что хорошо и даже очень).

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

kthxbye, какой у вас алгоритм? Или подождем еще решателей?

У меня вот такой:

b(n)=my(v1, v2, v3, v4); v1=vector(n+1, i, 1); v2=v1; v3=vector(n+1, i, 0); v3[1]=1; v4=vector(n, i, vector(i+1, j, binomial(i, j-1)*j)); for(i=1, n, for(q=0, n-i, v2[q+1]=sum(j=0, q+1, v4[q+1][j+1]*v1[j+1])); v1=v2; v3[i+1]=v1[1];); v3

Он базируется на простой и элегантной формуле

$$R(n,q)=\sum\limits_{j=0}^{q+1}\binom{q+1}{j}(j+1)R(n-1,j), R(0,q)=1$$

где мы имеем $R(n,0)=b(n+1)$. К сожалению как его генерализировать хотя бы на ту же A204264 совершенно не представляю.

Может кто-нибудь оценить сложность моего алгоритма с помощью нотации "О" большое? Еще желательно проверить его скорость на современных устройствах.

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

Хм... у вас небыстрый ноут, однако.

Он очень старый, но там можно (пусть и с лагами) смотреть ютуб. В основном я использую его для математических экспериментов и переписки.

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

2167 вышло.

У вас в результате просто нумерация начинается с двух единиц вместо одной. Имеем $b(1)=1, b(2)=3$. Но при этом, конечно, $b(0)=1$.

[wrest]

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

У меня вот такой:

Это немного медленнее (процентов на 10-15%) кода ув. maxal, по крайней мере на моём планшете.

[Null]

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

Может кто-нибудь оценить сложность моего алгоритма с помощью нотации "О" большое?

Тоже $O(n^3\times(n))$. Мы должны посчитать все $R(m,p)$ где $m+p\le n$ и на каждый $O(n)$ операций по формуле. И все та же длинная арифметика $O(n)$(на самом деле $O(n\log(n))$ как и у меня)

[Geen]

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

на самом деле $O(n\log(n))$

На самом деле меньше - у нас только сложение длинное, а умножение/деление везде только длинного числа на короткое.

[wrest]

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

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

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

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

Вот мой код на PARI/GP. На моем компе вычисляет 500 членов за 8 секунд.

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

-- 06.07.2023, 01:39 --