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