[Rd] pgamma discontinuity (PR#7307)

[Peter Dalgaard]

terra at gnome.org writes:

>   for (x = 99990; x <= 100009; x++)
>     printf ("%.14g --> %.14g\n", x, pgamma (1000000.0,  x, 1.0, 0, 1));
....
> and have no further local changes.
> 
> I get...
> src/nmath/standalone> LD_LIBRARY_PATH=. ./test 
...
> 99999 --> -669752.66592471
> 100000 --> -669750.36332851
> 100001 --> -599731.36192347       <--- Look at me jump!
> 100002 --> -599729.89768519


> i.e., the result changes 70000 orders of magnitude right here.  Ugh.
> This is affecting some erlang calculations of mine pretty badly.

Make that 30400 orders of magnitude (natural logs y'know)... On
something thats about 300000 orders of magnitude below 1, mind you!
What the devil are you calculating? The probability that a random
configuration of atoms would make up the known universe?
 
> Judging from comments in pgamma someone else might have gotten hit by this
> around R1.8.1.

Beginning to look like a method bug. Recipe is to read the original
paper, check for transcription errors from the Fortran algorithm, then
check for errors going from theory to Fortran, then check the theory...

PS: It happens on the R command line too:

>  pgamma(1000000,100000.001,1.0,,0,1) - pgamma(1000000,100000,1.0,,0,1)
[1] 70017.54


-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907