Трансформация RECORDS посл.-ти a(n)=isprime(n)+a(n\2)

[kthxbye]

Имеем последовательность (A078349)

$$a(n)=\operatorname{isprime}(n)+a\left(\left\lfloor\frac{n}{2}\right\rfloor\right), a(0)=0$$

где $\operatorname{isprime}(n)=1$, если $n$ - простое.

Если рассматривать трансформацию RECORDS, то бишь такие $n$, при которых $m$ встречается в последовательности впервые, имеем

$$2,5,11,23,47,191,383,1439,2879,11519$$

Это совсем ни о чем не говорит, но, что примечательно, эти первые 10 членов последовательности простые. Меня интересует продолжение, и, собственно, как эти члены емко извлекать.

Вот мой простейший код на pari:

default(parisizemax,10^9)
n=10
v=vector(2^n,i,0)
v[1]=0
for(i=1,2^n-1,v[i+1]=isprime(i)+v[i\2+1])
v1=vector(2^n,i,0)
for(i=1,2^n,v1[i]=if(sum(k=1,i,v[k]==v[i])==1,i-1,0))
print(v1)

Печатал до $2^{15}$, дальше все считается ну очень медленно.

[wrest]

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

дальше все считается ну очень медленно.

Так вы во втором цикле ещё считаете сумму на каждом проходе заново, время счета наверное квадратично растет. А зачем...
Надо избавиться от этого:
for(i=1,2^n,v1[i]=if(sum(k=1,i,v[k]==v[i])==1,i-1,0))

[kthxbye]

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

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

дальше все считается ну очень медленно.

Так вы во втором цикле ещё считаете сумму на каждом проходе заново, время счета наверное квадратично растет. А зачем...
Надо избавиться от этого:
for(i=1,2^n,v1[i]=if(sum(k=1,i,v[k]==v[i])==1,i-1,0))

Покопался в старых прогах OEIS, нашел вот такой вариант:

a(n)=if(n==0,0,isprime(n)+a(n\2))
b(n) = {my(k=0); while (a(k) != n, k++); k; }
for(i=0,15,print(isprime(b(i))))

Что примечательно, продолжает возвращать простые.

[lel0lel]

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

при которых $m$ встречается в последовательности впервые

Что такое $m$, это $a(n)$? Если это верно, то утверждение очевидно. Ведь если $n$ составное, то $a(n)=a\left(\left\lfloor\frac{n}{2}\right\rfloor\right)$

[Dmitriy40]

Не так уж и медленно:

? default(parisize,10^8)
  ***   Warning: new stack size = 100000000 (95.367 Mbytes).
? a=vectorsmall(2^22); a[1]=0; n=0; for(i=2,#a, a[i]=isprime(i)+a[i\2]; if(a[i]>n, print1(i,", "); n++));
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119,
time = 1,202 ms.

? default(parisize,5*10^9)\\GP64
  ***   Warning: new stack size = 5000000000 (4768.372 Mbytes).
? a=vectorsmall(2^28); a[1]=0; n=0; for(i=2,#a, a[i]=isprime(i)+a[i\2]; if(a[i]>n, print1(i,", "); n++));
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127,
time = 1min, 14,022 ms.

? default(parisize,10^8)
  ***   Warning: new stack size = 100000000 (95.367 Mbytes).
? a(n)=if(n<=#v,v[n], isprime(n)+a(n\2));\\n>1!
? v=vectorsmall(16*10^6); v[1]=0; for(i=2,#v, v[i]=isprime(i)+v[i\2]);
time = 3,714 ms.
? n=1; forprime(k=2,2*10^8, if(a(k)!=n, next); print1(k,", "); n++);
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127,
time = 20,124 ms.
? n=1; p=1; forprime(k=2,2*10^8, if(k<p || a(k)!=n, next); print1(k,", "); n++; p=2*k);
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127,
time = 14,867 ms.
? n=1; p=1; forprime(k=2,2*10^8, if(k<p || a(k\2)!=n-1, next); print1(k,", "); n++; p=2*k);
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127,
time = 7,646 ms.

? default(parisize,5*10^9)\\GP64
  ***   Warning: new stack size = 5000000000 (4768.372 Mbytes).
? a(n)=if(n<=#v,v[n], isprime(n)+a(n\2));\\n>1!
? v=vectorsmall(5*10^8); v[1]=0; for(i=2,#v, v[i]=isprime(i)+v[i\2]);
time = 1min, 57,406 ms.
? n=1; forprime(k=2,2*10^8, if(a(k)!=n, next); print1(k,", "); n++);
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127,
time = 3,198 ms.
? n=1; forprime(k=2,oo, if(a(k)!=n, next); print1(k,", "); n++; if(n>17, break));
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127, 1509921587,
time = 1min, 22,181 ms.

-- 26.11.2021, 16:23 --

Иллюстрация насколько к PARI/GP неприменимы обычные методы оптимизации, в данном случае исключение хвостовой рекурсии:

? default(parisize,10^8)
  ***   Warning: new stack size = 100000000 (95.367 Mbytes).
? v=vectorsmall(16*10^6); v[1]=0; for(i=2,#v, v[i]=isprime(i)+v[i\2]);
time = 3,714 ms.

? a(n)=if(n<=#v,v[n], isprime(n)+a(n\2));
? n=1; p=1; forprime(k=2,2*10^8, if(k<p || a(k\2)!=n-1, next); print1(k,", "); n++; p=2*k);
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127,
time = 8,002 ms.

? a(n)=my(x=0); while(n>#v, x+=isprime(n); n=n\2); x+v[n];
? n=1; p=1; forprime(k=2,2*10^8, if(k<p || a(k\2)!=n-1, next); print1(k,", "); n++; p=2*k);
2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127,
time = 30,749 ms.

[Dmitriy40]

Фактически исходная последовательность это количество простых в разложении числа по делению пополам с округлением вниз. Например:

? n=194224127; while(n>1, print1(n,":",isprime(n),",  "); n=n\2);
194224127:1,  97112063:1,  48556031:1,  24278015:0,  12139007:1,
6069503:0,  3034751:1,  1517375:0,  758687:1,  379343:1,  189671:1,
94835:0,  47417:1,  23708:0,  11854:0,  5927:1,  2963:1,  1481:1,
740:0,  370:0,  185:0,  92:0,  46:0,  23:1,  11:1,  5:1,  2:1,

[Dmitriy40]

Если перебор написать с использованием иерархии нескольких десятков решёт Эратосфена, то можно получить скорость перебора чисел порядка нескольких миллиардов в секунду (не на PARI/GP конечно). Т.е. за разумное время (несколько суток счёта) можно найти значения до сотен триллионов (или даже квадриллиона, $10^{15}$).
На PARI/GP тоже можно, но у него страшно медленный перебор простых (forprime) по сравнению с Эратосфеном.

За 2ч счёта на PARI/GP нашлось и следующее число:

2, 5, 11, 23, 47, 191, 383, 1439, 2879, 11519, 23039, 368633, 1474559, 2949119, 47161343, 194224127, 1509921587, 12078595559,

Вероятно следующие значения: ~193млрд, ~773млрд, ~3.1трлн, ~24.7трлн.

[Dmitriy40]

Dmitriy40 в сообщении #1540657 писал(а):

Вероятно следующие значения: ~193млрд, ~773млрд, ~3.1трлн, ~24.7трлн.

Нет, повезло, следующие два значения оказались меньше: 68438456063, 136876912127.

Что интересно, это первый случай когда разложение оканчивается на 3 вместо 2 (что и порушило мои предсказания):

? n=136876912127; while(n>1, print1(n,,":",isprime(n),",  "); n=n\2);
136876912127:1,  68438456063:1,  34219228031:1,  17109614015:0,
8554807007:1,  4277403503:1,  2138701751:1,  1069350875:0,
534675437:1,  267337718:0,  133668859:0,  66834429:0,  33417214:0,  16708607:0,
8354303:1,  4177151:0,  2088575:0,  1044287:1,  522143:0,  261071:1,  130535:0,
65267:1,  32633:1,  16316:0,  8158:0,  4079:1,  2039:1,  1019:1,
509:1,  254:0,  127:1,  63:0,  31:1,  15:0,  7:1,  3:1,

[Dmitriy40]

За три дня счёта нашлось и следующее число: 547507648511. Следующее будет точно не дальше 8.76трлн.

-- 30.11.2021, 15:58 --

Интересно почему нет решений в промежутке $2p+2\ldots4p$, кроме $1439$.

[Dmitriy40]

9 дней счёта понадобилось на нахождение a(22)=8760122376191. Следующее будет точно не дальше чем 70трлн.

[Dmitriy40]

Переписал программу и менее чем за час (вместо двух месяцев!) подтвердил что a(23)=70080979009529. Ещё за 11ч найдено и a(24)=791648291850239. Ещё за 86ч найдено и a(25)=6332662708454399.