The R Project SVN R

Rev

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

Rev 59167 Rev 59170
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) 1998 Ross Ihaka
 3 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 2000-2001 The R Core Team
 4 
 *  Copyright (C) 2000-2012 The R Core Team
5 
 *  Copyright (C) 2005	The R Foundation
 5 
 *  Copyright (C) 2005	The R Foundation
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
9 
 *  the Free Software Foundation; either version 2 of the License, or
 9 
 *  the Free Software Foundation; either version 2 of the License, or
Line 146... Line 146...
146 
	} else {
 146 
	} else {
147 
	    k = kk;
 147 
	    k = kk;
148 
	}
 148 
	}
149 
    }
 149 
    }
150 
    if (setup1 || setup2) {
 150 
    if (setup1 || setup2) {
151 
	m = (k + 1.0) * (n1 + 1.0) / (tn + 2.0);
 151 
	m = (int) ((k + 1.0) * (n1 + 1.0) / (tn + 2.0));
152 
	minjx = imax2(0, k - n2);
 152 
	minjx = imax2(0, k - n2);
153 
	maxjx = imin2(n1, k);
 153 
	maxjx = imin2(n1, k);
154 
    }
 154 
    }
155 
    /* generate random variate --- Three basic cases */
 155 
    /* generate random variate --- Three basic cases */
156 
 
 156 
 
Line 223... Line 223...
223 
	}
 223 
	}
224 
      L30:
 224 
      L30:
225 
	u = unif_rand() * p3;
 225 
	u = unif_rand() * p3;
226 
	v = unif_rand();
 226 
	v = unif_rand();
227 
	if (u < p1) {		/* rectangular region */
 227 
	if (u < p1) {		/* rectangular region */
228 
	    ix = xl + u;
 228 
	    ix = (int) (xl + u);
229 
	} else if (u <= p2) {	/* left tail */
 229 
	} else if (u <= p2) {	/* left tail */
230 
	    ix = xl + log(v) / lamdl;
 230 
	    ix = (int) (xl + log(v) / lamdl);
231 
	    if (ix < minjx)
 231 
	    if (ix < minjx)
232 
		goto L30;
 232 
		goto L30;
233 
	    v = v * (u - p1) * lamdl;
 233 
	    v = v * (u - p1) * lamdl;
234 
	} else {		/* right tail */
 234 
	} else {		/* right tail */
235 
	    ix = xr - log(v) / lamdr;
 235 
	    ix = (int) (xr - log(v) / lamdr);
236 
	    if (ix > maxjx)
 236 
	    if (ix > maxjx)
237 
		goto L30;
 237 
		goto L30;
238 
	    v = v * (u - p2) * lamdr;
 238 
	    v = v * (u - p2) * lamdr;
239 
	}
 239 
	}
240 
 
 240