Sat Jun 22 08:08:08 2013, comment #2:
I increased the limit on the number of iterations to a very large number, 1000000000. That allows using larger numbers as inputs:
but notice that a different type of error is beginning to surface. Also, near that point of failure the answer is starting to move away from 0.5.
ans = 0.39347
Without doing a mathematical convergence proof, my gut feeling is that the CDF for a Poisson r.v. when k = lambda should tend to 0.5 for large lambda. The Poisson r.v. begins to approximate a Gaussian r.v. with mean lambda and variance lambda.
So, at some point for large lambda this routine is experiencing numerical problems. Again, it isn't surprising because of the way gamma is defined. If one were to use the definition of CDF (accumulating PDF values), the probability would be spread so thin that I'm guessing that too would have numerical problems as well.
I think what we can do is increase the limit on that DO loop some, but at some point the user is going to have to realize there has to be reasonable bound on those input values or switch over to an approximation such as a Gaussian r.v., or scale their problem somehow.

Sat Jun 22 04:03:42 2013, comment #1:
Just as an initial investigation, I increased the allowable number of iterations of the core loop in the following way:
and then found that the largest value for similar integer inputs is 4241, i.e.:
Convergence obviously slows down with increasing magnitude of input, but running "poisscdf (4241,4241)" doesn't seem particularly longer. I'd say there might be an artificial limit on the number of iterations.
We could increase the allowable number of iterations. Typically, whenever in numerical routines I've experienced nonconvergence, the algorithm runs for a few seconds. So, I see no problems with extending the limit to be much higher if the algorithm doesn't take but a fraction of a second. That is sort of the simplest approach without having to do much analysis. I just increased the limit to 5000 and can get up to:
[We can see this is converging to 0.5, probably converging to some other distribution as argument tends to infinity.]
On the other hand, a more analytic approach would be to recognize some characteristics for the input in which the computation should be translated into a different problem which has faster convergence.
BUT, before doing any of that, I'm wondering if slatec:
http://www.netlib.org/slatec/index.html
is actively maintained. It doesn't appear to be. If a patch is made to the routine, it should probably be renamed to avoid confusion, or maybe even written is C.
One last thing is that the same problem exists in d9lgic.f, which apparently is ultimately called if the second input is larger than the first. (The "c" means complimentary.):
I will think this over a bit. I'm staring at the integral equation for the incomplete gamma function at the moment and don't see an easy way to translate for very large numbers. I can see why the equation becomes so sensitive. t^(largenumber  1) grows very fast. There is the normal gamma recursion translation, but that still doesn't get one very far.

Fri Jun 21 11:28:33 2013, original submission:
Hi,
this is an old bug report from Debian that I forward here in case someone wants to spend some time on Fortran or GSL. Please note that John already commented in the original bug report years ago.
http://bugs.debian.org/cgibin/bugreport.cgi?bug=418158
So, the issue is that
poisscdf (604,604)
works, but
poisscdf (605,605)
fails with a Fortran error, as the algorithm did not converge. This is ultimately due to the algorithm in libcruft/slatecfn/d9lgit.f.
There was the suggestion to use the gamma_inc_P function from the gsl package, but as far as I can see, this was never implemented.
Regards
Thomas
