Central moments of the Weibull distribution are given by combinations of gamma functions, often involving strong cancelations. For example, the variance is given by
$V = \Gamma(1 + 2/k) - [\Gamma(1+1/k)]^2$
For large k, both these values are close to $\Gamma(1)$ so there is siginificant cancelation. Of course there is a standard way to deal with this: expand each term in a series in $1/k$, see that the leading terms cancel, and evaluate the rest. And that works here, kind of. But there are some problems.
The first problem is that in this case the expansion is unusually laborious: (i) Expand $\log \Gamma$ near one in a Taylor series with coefficients $\psi^{(n)}(1)$, i.e. Euler's $\gamma$ and $\zeta(n)$. (ii) Use the exponential formula to transform that into a series of $\Gamma$ itself. This mixes up the coefficients that appear at order $n$ into a series of combinations of the $\psi^{(n)}(1)$ corresponding to the integer partitions of $n$. (iii) Transform that series into a series for this square, which of course does more coefficient mixing. In the end, the first two terms cancel and you end up with
$V = \frac{\zeta(2)}{k^2} - \frac{2[\gamma \zeta(2) + \zeta(3)]}{k^3} + \mathcal{O}(\frac{1}{k^4})$
Higher terms getter more complicated fast, but you can use Mathematica to get them, an numerically the result looks like
$V = (0.4112)(2/k)^2 - (0.5379)(2/k)^3 + (0.7324)(2/k)^4 + (0.8291)(2/k)^5 - (0.9012)(2/k)^6 + (0.9425)(2/k)^7 - (0.0.9674)(2/k)^8 + \cdots$
The next problem is, once you have this series, you find that it isn't all that good. Already at k=4 the difference of Gamma functions looses 4 digits of precision to cancelation, so you want to be switching to the series development around there. But that is way to early for this series: even for k=16, you need 16 terms to achieve full floating-point precision. Series acceleration techniques buy you a little, but not nearly enough to get you down to k=4.
So I'm looking for another approach to this problem. Is there some way to express differences of Gamma functions like this via some other function? Some other kind of expansion, e.g. a continued fraction? Is there some hint at a solution from the fact that coefficients in the series expansion eventually approach the values $\pm 1$ of the series expansion of $1/(1+x)$?
(By the way, the same problem occurs for higher central moments, and the problem gets even worse. In the expression for the third central moment, the first three terms in the even more complicated series expansion cancel, indicating that cancelations are even stronger so one wants to use the series even sooner, but the convergence properties of that series look even worse.)
