\\ Kevin Ryde, April 2024 \\ \\ This is some PARI/GP code to calculate an individual term of \\ \\ A105316 = concat( S(1), S(2), ... ) \\ where vectors \\ S(1) = [1] vector of a single 1 \\ S(k) = concat( S(k-1), expand(S(k-1)) ) \\ and expand(v) = replace each term in vector v by \\ 1 -> 2,3 \\ 2 -> 3,4 \\ 3 -> nothing \\ 4 -> 1,3 \\ \\ The method is a bit twiddle on the index n. This suits large and \\ small n. default(strictargs,1); default(recover,0); \\ numbers not power of 2, starting A057716(0) = 0 A057716_not_power_of_2(n) = { n++; my(k=logint(n,2)); n += k; n + bittest(n,k+1); } A105316_table = Vecsmall([ 3,4, 1,3, 2,3 ]); \\ for n>=1 A105316(n) = { my(t=A057716_not_power_of_2(n), h=hammingweight(t), r=bitand(t,1)); if(h==2 && r, 1, \\ exception for t=2^k+1 A105316_table[(2*h-r)%6 + 1]); } A105316_S_term_table = Vecsmall([ 1,3, 2,3, 3,4 ]); \\ A105316_S_term(n) for n>=1 returns term n in an S() part. \\ Each S(k) is a prefix of S(k+1), so term n in whatever S(k) \\ is long enough to have n terms. A105316_S_term(n) = { n>=1 || error(); if(n==1, 1, A105316_S_term_table[(2*hammingweight(n) - bitand(n,1))%6 +1]); } { \\ printed with spaces to see each concatenated S() print1("A105316 = "); my(k=1,space_at=1); for(n=1,15, print1(A105316(n),","); if(n==space_at, print1(" "); space_at += 2^(k++)-1)); print("..."); print1("and each part S() = "); for(n=1,7, print1(A105316_S_term(n),",")); print("..."); } \\----------------------------------------------------------------------------- \\ S(k) Parts \\ ---------- \\ \\ expand(v) in the definition makes pairs [1,3] [2,3] [3,4]. \\ Since 3 expands to nothing, one such pair expands to another \\ such pair. \\ \\ Each pair in S(k-1) expands to a new pair which is appended to \\ make S(k), except the initial "1" (from S(0) = [1]) which remains \\ unpaired and so \\ \\ length(S(k)) = 2^k - 1 \\ \\ Pairs are distinguished by which of 1,2,4 they contain. \\ These 1,2,4 expand in a cycle \\ \\ 1 -> 2 2 -> 4 4 -> 1 \\ \\ In S(k) use an index i=0 for the initial [1] and then indices \\ 1 <= i <= 2^(k-1) - 1 for the pairs. For example \\ \\ k=4 \\ S(k) = 1, 2,3, 2,3, 3,4, 2,3, 3,4, 3,4, 1,3 \\ i = 0 1 2 3 4 5 6 7 \\ \\ expand() means new pairs appended in S(k) are taken from S(k-1) by \\ \\ pair(i + 2^(k-1)) = expand(pair(i)) \\ \\ So the number of expand()s which result in pair i is \\ hammingweight(i). The cycle 1->2->4->1 in expand() means take \\ that hammingweight mod 3, \\ \\ pair(i) = [1,3] or [2,3] or [3,4] \\ according as hammingweight(i) == 0,1,2 mod 3 respectively \\ for i>=1 \\ \\ For A105316_S_term(n), indexed n>=1, the pair terms are by an \\ index into A105316_S_term_table[] \\ \\ i = floor(n/2) pair number \\ r = n mod 2 whether take the first or second in pair \\ = bitand(n,1) \\ table index = (2*hammingweight(i) + r) mod 6 + 1 \\ for n>=2, which is i>=1 \\ \\ The implementation uses hammingweight(n) instead of hammingweight(i), \\ since that's parameter n unchanged and is as easy as "-r" instead \\ of "+r" for the table index, \\ \\ hammingweight(i) = hammingweight(n) - r \\ table index = (2*hammingweight(n) - r) mod 6 + 1 \\ Whole Sequence \\ -------------- \\ \\ The whole A105316 sequence is concatenation of S(1),S(2),etc. \\ The index desired in those S() parts runs 1 .. 2^k-1 for \\ successive k. In terms of A-numbers, this is \\ \\ A105316(n) = A105316_S_term(A115218(n) + 1) \\ where A115218 = successive runs 0 .. 2^k-2 \\ \\ This kind of runs of length 2^k-1 means omit 1 element from each \\ 2^k block. A105316(n) is implemented by going to \\ \\ t = A057716_not_2pow(n) \\ numbers not a power of 2 \\ \\ Omitting powers of 2 means omit 1 from each block of 2^k. \\ The method in A057716_not_2pow() is a fairly usual \\ \\ A057716(n) = m+k if m+k < 2^(k+1) \\ = m+k+1 if m+k >= 2^(k+1) \\ where m = n+1, k = floor(log2(m)) \\ \\ If t = m+k >= 2^(k+1) then t has stepped over an additional \\ power of 2 (ie. 2^(k+1)) and so has to be 1 bigger. \\ Comparison >= 2^(k+1) is a bittest() since "+k" in m+k at most \\ pushes up into one additional bit position. \\ \\ t = A057716(n) has an additional most significant 1 bit over \\ the desired 1..2^k-1 runs. That would be removed to index into \\ an S(k) part, \\ \\ A105316(n) = A105316_S_term(t - 2^k) \\ where t = A057716_not_2pow(n) \\ k = floor(log2(t)) \\ so 1 <= t - 2^k <= 2^k - 1 \\ \\ A105316() code instead leaves the 2^k since it's a single \\ additional 1 bit and so merely increases hammingweight(t) by 1 \\ and that can be accounted for in A105316_table[]. \\ \\ The table index uses hammingweight(t) in the same way as \\ described above for A105316_S_term() (rather than floor(t/2)). \\ \\ The first term in each S(k) is an exception (not in a pair), \\ which is the n=1 case in A105316_S_term(). In A105316() that \\ becomes each t = 2^k + 1 and they're identified by h=2 and r=1. \\ Morphism \\ -------- \\ \\ Some early comments in A105316 (since dropped for unclearness) \\ suggested an aim was 4 symbols where the expansions include \\ each from -> to combination in one direction and exactly once, \\ and one sink symbol (symbol 3). It's not clear what merit \\ concatenating S(k) parts has, since each S is just a prefix of \\ the next. \\ \\ Both S(infinity) and the resulting whole A105316 are morphic, \\ meaning they're symbol to symbol mapping of a morphism. \\ \\ One such morphism is by introducing two additional types of "1". \\ \\ 1a -> 1a, 1b,2,3 variant 1 start of whole sequence \\ 1b -> 1b, 2,3 variant 1 start of S(k) part \\ 1 -> 1,3, 2,3 \\ 2 -> 2,3, 3,4 \\ 3 -> nothing \\ 4 -> 3,4, 1,3 \\ \\ A105316 = fixed point starting from 1a, \\ mapped down 1a -> 1 and 1b -> 1 \\ \\ S(infinity) = fixed point starting from 1b, \\ mapped down 1b -> 1 \\ \\ In this morphism, 1,2,4 expand the way hammingweight mod 3 does, \\ \\ A071858 = hammingweight mod 3 \\ 0 -> 0,1 morphism starting from 0 \\ 1 -> 1,2 \\ 2 -> 2,0 \\ \\ In a uniform morphism like this (each expansion the same length), \\ the first and second new terms are a new least significant 0 or 1 \\ bit (respectively) on an index >=0. For hammingweight mod 3, \\ the new terms are resulting new 1 bit counts mod 3. \\ \\ In the A105316 morphism, 1,2,4 are the corresponding two \\ hammingweight terms and an additional 3 each to make pairs. \\ \\ Symbol 1a for the whole sequence is the first part S(1)=[1]. \\ It expands to S(2) and preserves itself ready to do the same \\ again on next expansion. \\ \\ Symbol 1b is the initial 1 at the start of each S part. \\ It differs from an ordinary 1 by not being paired so it expands \\ to preserve itself as 1b rather than pair 1,3 (and still followed \\ by 2,3 the same as an ordinary 1). \\----------------------------------------------------------------------------- print("end");