[Octave-cvsupdate] octave-forge/FIXES rand.cc,1.17,1.18 randpoisson.c,1.1,1.2

[Paul K. <pki...@us...>]

Update of /cvsroot/octave/octave-forge/FIXES
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv1534

Modified Files:
	rand.cc randpoisson.c 
Log Message:
Test extreme values; set randp(inf) to NaN --- is this correct?


Index: rand.cc
===================================================================
RCS file: /cvsroot/octave/octave-forge/FIXES/rand.cc,v
retrieving revision 1.17
retrieving revision 1.18
diff -u -d -r1.17 -r1.18
--- rand.cc	22 Feb 2005 01:53:42 -0000	1.17
+++ rand.cc	22 Feb 2005 02:27:17 -0000	1.18
@@ -454,6 +454,7 @@
 %!test
 %! % statistical tests may fail occasionally
 %! x = rande(100000,1);
+%! assert(min(x)>0); % *** Please report this!!! ***
 %! assert(mean(x),1,0.01);
 %! assert(var(x),1,0.03);
 %! assert(skewness(x),2,0.06);
@@ -548,10 +549,12 @@
 #define LGAMMA xlgamma
 
 /*
+%!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report
 %!test
 %! % statistical tests may fail occasionally.
 %! for a=[5 15]
 %!   x = randp(a,100000,1);
+%!   assert(min(x)>=0); % *** Please report this!!! ***
 %!   assert(mean(x),a,0.03);
 %!   assert(var(x),a,0.2);
 %!   assert(skewness(x),1/sqrt(a),0.03);
@@ -561,6 +564,7 @@
 %! % statistical tests may fail occasionally.
 %! for a=[5 15]
 %!   x = randp(a*ones(100000,1),100000,1);
+%!   assert(min(x)>=0); % *** Please report this!!! ***
 %!   assert(mean(x),a,0.03);
 %!   assert(var(x),a,0.2);
 %!   assert(skewness(x),1/sqrt(a),0.03);
@@ -644,6 +648,7 @@
 }
 
 /*
+%!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report
 %!test
 %! % statistical tests may fail occasionally.
 %! a=0.1; x = randg(a,100000,1);

Index: randpoisson.c
===================================================================
RCS file: /cvsroot/octave/octave-forge/FIXES/randpoisson.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -d -r1.1 -r1.2
--- randpoisson.c	7 Oct 2004 03:27:31 -0000	1.1
+++ randpoisson.c	22 Feb 2005 02:27:17 -0000	1.2
@@ -5,6 +5,7 @@
  * RUNI: uniform generator on (0,1)
  * RNOR: normal generator
  * LGAMMA: log gamma function
+ * INFINITE: function to test whether a value is infinite
  */
 
 
@@ -333,7 +334,7 @@
 void fill_randp(double L, size_t n, double *p)
 {
   size_t i;
-  if (L < 0.0) {
+  if (L < 0.0 || INFINITE(L)) {
     for (i=0; i<n; i++) p[i] = NAN;
   } else if (L <= 10.0) {
     poisson_cdf_lookup(L, p, n);
@@ -367,6 +368,10 @@
   } else if (L <= 1e8) {
     /* numerical recipes */
     poisson_rejection(L, &ret, 1);
+  } else if (INFINITE(L)) {
+    /* XXX FIXME XXX R uses NaN, but the normal approx. suggests that as
+     * limit should be inf. Which is correct? */
+    ret = NAN;
   } else {
     /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
     ret = floor(RNOR*sqrt(L) + L + 0.5);