The R Project SVN R

Rev

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

Rev 77685 Rev 88401
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--2024 The R Core Team
3 
 *  Copyright (C) 1998 Ross Ihaka
 4 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 2000--2016 The R Core Team
 - 
 
5 
 *
 5 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 6 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
 7 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
 8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 9 
 *  (at your option) any later version.
Line 100... Line 100...
100 
	    p = e * unif_rand();
 100 
	    p = e * unif_rand();
101 
	    if (p >= 1.0) {
 101 
	    if (p >= 1.0) {
102 
		x = -log((e - p) / a);
 102 
		x = -log((e - p) / a);
103 
		if (exp_rand() >= (1.0 - a) * log(x))
 103 
		if (exp_rand() >= (1.0 - a) * log(x))
104 
		    break;
 104 
		    break;
105 
	    } else {
 105 
	    } else { /* p < 1  <==>  log(p) < 0 */
106 
		x = exp(log(p) / a);
 106 
		x = exp(log(p) / a); // exp(*) "often" underflows to 0 for small a (e.g. a = 1/1000)
107 
		if (exp_rand() >= x)
 107 
		if (exp_rand() >= x)
108 
		    break;
 108 
		    break;
109 
	    }
 109 
	    }
110 
	}
 110 
	}
111 
	return scale * x;
 111 
	return scale * x;