Possible bug in POIDEV

[blackmime123]

I noticed this whilst using POIDEV function, of the F90 NR version.
POIDEV relies on a subroutine called assert 1, given in the module nrutil. As part of this, assert 1 tests whether the value input into POIDEV is greater than 0, a condition needed to fulfill the mathematical operation of POIDEV. However in it's current form, there is a single situation in which this can fail. this produces the error message:

nrerror: an assertion failed with this tag:gammln_s arg


The expression includes the calculation: y=tan(PI*harvest).
In my simulation I found that if harvest, a random number, is equal to 0.5 then y takes a value of infinite (as tan(pi/2) = infinite). this will pass the test as to whether the value is larger than 0, but when it is converted to a integer later it can be assigned a negative value, as it cannot cope with infinite.In this situation the programme will stop, according to assert1.
There is a simple solution to this, which avoids the issue, though I have not tested to see if there are any other implications. this is given below:

Current code:

call ran2(harvest)
y=tan(PI*harvest)
em=sq*y+xm
if (em >= 0.0) exit
end do
em=int(em)

Edit to:

call ran2(harvest)
y=tan(PI*harvest)
em=sq*y+xm
em=int(em)
if ((em >= 0.0)) exit
end do

By adjusting em to an integer inside the loop the problem is solved. I suppose another simple solution would be to add a condition to draw another random number if 0.5 is drawn. Either way, I hope this stops some of you feeling the frustration I felt when I got this error!

[Saul Teukolsky]

This is a good fix for this problem.

By the way, the 3rd edition has much faster routines for Poisson and other deviates which don't have this problem. (C++ only.)

Saul Teukolsky