[pontus]

Another way to get to the same result is to use "Feynman's Trick" of differentiating inside a sum:

Consider the function f(x) = Sum_{n=1}^\infty c^(-xn)

Then differentiate this k times. Each time you pull down a factor of n (as well as a log(c), but that's just a constant). So, the sum you're looking for is related to the kth derivative of this function.

Now, fortunately this function can be evaluated explicitly since it's just a geometric series: it's 1 / (c^x - 1) -- note that the sum starts at 1 and not 0. Then it's just a matter of calculating a bunch of derivatives of this function, keeping track of factors of log(c) etc. and then evaluating it at x = 1 at the very end. Very labor intensive, but (in my opinion) less mysterious than the approach shown here (although, of course the polylogarithm function is precisely this tower of derivatives for negative integer values).

[jwmerrill]

Instead of differentiating c^(-xn) w.r.t. x to pull down factors of n (and inconvenient logarithms of c), you can use (z d/dz) z^n = n z^n to pull down factors of n with no inconvenient logarithms. Then you can set z=1/2 at the end to get the desired summand here. This approach makes it more obvious that the answer will be rational.

This is effectively what OP does, but it is phrased there in terms of properties of the Li function, which makes it seem a little more exotic than thinking just in terms of differentiating power functions.

[lupire]

And since it's discrete, you can use finite differences.

  a =  sum_1 n^3 / 2^n 
    = sum_0 (n+1)^3 / 2^(n+1)
    =  (1/2) (1 + sum_1 (3n^2 + 3n + 1)/2^n)
    
  b = sum_1 n^2 / 2^n  
    = (1/2) (1+ b + sum_1 (2n +1)/2^n)

  c = sum_1 n / 2^n  
    = (1/2) (1+ c + sum_1  1 / 2^n)

  d =  sum_1 1 / 2^n  
    = sum_0 1 / 2^(n+1)
    = (1/2) (1 + d)
    = 1


  d = 1               =             1
  c = d +  1          = 1+1      =  2
  b = d + 2c +  1     = 1+1+4    =  6
  a = d + 3c + 3b + 1 = 1+1+6+18 = 26

In general,

  f_k = sum n^k/2^n 
      = (*k*th row of Pascal's triangle)•(f_0, ..., f_{k-1},1)

https://oeis.org/A000629

Number of necklaces of partitions of n+1 labeled beads.

[amelius]

To be honest, the whole use of the Li function before defining it made me stop reading.

[mturmon]

Yeah, differentiating these infinite sums to pull down polynomial factors is a familiar trick.

It happens in basic moment generating function manipulations (e.g., higher moments of random variables). Or from z-transforms in signal processing (z transforms of integrals or derivatives). And (a little less obvious, but still the same) from Fourier analysis.

The concept applies to any moment generating function, z-transform, whatever. It’s clearest for the geometric distribution, where the distribution itself has the geometric form (https://mathworld.wolfram.com/GeometricDistribution.html, around equation 6).

I agree that the Li function seems like a detour, but maybe it can make some of the manipulation easier?