\\ Kevin Ryde, May 2024
\\
\\ This is some PARI/GP code to calculate an individual term of
\\ 
\\     A369989 = infinite word "y" of Berstel et al
\\     
\\ The method is my formula which is a bit twiddle on quotient
\\ and remainder of index n>=0 divided by 9.

default(recover,0);
default(strictargs,1);

A369989(n) =
{
  my(q,r); [q,r]=divrem(n,9);
  if(r==2||r==4,
     bitxor(bittest(r,1), bittest(q,valuation(q+1,2)+1)),
     bittest(135369, r + if(q%2==1,9)));
}

{
  \\ spaces between "1,1" to show blocks of expansion h
  my(a=A369989);
  print1("A369989 = ");
  for(n=0,24,
     if(n>=1 && a(n)==1 && a(n-1)==1, print1(" "));
     print1(a(n),","));
  print(" ...");
  print();
}


\\-----------------------------------------------------------------------------
\\ The h expansion which is the definition of A369989 is

A369989_h_table={[[1,0,1],
                  [1,0,0,1],
                  [1,0,0,0,1],
                  [1,0,0,0,0,1]]};

\\ A369989_h_expand(v) takes v a vector of terms 0..3 as from A372257.
\\ Return them expanded to runs 1,0,...,0,1 per the h map.
A369989_h_expand(v) = concat(apply(t->A369989_h_table[t+1], v));

A372257(n) =
{
  bitxor(bitand(n,3),
         bittest(n, valuation(n>>1+1,2)+2));
}

{
  \\ check A369989() vs expand
  my(v=vector(100,n,n--; A372257(n)));
  v = A369989_h_expand(v);
  vector(#v,n,n--; A369989(n)) == v  || error();
}


\\-----------------------------------------------------------------------------
\\ As described in the sequence, blocks of 4 terms of A372257 become
\\ 18 terms here.

{
  print("Blocks of 4 terms from A372257 expand under h");
  \\ it happens that i=4..7 are the 4 possible blocks which occur
  for(i=4,7,
     my(v=apply(A372257, [4*i..4*i+3]),
        e=A369989_h_expand(v));
     print("  "v" -> "e));
}

\\ Notice that in the blocks of 4 expanded to blocks of 18,
\\ only 4 of the terms vary with remainder s = n mod 18.
\\ The rest are fixed.
\\
\\     s =    2     4                   11    13
\\     [1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1]
\\     [1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1]
\\     [1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1]
\\     [1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1]
\\      ^        ^        ^  ^              ^              ^
\\     fixed 1s
\\
\\ The magic number 135369 in the A369989() code has 1 bits at the
\\ fixed 1 terms,
\\ 
\\     Vecrev(binary(135369)) == \
\\     [1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1]
\\      ^        ^        ^  ^              ^              ^
\\
\\ The sequence formulas and A369989() code above use a quotient and
\\ remainder mod 9 since it can be noticed that the cases in the
\\ underlying A372257 mean that remainders s=2 and s=11 can be
\\ treated together, and similarly s=4 and s=13 together.


\\-----------------------------------------------------------------------------
print("end");