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