\\ 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");