\\ Kevin Ryde, August 2023
\\
\\ This is some PARI/GP code to calculate an individual term
\\ of A026465 which is the length of the n'th run in the
\\ Thue-Morse sequence (A010060).
\\
\\ The method is some bit twiddling on 3*(n-1) using the
\\ interpretation of A026465(n) as number of ways to write n-1
\\ as a sum of Jacobsthal numbers.  This suits large and small n.

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


A026465(n) =
{
  my(m=3*n-2, ret=1);
  forstep(k=bitnegimply(logint(m,2)+1,1), 2, -2,
     if(bittest(m,k),   ret=2, m--);    \\ k even
     if(bittest(m,k-1), ret=1, m++));   \\ k-1 odd
  ret;
}

{
  print1("A026465 = ");
  for(n=1,20, print1(A026465(n),","));
  print("...");
}


\\-----------------------------------------------------------------------------

\\ Representations
\\ ---------------
\\
\\ The Jacobsthal numbers are
\\
\\     J(k) = (2^k - (-1)^k)/3    (A001045)
\\          = 1, 1, 3, 5, 11, 21, 43, 85, ...  for k>=1
\\
\\ The representations interpretation of A026465 is
\\
\\     A026465(n) = number of ways to write n-1 as sum of J's,
\\     n-1 = J(x) + J(y) + ... with distinct indices x,y,...
\\
\\ For n=1, n-1 = 0 has 1 representation: the empty sum.
\\
\\ J(1)=1 and J(2)=1 are different indices and are different ways
\\ to have 1 in the sum.  For example,
\\
\\     n = 5
\\     n-1 = 4 = J(3) + J(2)  = 3 + 1   \ two representations
\\             = J(3) + J(1)  = 3 + 1   /
\\
\\ The greedy representation takes the largest possible J(k) into
\\ the sum each time.
\\ 
\\     for J(k) <= n-1 <= J(k+1)-1, greedy takes J(k)
\\     subtract to n'-1 = n-1 - J(k) and repeat
\\
\\  Every n-1 has such a representation since the largest remaining
\\  n'-1 can be reached using smaller terms
\\ 
\\     all(k-1) = J(k-1) + J(k-2) + ... + J(1)
\\
\\     0 <= n'-1 <= J(k+1) - J(k) - 1
\\                = all(k-1) - (1 if k odd)
\\               <= all(k-1)
\\
\\ A lazier representation omits some J(k) which greedy would take.
\\ This is only possible if n-1 (no J(k) subtraction) is small
\\ enough that all(k-1) can sum to it.
\\ 
\\     for J(k) <= n-1 <= J(k+1)-1, don't take J(k)
\\     requires n-1 <= all(k-1) = J(k) - (1 if k odd)
\\ so
\\     J(k) <= n-1 <= J(k) - (1 if k odd)
\\     
\\ which is possible only when n-1 = J(k) with k even.  Hence
\\
\\     A026465(n) = /  2   if greedy representation of n-1 ends k even
\\                  \  1   if greedy representation of n-1 ends k odd
\\
\\
\\ Implementation
\\ -------------- 
\\ 
\\ The code has k as index of a prospective J(k) to take into the
\\ greedy representation.  k is large enough that n-1 < J(k+1) so
\\ J(k+1) cannot be taken.  Whether n-1 is big enough to take J(k),
\\ ie. n-1 >= J(k), will be tested.
\\
\\ It's convenient to use a working quantity
\\
\\     m = 3*(n-1) + (1 if k even)
\\
\\ The test n-1 >= J(k) is then
\\
\\     m >= 2^k + (1 if k odd)
\\     m >= 2^k    since m multiple of 3 and 2^k is not
\\     hence bittest(m,k) in the code
\\
\\ If J(k) is taken into the sum then it's subtracted from n-1,
\\ which is subtract 3*J(k) from m.  m also switches "1 if k even"
\\ ready for k-1 next,
\\
\\     m' = /  m - 3*J(k) - 1  if k even
\\          \  m - 3*J(k) + 1  if k odd
\\        = m - 2^k  
\\
\\ This m - 2^k is not performed as such, since bit positions k
\\ and above are not examined again.
\\
\\ If J(k) is not taken into the sum, then switch "1 if k even",
\\
\\     m' = /  m - 1    if k even
\\          \  m + 1    if k odd
\\
\\ The code sets "ret" to 2 or 1 according to k even or odd at
\\ the smallest J(k) taken into the sum.  Initial value ret=1
\\ is for n-1 = 0 which takes no terms into the sum.
\\ 
\\ Stepping k by 2 and separate action for k even and k-1 odd is
\\ quite clean.  The same could be done with some arithmetic.


\\-----------------------------------------------------------------------------
\\ Generating Function

\\ A generating function product follows from the Jacobsthal
\\ interpretation in the usual way for how many ways to partition
\\ n into a sum of certain values, and which is per formula by
\\ Paul Barry in the sequence.
{
  \\ while() loop continues until Jacobsthal(k) > desired length,
  \\ since such k does not change g.  Offset so that A026465(n)
  \\ is representations of n-1 is factor x on g.
  \\
  my(Jacobsthal(k) = (2^k - (-1)^k)/3);
  my(maxn = 1000,
     g = x + O(x^(maxn+1)),
     k = 1);
  while(Jacobsthal(k) <= maxn, 
        g *= 1 + x^Jacobsthal(k);
        k++);
  g == sum(n=1,maxn, A026465(n)*x^n) || error();
}

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