Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > Methods: Chapters 7, 8, and 20

Reply
 
Thread Tools Display Modes
  #1  
Old 04-11-2012, 07:11 AM
blackmime123 blackmime123 is offline
Registered User
 
Join Date: Apr 2012
Posts: 1
Possible bug in POIDEV

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!
Reply With Quote
  #2  
Old 04-11-2012, 09:22 AM
Saul Teukolsky Saul Teukolsky is offline
Numerical Recipes Author
 
Join Date: Dec 2001
Posts: 211
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
Reply With Quote
Reply

Tags
assertion fail, gammln_s arg, nrerror, poidev

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 11:16 PM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.