Генерация перестановки (PARI)

[kthxbye]

Пусть $P(n)$ последовательности $s(1),s(2),s(3),...$ находится через оставление $s(1),...,s(n)$ на своих позициях и сортировку в порядке убывания каждых $n$ последовательных членов следующих далее; применим $P(2)$ к $1,2,3,...$ чтобы получить $PS(2)$, затем применим $P(3)$ к $PS(2)$ чтобы получить $PS(3)$, затем применим $P(4)$ к $PS(3)$ и т.д. Предел $PS(n)$ есть $a(n)$.

Последовательность начинается так:

$$1, 2, 4, 6, 10, 12, 18, 22, 30, 34, 42$$

Примеры:

$$1, 2, (4, 3), (6, 5), (8, 7), (10, 9), (12, 11), (14, 13), (16, 15), (18, 17)$$

$$1, 2, 4, (6, 5, 3), (10, 8, 7), (12, 11, 9), (16, 14, 13), (18, 17, 15)$$

$$1, 2, 4, 6, (10, 8, 5, 3), (12, 11, 9, 7), (18, 16, 14, 13)$$

$$1, 2, 4, 6, 10, (12, 11, 8, 5, 3), (18, 16, 14, 9, 7)$$

$$1, 2, 4, 6, 10, 12, (18, 16, 11, 8, 5, 3)$$

Судя по всему $a(n)$ это A002491.

Можно ли реализовать процесс генерации перестановки на PARI? Можно ли также отдельно задавать сколько первых членов остается на своих позициях и каковы условия перестановки следующих?

[kthxbye]

На данный момент имеется всего одно вот такое громоздкое решение:

default(parisizemax,2^8*10^6)
n=8
size=vector(2^n,i,0)
size[2]=2^n
size[3]=floor(2^n/3)
for(i=4,2^n,size[i]=floor(size[i-1]*(i-1)/i))
ps2=vector(2^n,i,0)
ps2[1]=1
ps2[2]=2
for(i=3,2^n,ps2[i]=i-(-1)^i)
\\ блок #3
vs3=vector(size[3],i,vector(3,j,0))
for(i=1,size[3],for(j=1,3,vs3[i][j]=ps2[3*(i-1)+j]))
fs3(n,m)=my(s=Set(vs3[n]));s=setunion(s,Set(m));s[setsearch(s,m)-1]
vs31=vector(size[3],i,vector(3,j,0))
for(j=1,3,vs31[1][j]=vs3[1][j])
for(i=2,size[3],for(j=1,3,vs31[i][j]=if(j==1,vecmax(vs3[i]),fs3(i,vs31[i][j-1]))))
ps3=vector(size[3]*3,i,0)
for(i=1,size[3]*3,ps3[i]=vs31[(i-1)\3+1][(i-1)%3+1])
\\ блок #4
vs4=vector(size[4],i,vector(4,j,0))
for(i=1,size[4],for(j=1,4,vs4[i][j]=ps3[4*(i-1)+j]))
fs4(n,m)=my(s=Set(vs4[n]));s=setunion(s,Set(m));s[setsearch(s,m)-1]
vs41=vector(size[4],i,vector(4,j,0))
for(j=1,4,vs41[1][j]=vs4[1][j])
for(i=2,size[4],for(j=1,4,vs41[i][j]=if(j==1,vecmax(vs4[i]),fs4(i,vs41[i][j-1]))))
ps4=vector(size[4]*4,i,0)
for(i=1,size[4]*4,ps4[i]=vs41[(i-1)\4+1][(i-1)%4+1])
\\ блок #5
vs5=vector(size[5],i,vector(5,j,0))
for(i=1,size[5],for(j=1,5,vs5[i][j]=ps4[5*(i-1)+j]))
fs5(n,m)=my(s=Set(vs5[n]));s=setunion(s,Set(m));s[setsearch(s,m)-1]
vs51=vector(size[5],i,vector(5,j,0))
for(j=1,5,vs51[1][j]=vs5[1][j])
for(i=2,size[5],for(j=1,5,vs51[i][j]=if(j==1,vecmax(vs5[i]),fs5(i,vs51[i][j-1]))))
ps5=vector(size[5]*5,i,0)
for(i=1,size[5]*5,ps5[i]=vs51[(i-1)\5+1][(i-1)%5+1])
print(ps5)

Блоки надо вручную копировать и заменять $k$ на $k+1$ и в одном месте $k-1$ на $k$.

Нельзя ли обойтись без всего этого и сгенерировать перестановку чуть попроще?

[kthxbye]

Чуть-чуть упростил:

default(parisizemax,2^8*10^6)
n=10
m=20
size=vector(2^n,i,0)
size[2]=2^(n-1)
size[3]=floor(2^n/3)
for(i=4,2^n,size[i]=floor(size[i-1]*(i-1)/i))
ps=vector(m,i,vector(size[i]*i,j,0))
for(j=1,2,ps[2][j]=j)
for(j=3,size[2]*2,ps[2][j]=j-(-1)^j)
vs=vector(m,i,vector(size[i],j,vector(i,k,0)))
vs1=vector(m,i,vector(size[i],j,vector(i,k,0)))
f(q,n,m)=my(s=Set(vs[q][n]));s=setunion(s,Set(m));s[setsearch(s,m)-1]
\\ блок #3
for(i=3,3,for(j=1,size[i],for(k=1,i,vs[i][j][k]=ps[i-1][i*(j-1)+k])))
for(i=3,3,for(k=1,i,vs1[i][1][k]=vs[i][1][k]))
for(i=3,3,for(j=2,size[i],for(k=1,i,vs1[i][j][k]=if(k==1,vecmax(vs[i][j]),f(i,j,vs1[i][j][k-1])))))
for(i=3,3,for(j=1,size[i]*i,ps[i][j]=vs1[i][(j-1)\i+1][(j-1)%i+1]))

Изменения относительно незначительные.

Я пробовал запускать сразу

for(i=3,m,for(j=1,size[i],for(k=1,i,vs[i][j][k]=ps[i-1][i*(j-1)+k])))
for(i=3,m,for(k=1,i,vs1[i][1][k]=vs[i][1][k]))
for(i=3,m,for(j=2,size[i],for(k=1,i,vs1[i][j][k]=if(k==1,vecmax(vs[i][j]),f(i,j,vs1[i][j][k-1])))))
for(i=3,m,for(j=1,size[i]*i,ps[i][j]=vs1[i][(j-1)\i+1][(j-1)%i+1]))

Программа успешно выполняет циклы для $i=3$. Однако я предполагал, что при $i=4$ за значением ps[i-1][i*(j-1)+k] она полезет в последний цикл при $i=3$, откуда за значением vs1[i][(j-1)\i+1][(j-1)%i+1] она полезет в третий цикл при $i=3$ и т.д. Вместо этого она сначала отрабатывает весь первый цикл, затем весь второй и т.д.

Можно ли реализовать что-то вроде

for(i=3,m,
{
for(j=1,size[i],for(k=1,i,vs[i][j][k]=ps[i-1][i*(j-1)+k]))
for(k=1,i,vs1[i][1][k]=vs[i][1][k])
for(j=2,size[i],for(k=1,i,vs1[i][j][k]=if(k==1,vecmax(vs[i][j]),f(i,j,vs1[i][j][k-1]))))
for(j=1,size[i]*i,ps[i][j]=vs1[i][(j-1)\i+1][(j-1)%i+1])
})

Если да, то какой синтаксис для этого надо использовать?

[kthxbye]

Полазил по старым темам и нашел решение вот тут. По итогам программа упростилась до

default(parisizemax,2^8*10^6)
g(n,m)=if(n<2,0,if(n==2,2^(m-1),if(n==3,floor(2^m/3),floor(g(n-1,m)*(n-1)/n))))
e(q)={my(n=11,m=77,size=vector(m,i,g(i,n)),ps=vector(m,i,vector(size[i]*i,j,0)),vs=vector(m,i,vector(size[i],j,vector(i,k,0))),vs1=vector(m,i,vector(size[i],j,vector(i,k,0))));
ps=if(q<4,ps,e(q-1));
for(j=1,2,ps[2][j]=j);
for(j=3,size[2]*2,ps[2][j]=j-(-1)^j);
for(i=q,q,
for(j=1,size[i],for(k=1,i,vs[i][j][k]=ps[i-1][i*(j-1)+k]));
for(k=1,i,vs1[i][1][k]=vs[i][1][k]);
my(f(b,c,d)=my(s=Set(vs[b][c]));s=setunion(s,Set(d));s[setsearch(s,d)-1]);
for(j=2,size[i],for(k=1,i,vs1[i][j][k]=if(k==1,vecmax(vs[i][j]),f(i,j,vs1[i][j][k-1]))));
for(j=1,size[i]*i,ps[i][j]=vs1[i][(j-1)\i+1][(j-1)%i+1]);
return(ps))}
print(e(77)[77])

Есть идеи как все это оптимизировать, но это уже скорее всего на завтра.

[kthxbye]

Финальный вариант вот такой:

default(parisizemax,2^8*10^6)
g(n,m)=if(n<2,0,if(n==2,2^(m-1),if(n==3,floor(2^m/3),floor(g(n-1,m)*(n-1)/n))))
e1(q)={my(n=12,m=110,size=vector(m,i,g(i,n)),ps=vector(2^n,i,0));
for(j=1,2,ps[j]=j);
for(j=3,2^n,ps[j]=j-(-1)^j);
ps=if(q<4,ps,e1(q-1));
for(i=q,q,
my(vs=vector(size[i],j,vector(i,k,0)),vs1=vector(size[i],j,vector(i,k,0)));
for(j=1,size[i],for(k=1,i,vs[j][k]=ps[i*(j-1)+k]));
for(k=1,i,vs1[1][k]=vs[1][k]);
my(f(c,d)=my(s=Set(vs[c]));s=setunion(s,Set(d));s[setsearch(s,d)-1]);
my(f1(c,d)=my(s=Set(vs[c]));s=setunion(s,Set(d));s[setsearch(s,d)+1]);
for(j=2,size[i],for(k=1,i,vs1[j][k]=if(isprime(q)==0,if(k==1,vecmin(vs[j]),f1(j,vs1[j][k-1])),if(k==1,vecmax(vs[j]),f(j,vs1[j][k-1])))));
for(j=1,size[i]*i,ps[j]=vs1[(j-1)\i+1][(j-1)%i+1]);
return(ps))}
v=e1(110)
for(i=1,299,print(v[i]))

Дальше наверное уже не упрощается.

В приведенном выше примере я чуть-чуть поменял условия: сортировка по убыванию при простых и по возрастанию при составных. Тогда $a(\operatorname{prime}(n)+1)$ генерирует A073818 (с небольшой разницей в одном члене - $16$ вместо $15$).