The R Project SVN R

Rev

Rev 8431 | Rev 34956 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 8431 Rev 8628
Line 28... Line 28...
28 
double qhyper(double p, double NR, double NB, double n,
 28 
double qhyper(double p, double NR, double NB, double n,
29 
	      int lower_tail, int log_p)
 29 
	      int lower_tail, int log_p)
30 
{
 30 
{
31 
/* This is basically the same code as  ./phyper.c -- keep in sync! */
 31 
/* This is basically the same code as  ./phyper.c -- keep in sync! */
32 
    double N, xstart, xend, xr, xb, sum, term;
 32 
    double N, xstart, xend, xr, xb, sum, term;
- 
 
 33 
    int small_N;
- 
 
 34 
 
33 
#ifdef IEEE_754
 35 
#ifdef IEEE_754
34 
    if (ISNAN(p) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
 36 
    if (ISNAN(p) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
35 
	return p + NR + NB + n;
 37 
	return p + NR + NB + n;
36 
#endif
 38 
#endif
37 
    if(!R_FINITE(p) || !R_FINITE(NR) || !R_FINITE(NB) || !R_FINITE(n))
 39 
    if(!R_FINITE(p) || !R_FINITE(NR) || !R_FINITE(NB) || !R_FINITE(n))
Line 51... Line 53...
51 
    if(p == R_DT_1) return xend;
 53 
    if(p == R_DT_1) return xend;
52 
 
 54 
 
53 
    xr = xstart;
 55 
    xr = xstart;
54 
    xb = n - xr;
 56 
    xb = n - xr;
55 
 
 57 
 
- 
 
 58 
    small_N = (N < 1000); /* won't have underflow in product below */
- 
 
 59 
    /* if N is small,  term := product.ratio( bin.coef );
- 
 
 60 
       otherwise work with its logarithm to protect against underflow */
56 
    term = exp(lfastchoose(NR, xr) + lfastchoose(NB, xb) - lfastchoose(N, n));
 61 
    term = lfastchoose(NR, xr) + lfastchoose(NB, xb) - lfastchoose(N, n);
- 
 
 62 
    if(small_N) term = exp(term);
57 
    NR = NR - xr;
 63 
    NR -= xr;
58 
    NB = NB - xb;
 64 
    NB -= xb;
- 
 
 65 
 
59 
    sum = term;
 66 
    sum = term;
60 
    if(!lower_tail || log_p) {
 67 
    if(!lower_tail || log_p) {
61 
	p = R_DT_qIv(p);
 68 
	p = R_DT_qIv(p);
62 
    }
 69 
    }
63 
    p *= 1 - 64*DBL_EPSILON;
 70 
    p *= 1 - 64*DBL_EPSILON;
- 
 
 71 
    sum = small_N ? term : exp(term);
- 
 
 72 
 
64 
    while(sum < p && xr < xend) {
 73 
    while(sum < p && xr < xend) {
65 
	xr++;
 74 
	xr++;
66 
	NB++;
 75 
	NB++;
- 
 
 76 
	if (small_N) term *= (NR / xr) * (xb / NB);
67 
	term *= (NR / xr) * (xb / NB);
 77 
	else term += log((NR / xr) * (xb / NB));
68 
	sum += term;
 78 
	sum += small_N ? term : exp(term);
69 
	xb--;
 79 
	xb--;
70 
	NR--;
 80 
	NR--;
71 
    }
 81 
    }
72 
    return xr;
 82 
    return xr;
73 
}
 83 
}