(posting this here so the list can archive it, even though I've worked offline with Brandon on solving this)
Dear Brandon, Thank you for posting this bug. I looked in to it in detail and found that, as you had hinted with your guess, there are branching points that don't work when you have large similar values passed to gsl_cdf_gamma_Q(). I looked at the paper by Gil, Segura, and Temme: https://arxiv.org/abs/1306.1754 and by applying their criteria I was able to get the numbers in your sample program to give the correct incomplete gamma function. There is some detailed work before I can submit a patch to Patrick, but I'll give you a workaround for now.
