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.
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?
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.