\\ Kevin Ryde, February 2024

\\ This is some PARI/GP code to calculate an individual term of A116529,
\\ which is a recurrence
\\
\\     a(2*n+1) = a(n)
\\     a(2*n+2) = 2*a(n) + a(n-1)
\\     starting a(0..2) = 1, 1, 2
\\
\\ The method is double-and-step on a block of 3 values.
\\ This suits large and small n.

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


\\ A116529_triple(n) returns a vector of 3 values
\\     [ a(n-3), a(n-2), a(n-1) ]
\\ Must have n >= 0.
\\
\\ For n=0..2, imagined terms a(-3..-1) = -1,1,0 show in the return.
\\ They're not part of the sequence: only an artifact of the
\\ implementation.
\\
A116529_triple(n) =
{
  my(k = if(n,logint(n,2),-1),
     v = [-1, 1, 0]);
  forstep(i=k,0,-1,
     if(bittest(n,i),
        v = [      v[1]+2*v[2], v[3], v[2]+2*v[3]],
        v = [v[2], v[1]+2*v[2], v[3]             ]));
  v;
}

A116529(n) =
{
  n >>= valuation(n+1,2);   \\ strip low 1 bits
  A116529_triple(n+1)[3];
}

{
  print1("A116529 = ");
  for(n=1,16, print1(A116529(n),", "));
  print("...");
}


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

\\ Blocks of 3 Terms
\\ -----------------
\\
\\ The recurrence cases can be written
\\
\\     a(2n-3) =             a(n-2)
\\     a(2n-2) =  a(n-3) + 2*a(n-2)
\\     a(2n-1) =                        a(n-1)
\\     a(2n)   =             a(n-2) + 2*a(n-1)
\\
\\ Knowing 3 values a(n-3),a(n-2),a(n-1) allows a corresponding
\\ 3 values to be calculated at either 2n or 2n+1,
\\
\\     n-3,n-2,n-1  ->  2n-3,2n-2,2n-1
\\                   or      2n-2,2n-1,2n
\\
\\ So bits of a target n can be considered from most to least
\\ significant, successively doubling a working index n to either
\\ 2n or 2n+1 according as the bits of the target n.
\\
\\ Double and step is a usual approach and in a recurrence of this
\\ type the only thing needed is to identify what set of terms can
\\ double out to 2n and 2n+1.


\\ Initial Values
\\ --------------
\\
\\ The starting 3 values are
\\
\\     starting n=0, no bits
\\     using a(-3) = -1
\\           a(-2) = 1
\\           a(-1) = 0
\\
\\ a(-1) and a(-2) are by extending the recurrence cases back.
\\ a(-3) is not uniquely determined, since even and odd cases would be
\\
\\     a(-2) = 2*a(-2) + a(-3)   even case, must a(-3) = -1
\\     a(-3) = a(-2)             odd case,  a(-3) = 1
\\
\\ But A116529_triple() only uses a(-3) for the even case going
\\ n=0 to n=1, hence take a(-3) = -1.


\\ Trailing 1 Bits
\\ ---------------
\\
\\ The odd case a(2*n+1) = a(n) allows any least significant 1 bits
\\ to be discarded from n since they're no change to the result.
\\ A116529() does this as a small optimisation.
\\
\\     a(n) = a(n sans trailing 1 bits)
\\
\\ This includes n = 2^k-1 which is all 1 bits,
\\
\\     a(2^k-1) = a(0) = 1


\\ Matrix Multiply
\\ ---------------
\\
\\ The double or double and step cases in A116529_triple() are
\\ linear combinations of the v[] terms and can be expressed by
\\ matrix multiplies.
\\
\\ The following M0 and M1 are for 0 bit or 1 bit multiplying on
\\ the left of a column vector v = [a(n-3),a(n-2),a(n-1)]~.
{
  M0 = [0,1,0;
        1,2,0;
        0,0,1];
  M1 = [1,2,0;
        0,0,1;
        0,1,2];
}
A116529_triple_by_matrices(n) =
{
  \\ Bits of n from least to most significant
  \\ become product M0 or M1 left to right.
  my(b=Vecrev(binary(n)));
  Vec(prod(i=1,#b, if(b[i]==0,M0,M1)) * [-1, 1, 0]~);
}
{
  for(n=0,256,
     A116529_triple(n) == A116529_triple_by_matrices(n)  || error());
}

\\ M0 and M1 have the same characteristic polynomial,
\\
charpoly(M0) == (x-1) * (x^2-2*x-1)   || error();
charpoly(M1) == (x-1) * (x^2-2*x-1)   || error();

\\ The largest eigenvalue is 1+sqrt(2) = 2.414... which allows
\\ terms to grow by that factor per bit of n.  For example,
\\
\\     a(2^k-2) = Pell(k)
\\     where Pell(k) = A000129(k) is the Pell numbers,
\\     which grow as (1+sqrt(2))^k

my(c=1+quadgen(8)); /* 1+sqrt(2) */ \
Pell(n) = imag(c^n);

for(k=1,100, A116529(2^k-2) == Pell(k) || error());

\\ It seems record high values might be those terms and a similar
\\ one in between,
\\
\\     a(2*2^k-2) = Pell(k+1)
\\     a(3*2^k-2) = Pell(k+2) - Pell(k+1)


\\ Matrix Powers
\\ -------------
\\
\\ n with a repeating bit pattern can usually be treated easily since
\\ matrix powers like M0^k or M1^k have matrix elements which are
\\ linear recurrences.  The powers of M0 and M1 are
{
  for(k=0,100,
    M0^k == [Pell(k-1), Pell(k),   0;
             Pell(k),   Pell(k+1), 0;
             0,         0,         1]  || error();

    M1^k == [1,  Pell(k)-Pell(k-1)+1, Pell(k+1)-Pell(k)-1;
             0,  Pell(k-1),           Pell(k);
             0,  Pell(k),             Pell(k+1)]  || error());
}
\\
\\ In M0^k, the 3rd row 0,0,1 is the odd recurrence case which is no
\\ change,
\\
\\      a(2*n-1) = a(n-1)
\\      M0 maps n-3,n-2,n-1  ->  2*n-3,2*n-2,2*n-1
\\      3rd row 0,0,1 so 3rd entry of column vector v unchanged
\\

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