The R Project SVN R

Rev

Rev 87974 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 87974 Rev 90223
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 2000-2025 The R Core Team
 3 
 *  Copyright (C) 2000-2026 The R Core Team
4 
 *  Copyright (C) 2005-2020 The R Foundation
 4 
 *  Copyright (C) 2005-2020 The R Foundation
5 
 *  Copyright (C) 1998 Ross Ihaka
 5 
 *  Copyright (C) 1998 Ross Ihaka
6 
 *
 6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
 7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
 8 
 *  it under the terms of the GNU General Public License as published by
Line 80... Line 80...
80 
}
 80 
}
81 
 
 81 
 
82 
//     rhyper(NR, NB, n) -- NR 'red', NB 'blue', n drawn, how many are 'red'
 82 
//     rhyper(NR, NB, n) -- NR 'red', NB 'blue', n drawn, how many are 'red'
83 
double rhyper(double nn1in, double nn2in, double kkin)
 83 
double rhyper(double nn1in, double nn2in, double kkin)
84 
{
 84 
{
85 
    /* extern double afc(int); */
 - 
 
86 
 
 - 
 
87 
    /* check parameter validity */
 85 
    /* check parameter validity */
88 
 
 86 
 
89 
    if(!R_FINITE(nn1in) || !R_FINITE(nn2in) || !R_FINITE(kkin))
 87 
    if(!R_FINITE(nn1in) || !R_FINITE(nn2in) || !R_FINITE(kkin))
90 
	ML_WARN_return_NAN;
 88 
	ML_WARN_return_NAN;
91 
 
 89 
 
92 
    nn1in = R_forceint(nn1in);
 90 
    nn1in = R_forceint(nn1in);
93 
    nn2in = R_forceint(nn2in);
 91 
    nn2in = R_forceint(nn2in);
94 
    kkin  = R_forceint(kkin);
 92 
    kkin  = R_forceint(kkin);
- 
 
 93 
    double N = nn1in + nn2in;
95 
 
 94 
 
96 
    if (nn1in < 0 || nn2in < 0 || kkin < 0 || kkin > nn1in + nn2in)
 95 
    if (nn1in < 0 || nn2in < 0 || kkin < 0 || kkin > nn1in + nn2in)
97 
	ML_WARN_return_NAN;
 96 
	ML_WARN_return_NAN;
98 
    if (nn1in >= INT_MAX || nn2in >= INT_MAX || kkin >= INT_MAX) {
 97 
    if (N > INT_MAX || kkin >= INT_MAX) {
99 
	/* large n -- evade integer overflow (and inappropriate algorithms)
 98 
	/* large n -- evade integer overflow (and inappropriate algorithms)
100 
	   -------- */
 99 
	   -------- */
101 
#ifdef DEBUG_rhyper
 100 
#ifdef DEBUG_rhyper
102 
	REprintf("rhyper(nn1=%.0f, nn2=%.0f, kk=%.0f): 'large n' case\n", nn1in, nn2in, kkin);
 101 
	REprintf("rhyper(nn1=%.0f, nn2=%.0f, kk=%.0f): 'large n' case\n", nn1in, nn2in, kkin);
103 
#endif
 102 
#endif
Line 119... Line 118...
119 
 
 118 
 
120 
    /* These should become 'thread_local globals' : */
 119 
    /* These should become 'thread_local globals' : */
121 
    static int ks = -1, n1s = -1, n2s = -1;
 120 
    static int ks = -1, n1s = -1, n2s = -1;
122 
    static int m, minjx, maxjx;
 121 
    static int m, minjx, maxjx;
123 
    static int k, n1, n2; // <- not allowing larger integer par
 122 
    static int k, n1, n2; // <- not allowing larger integer par
124 
    static double N;
 - 
 
125 
 
 123 
 
126 
    bool setup1, setup2;
 124 
    bool setup1, setup2;
127 
    /* if new parameter values, initialize */
 125 
    /* if new parameter values, initialize */
128 
    if (nn1 != n1s || nn2 != n2s) { // n1 | n2 is changed: setup all
 126 
    if (nn1 != n1s || nn2 != n2s) { // n1 | n2 is changed: setup all
129 
	setup1 = true;	setup2 = true;
 127 
	setup1 = true;	setup2 = true;
Line 132... Line 130...
132 
    } else { // all three unchanged ==> no setup
 130 
    } else { // all three unchanged ==> no setup
133 
	setup1 = false;	setup2 = false;
 131 
	setup1 = false;	setup2 = false;
134 
    }
 132 
    }
135 
    if (setup1) { // n1 & n2
 133 
    if (setup1) { // n1 & n2
136 
	n1s = nn1; n2s = nn2; // save
 134 
	n1s = nn1; n2s = nn2; // save
137 
	N = nn1 + (double)nn2; // avoid int overflow
 - 
 
138 
	if (nn1 <= nn2) {
 135 
	if (nn1 <= nn2) {
139 
	    n1 = nn1; n2 = nn2;
 136 
	    n1 = nn1; n2 = nn2;
140 
	} else { // nn2 < nn1
 137 
	} else { // nn2 < nn1
141 
	    n1 = nn2; n2 = nn1;
 138 
	    n1 = nn2; n2 = nn1;
142 
	}
 139 
	}