The R Project SVN – R – Diff – /trunk/src/nmath/rhyper.c

 /*
  *  Mathlib : A C Library of Special Functions
  *  Copyright (C) 1998 Ross Ihaka
  *  Copyright (C) 2000-2001 The R Development Core Team
+ *  Copyright (C) 2005	The R Foundation
+ *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation; either version 2 of the License, or
  *  (at your option) any later version.
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
+ *
- *  You should have received a copy of the GNU General Public License
+ *  A copy of the GNU General Public License is available via WWW at
+ *  http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
- *  along with this program; if not, write to the Free Software
+ *  writing to the Free Software Foundation, Inc., 59 Temple Place,
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
+ *  Suite 330, Boston, MA  02111-1307  USA.
+ *
  *  SYNOPSIS
+ *
  *    #include <Rmath.h>
  *    double rhyper(double NR, double NB, double n);
  *  REFERENCE
+ *
  *    V. Kachitvichyanukul and B. Schmeiser (1985).
  *    ``Computer generation of hypergeometric random variates,''
  *    Journal of Statistical Computation and Simulation 22, 127-145.
+ *
+ *    The original algorithm had a bug -- R bug report PR#7314 --
+ *    giving numbers slightly too small in case III h2pe
+ *    where (m < 100 || ix <= 50) , see below.
  */
 #include "nmath.h"
 /* afc(i) :=  ln( i! )	[logarithm of the factorial i.
     /* extern double afc(int); */
     int nn1, nn2, kk;
     int i, ix;
     Rboolean reject, setup1, setup2;
     double e, f, g, p, r, t, u, v, y;
     double de, dg, dr, ds, dt, gl, gu, nk, nm, ub;
     double xk, xm, xn, y1, ym, yn, yk, alv;
     /* These should become `thread_local globals' : */
 	/* acceptance/rejection test */
 	if (m < 100 || ix <= 50) {
 	    /* explicit evaluation */
-            /* The original algorithm (and TOMS 668) have
+	    /* The original algorithm (and TOMS 668) have
-                   f = f * i * (n2 - k + i) / (n1 - i) / (k - i);
+		   f = f * i * (n2 - k + i) / (n1 - i) / (k - i);
 	       in the (m > ix) case, but the definition of the
-               recurrence relation on p134 shows that the +1 is
+	       recurrence relation on p134 shows that the +1 is
-               needed. */
+	       needed. */
 	    f = 1.0;
 	    if (m < ix) {
 		for (i = m + 1; i <= ix; i++)
 		    f = f * (n1 - i + 1) * (k - i + 1) / (n2 - k + i) / i;
 	    } else if (m > ix) {
 		if (e < 0.0)
 		    de /= (1.0 + e);
 		if (alv < ub - 0.25 * (dr + ds + dt + de)
 		    + (y + m) * (gl - gu) - deltal) {
 		    reject = FALSE;
 		}
 		else {
 		    /* * Stirling's formula to machine accuracy
 		     */
 		    if (alv <= (a - afc(ix) - afc(n1 - ix)
 				- afc(k - ix) - afc(n2 - k + ix))) {

The R Project SVN R

(root)/trunk/src/nmath/rhyper.c – Rev 33687 → 33718