Счастливые числа (PARI)

[kthxbye]

Для начала коротенький ликбез: Lucky number (Википедия). Там есть страничка и на русском, только вот ссылка у меня что-то не вставляется.

Решил я значит написать простенькую программку без использования сита. Вот что у меня получилось:

a(n)=if(n<4,2^n-1,{my(b(k)=local(ok = 0, m = 1, A=k, B=0); while (!ok, if ((A-B == 0), ok = 1, B=A; A=if(m==1,A\2+1,A*(a(m)-1)\a(m)+1); m++); ); B;);
my(c(k)=local(ok = 0, m = a(n-1)+3); while (!ok, if ((b(m) == k+1), ok = 1, ok = 0; m++); ); m-1;);
return(c(n))})
for(i=1,299,print(a(i)))

Поначалу выдает неплохо, но потом постепенно ступорится, вероятно из-за повторных вычислений $a(n)$. Можно ли где-нибудь сохранять ее значения, чтобы хоть как-то ускорить процесс? Что еще здесь можно улучшить?

[wrest]

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

Можно ли где-нибудь сохранять ее значения, чтобы хоть как-то ускорить процесс?

По мемоизации см. тут
post1403201.html#p1403201
И тут
post1405301.html#p1405301

Ну и A000959

[kthxbye]

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

Ну и A000959

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

программку без использования сита

А там таких как раз-таки и нету.

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

По мемоизации см. тут post1403201.html#p1403201
И тут post1405301.html#p1405301

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

Единственная функция, которая у меня вызывается повторно - это $a(n)$. Вы не могли бы хотя бы намекнуть какую часть из вашего конечного решения мне можно просто прилепить к своей? Или, например, как вызывать функцию не рекурсивно, а из словаря?

[wrest]

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

Вы не могли бы хотя бы намекнуть какую часть из вашего конечного решения мне можно просто прилепить к своей?

Ну идея такая, что при вычислении a(n) надо проверить не вычислялось ли это раньше.
Для этого создается сперва пустая карта (словарь) M куда будут запоминаться ранее вычисленные значения в виде пар n,a(n)
Создаем карту M=map()
Проверка вычислялась ли уже a(n) делается так.
Внутри функции a() определяем переменную res -- в неё будем записывать результат.
Проверяем
if(mapisdefined(M,n,&res),return(res)) -- да, вычисляли. Тогда сразу же аозвращаем ранее вычисленное значение
Дальше собсно вычисляем a(n), результатом будет наша res
Когда вычислили, перед возвратом запоминаем результат в карту:
mapput(M,n,res)
После чего возвращаем результат вычисления a(n).
return(res)

Недостаток тут в том, что в функции a() используется глобальная переменная M. За этим надо следить -- создать её например, обнулить где-то потом и т.п.

О том как этого избежать, вторая ссылка. Это акттуально для рекурсивных функций, для нерекурсивных вполне можно следить за словарём в стороне.

[kthxbye]

wrest, благодарю за объяснение, идея стала предельно ясна, но как воплотить ее в свой программе так и не додумался. Зато нашел вот такой примитивный вариант:

v1=vector(2^10,i,0)
for(i=1,3,v1[i]=2^i-1)
f(n)=if(n<4,2^n-1,{for(i=n,n,
my(b(k)=local(ok = 0, m = 1, A=k, B=0); while (!ok, if ((A-B == 0), ok = 1, B=A; A=if(m==1,A\2+1,A*(v1[m]-1)\v1[m]+1); m++); ); B;);
my(c(k)=local(ok = 0, m = v1[n-1]+3); while (!ok, if ((b(m) == k+1), ok = 1, ok = 0; m++); ); m-1;);
v1[i]=c(i);
return(v1[i]));})
v2=vector(2^10,i,f(i))
a(n)=v2[n]

[wrest]

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

благодарю за объяснение, идея стала предельно ясна, но как воплотить ее в свой программе так и не додумался.

Я в ваши пргораммы не вникаю и не понимаю :)
Вы писали, что у вас многократно производится вычисление a(n) для одних и тех же n и это замедляет вычисления. Если это так (я этого не вижу), то мемоизация вам поможет. Нет так нет.