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