Организация нехитрых вычислений

[kthxbye]

Пусть $a(n)$ - это A329369, число перестановок, получаемых из множества $\left\lbrace1,2,\cdots,m\right\rbrace$, с превышающим множеством, которое строится через включение $m-i$ ($0<i<m$) если $b_{i-1}=1$, где $b_kb_{k-1}\cdots b_1b_0$ ($0\leqslantk<m-1$) это двоичное представление $n$.

Превышающее множество перестановки $\pi$, получаемой из множества $\left\lbrace1,2,\cdots,m\right\rbrace$, это множество индексов $i$, таких, что $\pi_i>i$. Смотрите также секцию примеры на страничке A329369 (там вроде все понятно, так что переводить не буду).

Последовательность добавлена в энциклопедию мною, а нашел я ее исходя из абсолютно иного определения, о чем писал на этом форуме ранее.

На данный момент лучшая программа для вычисления членов последовательности без рекурсии базируется на формуле, эквивалентной обратной биномиальной трансформации по модулю $2$ последовательности A284005.

Первый значимый член у $a(n)$ - это $a(0)$. Однако если сопоставлять $a(n)$ и двоичное представление $n+1$ можно обнаружить много интересного.

Важно отметить, что работать мы будем с двоичным представлением исключительно нечетных чисел.

Например, если известен член с номером в двоичной $101$, то можно без труда с помощью одной формулы (использующей две легко генерируемые вспомогательные последовательности) найти значение члена последовательности, сопоставимого с $11(0\cdots0)101$, где число нулей в скобках может быть любым.

Аналогично можно получить $1(0\cdots0)101$ если известны члены сопоставимые с $101$ и $1101$.

Отсюда вытекает нехитрая задача. Например я хочу узнать номер члена, сопоставимого с $100110001001$. Как это сделать?

Чтобы получить $100110001001$, надо знать $110001001$ и $1110001001$.
Чтобы получить $110001001$, надо знать $1001$.
Чтобы получить $1110001001$ надо знать $10001001$.
Чтобы получить $1001$, надо знать $1$ и $11$.
Чтобы получить $10001001$ надо знать $1001$ и $11001$.
Как получать $1001$ мы уже знаем.
Чтобы получить $11001$, надо знать $1$.

Здесь надо отметить, что члены сопоставимые со сплошной последовательностью единиц в двоичной имеют очень простую формулу.

Встает вопрос: как организовать вычисления, чтобы сначала считались самые маленькие нужные члены, а уже из них потом, соответственно, те, что постарше? Кроме того, как мы видим, некоторые члены могут вызываться повторно, т.е. необходимо также как-то избегать повторных вычислений.

Все это желательно иметь на PARI, чтобы можно было интегрировать непосредственно с формулами.

[Dmitriy40]

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

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

forstep(i=1,10^6,2, SolveNumber(i));

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

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

А это я, и не только я, Вам уже несколько раз показывал: можно сделать табличку/вектор и хранить в ней уже вычисленное, а ранее невычисленное вычислить и положить в неё, можно вместо вектора применить map (wrest тоже неоднократно показывал Вам как) и хранить в ней только уже вычисленные значения.
Т.е. задача банальна и уже много раз решалась, даже конкретно для Вас.

[kthxbye]

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

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

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

forstep(i=1,10^6,2, SolveNumber(i));

Благодарю за комментарий! Признаться, я ожидал его получить либо от вас, либо от wrest'a.

Вы абсолютно верно поняли, что мы работаем с нечетными числами, ввиду чего предложили мне вычисления с шагом $2$. Но я, вероятно, задал неправильный вопрос.

Попробую еще раз. Вот самое важное:

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

Отсюда вытекает нехитрая задача. Например я хочу узнать номер члена, сопоставимого с $100110001001$. Как это сделать?

Чтобы получить $100110001001$, надо знать $110001001$ и $1110001001$.
Чтобы получить $110001001$, надо знать $1001$.
Чтобы получить $1110001001$ надо знать $10001001$.
Чтобы получить $1001$, надо знать $1$ и $11$.
Чтобы получить $10001001$ надо знать $1001$ и $11001$.
Как получать $1001$ мы уже знаем.
Чтобы получить $11001$, надо знать $1$.

Здесь надо отметить, что члены сопоставимые со сплошной последовательностью единиц в двоичной имеют очень простую формулу.

Как для всех нечетных чисел получать вышеуказанные цепочки того из чего они будут вычисляться, и как потом уже, собственно, в обратной последовательности их вычислять? Я предполагаю, что здесь нужно дерево, например двоичное (где в некоторых случаях последующие ветки будут просто нулевые), но как это организовать просто не представляю.

Я кое-как интуитивно написал программки для двух вспомогательных последовательностей, а тут ну просто увольте.

Буду крайне признателен за помощь.

[Dmitriy40]

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

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

Я по прежнему не понимаю проблемы. Вы же знаете как любую цепочку (т.е. число) превратить в меньшие цепочки/числа, значит это обычный рекурсивный спуск.
А чтобы не спускаться каждый раз вплоть до 1 можно заранее вычислить (уже любым способом) несколько тысяч/миллионов первых значений и положить в вектор, а при спуске проверять не попали ли ниже вычисленной ранее границы и если попали то сразу возвращать результат из вектора вместо продолжения спуска, конечно если он там реально есть, если же нету, то продолжить спуск, а при возвратном поднятии положить новый результат на его место в вектор.
Как Вы там делите одно число на несколько меньших я даже не смотрел, это не принципиально и на организацию спуска и исключения повторных вычислений не влияет.

Примерная реализация (банальна):

v=vector(1000,i,-1);\\Сначала ничего не вычислено
SolveNumber(x)=my(t); if(x<=#v && v[x]!=-1, return v[x]); t=SolveNumber(x-1)+SolveNumber(x-2); if(x<=#v, v[x]=t); return t;

Всё!
Как именно будете получать t из x меня не касается, тут написал просто от фонаря для примера, главное положить его в вектор на будущее.
Это исключение повторных вычислений.
Ну а расчёт от меньших к большим реализуется циклом forstep как писал выше. Например сначала для всех помещающихся в v[].

[wrest]

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

Признаться, я ожидал его получить либо от вас, либо от wrest'a.

Я не понимаю чего вы хотите, поэтому мне нечего советовать. Чтобы что-то понять, надо читать ваш пост от царя Гороха (пройти ссылки и т.п.), если бы вы могли внятно формулировать именно вычислительную (алгоритмическую) задачу, можно было бы и подумать. А что значат ваши нули и единички и например какой смысл имеет выражение "Чтобы получить 11001, надо знать 1", находится за пределами моего понимания..

[kthxbye]

Dmitriy40, еще раз благодарю! Буду разбираться.

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

А что значат ваши нули и единички и например какой смысл имеет выражение "Чтобы получить 11001, надо знать 1", находится за пределами моего понимания..

Допустим мы ищем член с номером $25$ ($11001$ в двоичной). Тогда мы смотрим на двоичное представление. Если у нас слева есть две единицы, а дальше нули, то мы все это обрасываем. Остается $1$ в двоичной, т.е. чтобы вычислить $a(25)$ нам достаточно знать $a(1)$ и применить формулу номер один. Если же мы вычисляем член с номером $9$ ($1001$ в двоичной), то мы видим, что в двоичном представлении слева всего одна единица. Тогда нам нужно уже не один, а два члена. Первый мы получаем отбрасывая единицу и нули, т.е. остается $1$ в двоичной, а для второго мы дописываем к первому справа одну единицу, т.е. $11$ в двоичной. Это значит, что чтобы вычислить $a(9)$, нам достаточно знать $a(1)$ и $a(3)$ и применить формулу два.

Все то же самое вроде бы внятно описано вот тут:

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

Например, если известен член с номером в двоичной $101$, то можно без труда с помощью одной формулы (использующей две легко генерируемые вспомогательные последовательности) найти значение члена последовательности, сопоставимого с $11(0\cdots0)101$, где число нулей в скобках может быть любым.

Аналогично можно получить $1(0\cdots0)101$ если известны члены сопоставимые с $101$ и $1101$.

Здесь важно отметить, что конечные члены, из которых мы будем вычислять, должны представлять из себя сплошные последовательности единиц в двоичной, т.е. это $a(1), a(3), a(7), \cdots, a(2^k-1)$.

Теперь нам дано любое нечетное число. Надо составить какую-то карту того, из каких членов оно будет вычисляться, а потом с конца начать ее вычислять.

[Dmitriy40]

kthxbye
Вот я возьмусь и детально прокомментирую что нифига непонятно из Вашего текста. Не с целью получить ответы, роазбираться я всё равно не буду, а показать что надо гораздо аккуратнее описывать задачу.

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

Если у нас слева есть две единицы, а дальше нули, то мы все это обрасываем.

Что именно "все это"?

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

т.е. чтобы вычислить $a(25)$ нам достаточно знать $a(1)$ и применить формулу номер один.

Формулы номер один в текстах выше нет.

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

то мы видим, что в двоичном представлении слева всего одна единица. Тогда нам нужно уже не один, а два члена.

И снова непонятно почему два.

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

Первый мы получаем отбрасывая единицу и нули,

Слева или справа или из середины? А если там не две единицы, а 55 и 1100 нулей, то что докуда отбрасывать?

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

и применить формулу два.

Формулы два в тексте тоже нет.

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

Все то же самое вроде бы внятно описано вот тут:

Похоже это внятно только для Вас.

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

то можно без труда с помощью одной формулы (использующей две легко генерируемые вспомогательные последовательности) найти значение члена последовательности, сопоставимого с $11(0\cdots0)101$, где число нулей в скобках может быть любым.

Это ну совершенно не "внятно". Какие-то неизвестные две новые последовательности в новой формуле, непонятно что значит "сопоставимого с".

Ну и ходить по OEIS, да ещё и по примерам в них, чтобы понять что же Вы хотите ... Можно было всё важное и необходимое сформулировать сразу здесь (кажется это даже правилами требуется).
Тем более что Ваш вопрос не как что-то там считать, а как организовать цикл и рекурсию. Это вообще не про математику и OEIS и последователньости. И они тут в общем лишние.

[kthxbye]

Хорошо, по всей видимости объяснить все у меня действительно не получилось.

Давайте еще раз с самого начала. У нас есть $a(n)$. Что она из себя представляет описано подробно в первом посте. Первый значимый член у $a(n)$ - это $a(0)$. Наипримитивнейшая рекурсия базируется на том, что задано $a(0)=1$.

Введем дополнительно $b(n)=a(n-1)$, чтобы работать с последовательностью, у которой первый значимый член это $b(1)$. Это позволит нам увидеть новые закономерности.

Пусть также $a_1(n)$ - это A284005. Не суть важно что это за последовательность. Нам нужно знать только, что $a(n)$ можно получать через обратную биномиальную трансформацию по модулю $2$ последовательности $a_1(n)$:

$$a(n)=\sum\limits_{k=0}^{n}(-1)^{\operatorname{wt}(n)-\operatorname{wt}(k)}\left[\binom{n}{k}\bmod 2 \right]a_1(k)$$

где $\operatorname{wt}(n)$ - это A000120, двоичный вес $n$ или же число единиц в двоичной записи $n$.

Перебирать все $k$, проверяя остаток от деления бинома на $2$ - это долго и неудобно. На этот случай существует эквивалентная формула, где слагамые $a_1(k)$ генерируются сразу для $\binom{n}{k}\bmod 2 = 1$. Вот как она выглядит на страничке A329369:

Цитата:

a(n) = Sum_{j=0..2^wt(n) - 1} (-1)^(wt(n) - wt(j)) Product_{k=0..wt(n) - 1} (1 + wt(floor(j/2^k)))^T(n,k) for n > 0 with a(0) = 1 where wt(n) = A000120(n), T(n,k) = T(floor(n/2), k - n mod 2) for k > 0 with T(n,0) = A001511(n) (equivalent to theorem 6.3 at page 296, see R. Ehrenborg and E. Steingrimsson link). Here T(n, k) - 1 for k > 0 is the length of the run of zeros between k-th pair of ones from the right side in the binary expansion of n. Also this formula is equivalent to inverse modulo 2 binomial transform of A284005.

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

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

$$b(2n-1+3\cdot 2^{\ell(2n-1)+k})=3^kc_1(n)-2^kc_2(n)+b(2n-1)$$

Т.е. если $2n-1=1xyz1_2$ (двоичная запись), то

$$2n-1+3\cdot 2^{\ell(2n-1)+k}=11\underbrace{0\cdots 0}_{k-1}1xyz1_2$$

Т.е. мы можем как бы "нарастить" в двоичной записи слева две единицы и $k-1$ нулей.

Аналогично существует вторая формула (про которую также нигде ничего не было сказано):

$$b(2n-1+2^{\ell(2n-1)+k})=2^{k-1}\left(b(2n-1)+b(2n-1+2^{\ell(2n-1)+1})\right)-b(2n-1)$$

Здесь все также $2n-1=1xyz1_2$, но уже

$$2n-1+2^{\ell(2n-1)+k}=1\underbrace{0\cdots 0}_{k-1}1xyz1_2$$

Т.е. "наращивается" всего одна единица и все также $k-1$ нулей.

Обе формулы справедливы для $n>0, k>0$, однако вторая при $k=1$ не имеет смысла (легко догадаться почему).

В обоих формулах $\ell(n)=\left\lfloor\log_2 n\right\rfloor$. Кроме того, $c_1(n)$ и $c_2(n)$ это две вспомогательные последовательности, которые генерятся относительно несложно, если только понимать этот механизм.

Теперь, собственно, программки:

default(parisizemax,2^8*10^6)
l(n)=logint(n,2)
a(n)=my(n=(n+1)/2^valuation(n+1, 2)-1, wt=hammingweight(n), v=[], A=0); if(n==0, 1, for(i=0, logint(n, 2), if(bittest(n, i), v=concat(v, A); A=0, A++)); sum(j=0, 2^wt-1, (-1)^(wt - hammingweight(j))*prod(k=0, wt-1, (1 + hammingweight(j\2^k))^(v[k+1] + 1))))
b(n)=a(n-1)
f1(n)=my(m=2*n-1, A=0, v=[], v1); for(i=0,l(m),if(bittest(m,i),v=concat(v,A); A=0, A++)); v=Vecrev(v); if(n==1, [8,6], v1=vector(#v,i,vector(i,j,0)); for(j=1,#v,v1[#v][j]=[2*(#v-j+4),2*(#v-j+3)]); for(i=1,#v-1,for(j=1,#v-i,v1[#v-i][j]=[(#v-j+4-i+1)^(v[#v-i]+1)*v1[#v-i+1][j][1]-(#v-j+3-i+1)^(v[#v-i]+1)*v1[#v-i+1][j][2],(#v-j+3-i+1)^(v[#v-i]+1)*v1[#v-i+1][j+1][1]-(#v-j+2-i+1)^(v[#v-i]+1)*v1[#v-i+1][j+1][2]])); v1)
f2(n)=my(m=2*n-1, A=0, v=[], v1); for(i=0,l(m),if(bittest(m,i),v=concat(v,A); A=0, A++)); v=Vecrev(v); if(n==1, [6,6], v1=vector(#v,i,0); v1[#v]=[9,6]; my(B=f1(n), v2=[]); for(i=1,#B,v2=concat(v2,B[i][i][2])); for(i=1,#v-1,v1[#v-i]=[v2[#v-i]*3/2,3^(v[#v-i]+1)*v1[#v-i+1][1]-2^(v[#v-i]+1)*v1[#v-i+1][2]]); [v1[1][1]*2/3,v1[1][2]])
f3(n,k)=my(A=b(2*n-1), B=f2(n)); 3^k*B[1]-2^k*B[2]+A
f4(n,k)=my(A=b(2*n-1), B=b(2*n-1+2^(l(2*n-1)+1))); 2^(k-1)*(A+B)-A
m=5
\\ for(i=1,299,print(f3(i,m)-b(2*i-1+3*2^(l(2*i-1)+m)))) \\ проверка первой формулы где m может быть любым натуральным
\\ for(i=1,299,print(f4(i,m)-b(2*i-1+2^(l(2*i-1)+m)))) \\ проверка второй формулы где m может быть любым натуральным
f5(n)=my(v=[1,1]); for(i=1,n-1,v=concat(v,[0,1,1])); v=Vecrev(v); fromdigits(v,2)
\\ for(i=1,100,print(b(f5(i)))) \\ иллюстрация того, как медленно считается по исходной формуле
v=vector(100,i,0)
v[1]=3
\\ for(i=2,100,my(B=f2((f5(i-1)+1)/2)); v[i]=3^2*B[1]-2^2*B[2]+v[i-1]) \\ иллюстрация выгоды от использования первой формулы (считает 99 членов за 6 секунд)

Эффективность первой формулы проиллюстрирована на примере вычисления члена с номером в двоичной

$$\underbrace{(110)(110)\cdots(110)}_{k}11_2$$

Теперь еще раз о том, что мне нужно. У меня есть две формулы. Если все пробеги единиц в двоичной записи $2n-1$ имеют четную длину (кроме крайнего справа, там это не важно), то мы просто применяем многократно первую формулу и имеем результат. В процессе ее применения мы вычисляем все последовательно, наращивая каждый раз пару единиц слева (и нули если нужно). Это рекурсия в определенном смысле, но не такая, что мы создаем массив из $1000$ элементов и проверяем каждое на то, вычислено оно или нет (и если нет, то вычисляем), а как бы последовательно наращиваем значения не для всех, а лишь для выборочных членов последовательности. Как они будут выбираться, очевидно из двоичной записи $2n-1$.

И теперь, собствено, вопрос (в который раз уже задаю его): как организовать вычисления, чтобы применяя либо первую, либо вторую формулу нарастить к сплошной последовательности единиц справа все необходимые нули и единицы слева?

[Dmitriy40]

Правильно ли я понял что для любого числа $xxx_2$ и любого $k>0$ есть всего два варианта трансформации исходного числа в новое, это $xxx_2 \to 110xxx_2$ и $xxx_2 \to 10xxx_2$ (где количество нулей равно $k$ или возможно $k-1$ в Вашем понимании числа $k$)? Тогда весь Ваш текст до последнего вопроса с кучей формул и дополнительными последовательностями и отсылками к OEIS вообще не нужен, достаточно всего лишь одного этого предложения. Предположим что я понял правильно.

Т.е. задача состоит в получении всех вариантов дополнения слева чисел и желательно в порядке возрастания. ОК.
Рассмотрим все возможные варианты дополнения слева для $0<k<5$ (далее приму для простоты что $k$ это количество нулей) в порядке возрастания величины получаемых чисел:
$xxx_2 \to 10xxx_2, k=1$
$xxx_2 \to 100xxx_2, k=2$
$xxx_2 \to 110xxx_2, k=1$
$xxx_2 \to 1000xxx_2, k=3$
$xxx_2 \to 1100xxx_2, k=2$
$xxx_2 \to 10000xxx_2, k=4$
$xxx_2 \to 11000xxx_2, k=3$
$xxx_2 \to 110000xxx_2, k=4$
Легко заметить что для добавления слева ровно $n>2$ битов/цифр есть ровно два варианта:
$xxx_2 \to 100...0xxx_2, k=n-1$
$xxx_2 \to 110...0xxx_2, k=n-2$

Соответственно для получения из любого числа $xxx_2$ всех возможных однократных дополнений слева достаточно перебирать в порядке возрастания $n>2$ (случай $n=2$ проще обработать отдельно) и дописывать слева числа $\{2,3\}\cdot 2^{n-2}$ именно в таком порядке.

Для многократного дописывания, например для $n\ge4$, кроме двух вариантов однократного дописывания будут ещё и подварианты многократного. Если нужно получить их все, то проще всего рекурсивно (да, можно и без рекурсии, но "вам шашечки или ехать"?). Вопрос как их получать строго в порядке возрастания величины чисел вполне решаем, хоть и не совсем прост.
Но я бы задался вопросом а так ли вообще необходимо получать значения строго в порядке возрастания? Нельзя ли их сначала получить (например для конкретного $n$), а потом отсортировать?

Ну и дополнительный вопрос: а до скольких битов вообще надо считать? 10, 100, 1000, миллион, миллиард, пока не надоест? Зависимость грубо говоря вроде бы квадратичная, т.е. при увеличении длины ($k$ или $n$ без разницы) в битах вдвое потребное время растёт в 4 раза, увеличение длины в 10 раз требует в 100 раз больше времени.


Если же я понял неправильно и надо разбираться в этом нагромождении малопонятных формул и последовательностей ... то я пас. По опыту почти наверняка всё сведётся к чему-то довольно простому (излишняя математичность записи хороша для анализа, но понимание только запутывает).

[wrest]

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

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

Ну вот это примерно поняно (а все что до этого - не понятно и мной пропущено мимо ушей). Нам нужны "пробеги единиц" четной длины. Т.е. мы смотрим на двоичную запись и считаем сколько у нас единиц идут подряд между нулями. Если все эти количества чётные - нам повезло.
Вот вам функция lucky(n) которая вычисляет "повезло" или нет.

(Оффтоп)

lucky(n)=my(dg=binary(n),Z=List());print(n, " converted to binary:",dg);for(i=1,#dg,if(dg[i]==0,listput(Z,i)));print("found ", #Z," zeroes total");print("zero positions:",Z);listinsert(Z,0,1);listput(Z,#dg+1);for(i=2,#Z-1,print("testing ",i-1,"th zero");gap=Z[i]-Z[i-1];if(gap==1,print("adjacent zeroes, go next");next());if( gap%2==1,print("gap is ", gap-1, " even"),print("gap is ", gap-1, " odd");print("result: odd number of ones sequence found, no luck!");return(0)));print("result: bingo!");return(1)

Эта функция возвращает 1 (т.е. "повезло!"), если в двоичной записи n количество единиц между нулями везде четное, и 0 (т.е. "не повезло") если есть и нечетное количество единиц между нулями.

Условие "кроме крайнего справа" учитывается, то есть просматривается вся двоичная запись кроме единиц после последнего нуля. Сответственно, если в двоичной записи только единицы (т.е. число вида $2^n-1$), то они все "крайние справа", и возвращается "повезло".
Функция по ходу дела печатает что и как, соответственно все print(); можно убрать, сильно сократив таким образом текст.

В общем, чем мог - помог :)
Знач там обратите внимание. Применяется список (List()). Это структура в которую можно что-нибудь вставить или удалить. Мне кажется из вашего текста что вы это хотите делать с двоичной записью ("наращивать"). Я использую List чтобы в нём хранить позиции нулей в двоичной записи. Не использую вектор т.к. изначально неясно сколько элементов там будет. После того как список сделал, в служеных целях дополняю его нулями справа и слева (чтоб лишней логики не городить).

Ну может вам это и не нужно, я не знаю.

[kthxbye]

wrest, спасибо за помощь! Буду разбираться.

Dmitriy40, Вы поняли все в принципе правильно. Надо отметить только, что

• $xxx_2$ исключительно нечетное
• подставлять можно любые $k>0$, а число нулей равно $k-1$
• вспомогательные последовательности для больших $\operatorname{wt}(n)$ требуют (пусть и на время) заполнение большого числа массивов, так что мы упираемся в предел памяти, выделяемой для этого

Было бы классно, если бы мы могли наращивать одну единицу без нулей. Но это невозможно, поскольку по второй формуле мы получаем два одинаковых неизвестных члена с обоих сторон.

Я пока слепил вот такую хреновину:

f6(n)=my(n=n/2^valuation(n,2), A=0, B, C, D, E, F, v, v1); while(bittest(n,A), A++); v=[2^A-1]; v1=[2^A-1]; if(A<l(n), v=concat(v,2^(A+1)-1); v1=concat(v1,2^(A+1)-1)); while(A<l(n), A++; while(!bittest(n,A), A++); B=sum(k=0,A,bittest(n,k)*2^k); C=B-2^l(B); D=C+2^(l(C)+1); if(!bittest(n,A-1), if(setsearch(v,C)>0 && setsearch(v,D)>0, v=concat(v,B); v1=concat(v1,2^(l(B)-l(C)-1)*(v1[setsearch(v,C)]+v1[setsearch(v,D)])-v1[setsearch(v,C)]), E=C-2^l(C); F=l(C)-l(E)-1; v=concat(v,D); v1=concat(v1,3^F*f2((E+1)/2)[1]-2^F*f2((E+1)/2)[2]+v1[setsearch(v,E)]); v=concat(v,B); v1=concat(v1,2^(l(B)-l(C)-1)*(v1[setsearch(v,C)]+v1[setsearch(v,D)])-v1[setsearch(v,C)])), C=C-2^l(C); E=C; F=l(B)-l(C)-1; v=concat(v,B); v1=concat(v1,3^F*f2((E+1)/2)[1]-2^F*f2((E+1)/2)[2]+v1[setsearch(v,E)]))); v1[#v1]

Дает много верных результатов, но и ошибочных тоже не мало.

Я решил на всякий случай (а точнее не имея альтернативы) вычислять наращивание для каждой отдельной единицы. Все это я разбил на три вида вычислений:

• вычисление по второй формуле, если оба члена, необходимых для вычисления известны
• предвычисление большего (второго) неизвестного члена, необходимого для вычисления, и потом уже вычисление по второй формуле
• вычисление по первой формуле

Формулы для вспомогательных последовательностей верны. Значит в каком-то из трех видов вычислений есть ошибка. Это надо вероятно смотреть на двоичное представление номеров членов, вычисленных с ошибкой, пошагово вычислять все вручную на бумажечке, и потом уже соображать что к чему.

[kthxbye]

Действительно, теперь все верно:

f6(n)=my(n=n/2^valuation(n,2), A=0, B, C, D, E, F, v, v1); while(bittest(n,A), A++); v=[2^A-1]; v1=[2^A-1]; if(A<l(n), v=concat(v,2^(A+1)-1); v1=concat(v1,2^(A+1)-1)); while(A<l(n), A++; while(!bittest(n,A), A++); B=sum(k=0,A,bittest(n,k)*2^k); C=B-2^l(B); D=C+2^(l(C)+1); if(!bittest(n,A-1), if(setsearch(v,C)>0 && setsearch(v,D)>0, v=concat(v,B); v1=concat(v1,2^(l(B)-l(C)-1)*(v1[setsearch(v,C)]+v1[setsearch(v,D)])-v1[setsearch(v,C)]), E=C-2^l(C); F=l(C)-l(E); v=concat(v,D); v1=concat(v1,3^F*f2((E+1)/2)[1]-2^F*f2((E+1)/2)[2]+v1[setsearch(v,E)]); v=concat(v,B); v1=concat(v1,2^(l(B)-l(C)-1)*(v1[setsearch(v,C)]+v1[setsearch(v,D)])-v1[setsearch(v,C)])), C=C-2^l(C); E=C; F=l(B)-l(C)-1; v=concat(v,B); v1=concat(v1,3^F*f2((E+1)/2)[1]-2^F*f2((E+1)/2)[2]+v1[setsearch(v,E)]))); v1[#v1]

Надо было просто поменять F=l(C)-l(E)-1 на F=l(C)-l(E), и все заработало.

До сих пор не понимаю, как я написал эту и предыдущие программы. Я делаю это как-то интуитивно и если меня даже сейчас заставить написать программу сначала, это все равно отнимет много времени.

Осталось только сделать оптимизацию. Это надо разбираться, какие именно наращивания одной единицы можно не делать.

[kthxbye]

Все оказалось куда тривиальнее, чем я ожидал.

Используя метод генерации первой вспомогательной последовательности, я случайно дошёл до того, что $a(2n)=T(n,1)$, где

$$T(n,k)=(k+1)^{f(n)+1}T(g(n),k+1)-k^{f(n)+1}T(g(n),k), T(0,k)=k$$

Здесь $f(n)$ это A290255 и $g(n)$ это A053645.