@@ -35,6 +35,7 @@ pub struct Poisson {
35
35
// precalculated values
36
36
exp_lambda : f64 ,
37
37
log_lambda : f64 ,
38
+ sqrt_2lambda : f64 ,
38
39
magic_val : f64 ,
39
40
}
40
41
@@ -43,11 +44,13 @@ impl Poisson {
43
44
/// `lambda`. Panics if `lambda <= 0`.
44
45
pub fn new ( lambda : f64 ) -> Poisson {
45
46
assert ! ( lambda > 0.0 , "Poisson::new called with lambda <= 0" ) ;
47
+ let log_lambda = lambda. ln ( ) ;
46
48
Poisson {
47
49
lambda : lambda,
48
50
exp_lambda : ( -lambda) . exp ( ) ,
49
- log_lambda : lambda. ln ( ) ,
50
- magic_val : lambda * lambda. ln ( ) - log_gamma ( 1.0 + lambda) ,
51
+ log_lambda : log_lambda,
52
+ sqrt_2lambda : ( 2.0 * lambda) . sqrt ( ) ,
53
+ magic_val : lambda * log_lambda - log_gamma ( 1.0 + lambda) ,
51
54
}
52
55
}
53
56
}
@@ -68,8 +71,6 @@ impl Distribution<u64> for Poisson {
68
71
}
69
72
// high expected values - rejection method
70
73
else {
71
- // some cached values
72
- let tmp = ( 2.0 * self . lambda ) . sqrt ( ) ;
73
74
let mut int_result: u64 ;
74
75
75
76
loop {
@@ -82,7 +83,7 @@ impl Distribution<u64> for Poisson {
82
83
// draw from the lorentzian distribution
83
84
comp_dev = ( PI * rng. gen :: < f64 > ( ) ) . tan ( ) ;
84
85
// shift the peak of the comparison ditribution
85
- result = tmp * comp_dev + self . lambda ;
86
+ result = self . sqrt_2lambda * comp_dev + self . lambda ;
86
87
// repeat the drawing until we are in the range of possible values
87
88
if result >= 0.0 {
88
89
break ;
0 commit comments