\\ Kevin Ryde, November 2021 \\ \\ Usage: gp <a347204.gp \\ \\ Calculate an individual term of A347204 from its index n using products \\ of 1-bit positions. This method suits large and small n and the work \\ done grows roughly as the number of 1-bits in n and to a lesser extent \\ the total bit length of n. default(recover,0); default(strictargs,1); a(n) = { my(s=Vec([1],hammingweight(n)+1), \\ s[1]=1 then s[2..h+1] = 0 end=0); for(i=0,if(n,logint(n,2)), if(bittest(n,i), forstep(j=end++,1,-1, s[j+1] += (i-j+2)*s[j]))); vecsum(s); } { print1("A347204 = "); for(n=0,17, print1(a(n),",")); print("..."); } \\ This code is based on the following formula from A347204 (among various \\ equivalents there). \\ \\ a(n) = a(C) + V*a(C/2) for n>=1, and a(0)=1 \\ where \\ C = ClearLowOne(n) = A129760(n) = n with lowest 1-bit cleared to 0 \\ V = LowOnePosition_1based(n) = A001511(n) = bit position of lowest \\ 1-bit in n, where position 1 is the least significant LowOnePosition_1based(n) = valuation(n,2) + 1; ClearLowOne(n) = n - 1 << valuation(n,2); { for(n=1,4096, my(C = ClearLowOne(n), V = LowOnePosition_1based(n)); a(n) == a(C) + V*a(C/2) || error()); } \\ This formula has a(C) and a(C/2) and they can be calculated together. \\ The 1-bits in C/2 are the same as the 1-bits of C except one bit position \\ lower, so the next 1-bit position V' in C will be V'-1 in C/2, or V'-2 in \\ C/4, and so on. \\ \\ Repeatedly expanding the formula gives products which can be illustrated \\ with a few 1-bits, \\ \\ n = 2^(b-1) + 2^(c-1) + 2^(d-1) + 2^(e-1) \\ ascending bit positions 1 <= b < c < d < e \\ \\ a(n) = 1 \\ + b + c + d + e \\ + b*(c-1) + b*(d-1) + b*(e-1) + c*(d-1) + c*(e-1) + d*(e-1) \\ + b*(c-1)*(d-2) + b*(c-1)*(e-2) + b*(d-1)*(e-2) + c*(d-1)*(e-2) \\ + b*(c-1)*(d-2)*(e-3) \\ \\ The two-terms products are all pairs of b,c,d,e. In each pair, the \\ second one (the larger) has "-1", eg. b*(c-1). \\ \\ The three-terms products are all triples chosen from from b,c,d,e. \\ In each triple, the second has a "-1" and the third has a "-2", again \\ in ascending order of the b,c,d,e values. \\ { \\ four b,c,d,e ascending my(limit=25); for(b=1,limit, for(c=b+1,limit, for(d=c+1,limit, for(e=d+1,limit, my(n=2^(b-1) + 2^(c-1) + 2^(d-1) + 2^(e-1)); a(n) == 1 + b + c + d + e + b*(c-1) + b*(d-1) + b*(e-1) + c*(d-1) + c*(e-1) + d*(e-1) + b*(c-1)*(d-2) + b*(c-1)*(e-2) + b*(d-1)*(e-2) + c*(d-1)*(e-2) + b*(c-1)*(d-2)*(e-3) || error())))); } \\ This kind of product can be taken on any vector of values b,c,etc. \\ In a(n), the terms are in ascending order but when just making products \\ there's no restrictions. \\ decprod(v) = { my(s=Vec([1],#v+1)); \\ s[1]=1 then s[2..#v+1] = 0 for(i=1,#v, forstep(j=i+1,2,-1, my(t = v[i] - (j-2)); s[j] += t*s[j-1])); vecsum(s); } \\ decprod() doesn't enquire into the nature of its v[] elements, \\ allowing the above b,c,d,e formula to be checked symbolically. { decprod([b,c,d,e]) == 1 + b + c + d + e + b*(c-1) + b*(d-1) + b*(e-1) + c*(d-1) + c*(e-1) + d*(e-1) + b*(c-1)*(d-2) + b*(c-1)*(e-2) + b*(d-1)*(e-2) + c*(d-1)*(e-2) + b*(c-1)*(d-2)*(e-3) || error(); } \\ return a vector which is the positions of the 1-bits within n, \\ in ascending order and least significant bit as position 1 n_1bit_positions_1based(n) = \ select(i->bittest(n,i-1), [1..if(n,logint(n,2)+1)]); { \\ a(n) using decprod() n_1bit_positions_1based(9) == [1,4] || error(); for(n=0,4096, a(n) == decprod(n_1bit_positions_1based(n)) || error()); } \\ In decprod() and a(), s[j] is the sum of products of j-1 many terms \\ (seen so far), \\ \\ s[1] = sum of products of no terms, so always 1 for empty product \\ \\ s[2] = sum of products of one term, so simply sum of terms so far \\ \\ s[3] = sum of products of two terms, such as b*(c-1), in all \\ combinations so far and with "-1" on the second (later position \\ in v[], bigger value in a()) \\ \\ etc \\ \\ With terms b,c so far, the way new term d accumulates into s[] can be \\ illustrated, \\ \\ existing unchanged additional \\ s[1] = 1 ==> s[1] = 1 \\ s[2] = b + c s[2] = b + c + ( 1 ) * d \\ s[3] = b*(c-1) s[3] = b*(c-1) + ( b + c ) * (d-1) \\ s[4] = ( b*(c-1) ) * (d-2) \\ \\ d or d-1 or d-2 is a new term on the products, so the preceding s[j-1] \\ products multiplied by d make products with one more term and so add \\ to s[j]. The existing products in s[j] remain, being combinations \\ without d. The d is d or d-1 or d-2 etc decreasing according to how many \\ existing terms it's multiplying onto, \\ \\ new s[j] = s[j] + t*s[j-1] where t = d - (j-2) \\ \\ s[] is modified in-place by working from largest j down so that s[j-1] is \\ not replaced until after s[j] has used its value. \\----------------------------------------------------------------------------- \\ Derivatives \\ \\ The factors of d, d-1, d-2 onto successive elements of s[] holding s[] as \\ coefficients of a polynomial of degree d and taking the derivative, \\ \\ d/dx x^d = d * x^(d-1) \\ d/dx x^(d-1) = (d-1) * x^(d-2) etc \\ \\ a_by_derivative() is a compact form for a(n) using this idea. it makes \\ good use of PARI primitives and runs a touch faster than the vector a() \\ presumably for that reason. However when n is very large, something \\ happens perhaps copying polynomials around so it becomes a touch slower \\ and uses more peak stack space. a_by_derivative(n) = { my(p=1); if(n, for(i=0,logint(n,2), p*=x; if(bittest(n,i), p+=p'))); vecsum(Vec(p)); } for(n=0,4096, a_by_derivative(n) == a(n) || error()); \\ The p+=p' steps can also be expressed by nested "plusderiv". \\ p*=x at each bit becomes a multiplier x^3 etc which is run length \\ of bits 100 etc in n, with the least significant run first \\ (deepest nested). \\ plusderiv(p) = p + deriv(p); { my(n = fromdigits([1,0,0, 1,0,0,0,0, 1,0,0,0], 2)); \\ \---/ \-------/ \-----/ \\ 3 5 4 run lengths a(n) == vecsum(Vec( plusderiv( x^3* plusderiv( x^5* plusderiv( x^4 ))))) || error(); } \\----------------------------------------------------------------------------- print("end");