\\ Kevin Ryde, August 2022 \\ \\ This is some PARI/GP code to calculate terms of A007843, the least k for \\ which k! is divisible by 2^n, meaning k! ends with n or more 0 bits. \\ \\ a(n) is an individual term by an attractive bit twiddle which deserves \\ some explanation as to why it works. a_vector(len) is repetitions per \\ Amarnath Murthy's comment in the sequence. default(strictargs,1); default(recover,0); a(n) = { if(n==0,1, forstep(i=logint(n,2),0,-1, n += bittest(n,i)); n); } { print1("A007843 = "); for(n=0,15, print1(a(n),", ")); print("..."); } \\ a_vector() returns a vector of the first "len" many terms. \\ The first term in the vector is a(0)=1. \\ a_vector(len) = { my(v=vector(len)); if(len, v[1]=1; my(n=1); \\ last stored to v[n] forstep(k=2,oo,2, for(i=1,valuation(k,2), \\ run of terms if(n>=len, break(2)); v[n++] = k))); v; } a_vector(1024) == vector(1024,n,n--; a(n)) || error(); \\ Runs of Terms \\ ------------- \\ \\ a_vector() is per comment by Amarnath Murthy in the sequence. \\ In terms of the sequence definition, k! has an additional valuation(k,2) \\ many low 0-bits over what (k-1)! had. So k is used for the next \\ valuation(k,2) many n. \\ \\ The initial a(1)=1 is an exception to terms otherwise all even. \\ Variations on the loops are possible, including for instance an \\ expression within vector(). \\ \\ \\ Recurrence \\ ---------- \\ \\ The a(n) code is a simple bit twiddle. It uses \\ \\ a(2^i + r) = / 2^i + a(r+1) for 0 <= r <= 2^i-2 \\ \ 2^i + 2^i for r = 2^i-1 { for(i=0,10, for(r=0,2^i-2, a(2^i + r) == 2^i + a(r+1) || error()); a(2^i + 2^i-1) == 2^i + 2^i || error()); } \\ \\ Do these formulas follow from some of the meta-Fibonacci considerations \\ in the related A046699? In any case they can be had here by considering \\ how runs of terms repeat. Runs can be illustrated, \\ \\ 2^i \\ n = 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 \\ a(n) = 1, 2, 4,4, 6, 8,8,8, 10, 12,12, 14, 16,16,16,16, \\ reps 1 2 1 3 1 2 1 4 \\ \--------------/ \---------------------/ \\ \\ Terms a(1..7) = 2 .. 8 are followed by terms 10 onwards. \\ Each value 10..14 = 8 + 2..6 which is an extra high 1-bit and which is no \\ change to their 2-adic valuation. So the run lengths of terms 10..14 are \\ the same as the run lengths of terms 2..6. \\ \\ Terms a(5..7) = 8 are run length 3, and a(12..14) = 8+8 = 16 are per the \\ recurrence 2^i + a(r+1) too, which completes the 0 <= r < 2^i-1 case. \\ \\ a(15) = 16 is the extra a(2^i + r) = 2^i + 2^i when r = 2^i-1. \\ \\ \\ Code a(n) \\ --------- \\ \\ In the a(n) code, the working n holds result bits so far at the high end \\ and a reduced index at the low end. \\ \\ high i low \\ n = binary [result] b [remaining r] \\ \-------------/ \\ to become a(b*2^i + r) \\ \\ The result bits are the "2^i" parts of "a(etc) = 2^i + etc". \\ \\ At the start of loop i, the aim is to calculate a(b*2^i + r). \\ If b=0 then no action, wait until reach b=1. \\ \\ If b=1 and if r < 2^i-1, then recurrence a(2^i + r) = 2^i + a(r+1) means \\ b=1 at i is unchanged to go into the result, and the remaining r part is \\ incremented. This increment does not cause a carry into bit position i \\ since r < 2^i-1. \\ \\ If b=1 and if r = 2^i-1, meaning r is all 1-bits, then increment gives \\ \\ a(2^i + 2^i-1) = 2^i + 2^i-1 + 1 \\ = 2^i + 2^i per formula \\ \\ This leaves all bit positions < i as 0 bits, so the rest of the i loop is \\ no change from n+=bittest, again as desired for the formula. \\ \\ \\ Initial a(0) \\ ------------ \\ \\ n=0 is a special case in the a(n) code since the loop would not produce \\ a(0)=1. It would make 0 by no iterations, and actually that could have \\ suited the definition if 0! = 1 were taken for n=0 many low 0-bits. \\ In any case, a(0) doesn't affect the recurrence so its value doesn't \\ matter. \\----------------------------------------------------------------------------- print("end");