\\ Kevin Ryde, April 2023 \\ This is some code to calculate an individual term of A359368, \\ which is a certain recurrence taking a(n-1)+a(n+1) or -a(n). \\ The strategy is double-and-step on a block of 5 values. \\ This suits large and small n. default(recover,0); default(strictargs,1); my(initial=[1, 1, 1, 2, -1, 3, -1, 0, -2]); \ A359368(n) = { if(n<=9, initial[n], my(w=logint(n,2), h=n>>(w-2), \\ high three bits, value 4..7 v=initial[h-2 .. h+2]); forstep(i=w-3,0,-1, v=if(bittest(n,i), [ -v[2], v[2]+v[4], -v[3], v[3]+v[5], -v[4]], [v[1]+v[3], -v[2], v[2]+v[4], -v[3], v[3]+v[5]])); v[3]); } { print1("A359368 = "); for(n=1,16, print1(A359368(n),", ")); print("..."); } \\----------------------------------------------------------------------------- \\ Implementation Notes \\ A359368 is defined by \\ \\ a(1..3) = 1 and thereafter \\ a(2n) = a(n-1) + a(n+1) \\ a(2n+1) = - a(n) \\ \\ A block of 5 consecutive values a(n-2 .. n+2) suffice to determine \\ the corresponding 5 consecutive values surrounding either 2n or 2n+1. \\ \\ around 2n around 2n+1 \\ 2n-2..2n+2 2n-1..2n+3 \\ a(2n-2) = a(n-2) + a(n) \ \\ a(2n-1) = - a(n-1) | \ \\ a(2n) = a(n-1) + a(n+1) | | \\ a(2n+1) = - a(n) | | \\ a(2n+2) = a(n) + a(n+2) / | \\ a(2n+3) = - a(n+1) / \\ \\ In the code, vector v = [a(n-2 .. n+2)] for a working "n" and a new least \\ significant bit uses the recurrence cases above, \\ \\ new bit 0 to new vector [a(2n-2..2n+2)] \\ new bit 1 to new vector [a(2n-1..2n+3)] \\----------------------------------------------------------------------------- \\ Matrix Powers \\ \\ The recurrence cases are linear combinations of terms from the "v" block. \\ They can be expressed by matrix multiplies on the right of a column \\ vector v. M0 = {[1, 0, 1, 0, 0; 0, -1, 0, 0, 0; 0, 1, 0, 1, 0; 0, 0, -1, 0, 0; 0, 0, 1, 0, 1]}; M1 = {[0, -1, 0, 0, 0; 0, 1, 0, 1, 0; 0, 0, -1, 0, 0; 0, 0, 1, 0, 1; 0, 0, 0, -1, 0]}; { my(column_of_5(n) = apply(A359368, [n-2 .. n+2]~)); for(n=4,100, M0 * column_of_5(n) == column_of_5(2*n) || error(); M1 * column_of_5(n) == column_of_5(2*n+1) || error()); } \\ It can be noted that powers of M0 and M1 satisfy \\ \\ M0^4 = id the 5x5 identity matrix \\ M1^5 = M1 M0^4 == matid(5) || error(); M1^5 == M1 || error(); \\ So any run of four 0-bits in n can be discarded, and any run \\ of five 1-bits in n can be collapsed to a single 1-bit. \\ \\ In program code, reducing run lengths mod 4 could be good \\ if n is already in run lengths form, or have a good bit-scan \\ operation, but otherwise those runs may be too unlikely \\ (in random input) to be worth watching for. \\ \\ When wanting just an individual A359368(n), like the A359368() \\ function above does, the defining recurrence a(2n+1) = -a(n) \\ means a least significant final run of 1-bits in n are just \\ sign changes. In the M1 matrix powers this is row alternating, M1 [3,] == [0,0,-1,0,0] || error(); (M1^2)[3,] == [0,0, 1,0,0] || error(); (M1^3)[3,] == [0,0,-1,0,0] || error(); (M1^4)[3,] == [0,0, 1,0,0] || error(); \\ and M1^5 == M1 back to M1 again \\----------------------------------------------------------------------------- print("end");