1 // Copyright 2018 Developers of the Rand project.
2 //
3 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4 // https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
5 // <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
6 // option. This file may not be copied, modified, or distributed
7 // except according to those terms.
8 
9 use Rng;
10 use distributions::{Distribution, Uniform};
11 
12 /// Samples uniformly from the edge of the unit circle in two dimensions.
13 ///
14 /// Implemented via a method by von Neumann[^1].
15 ///
16 ///
17 /// # Example
18 ///
19 /// ```
20 /// use rand::distributions::{UnitCircle, Distribution};
21 ///
22 /// let circle = UnitCircle::new();
23 /// let v = circle.sample(&mut rand::thread_rng());
24 /// println!("{:?} is from the unit circle.", v)
25 /// ```
26 ///
27 /// [^1]: von Neumann, J. (1951) [*Various Techniques Used in Connection with
28 ///       Random Digits.*](https://mcnp.lanl.gov/pdf_files/nbs_vonneumann.pdf)
29 ///       NBS Appl. Math. Ser., No. 12. Washington, DC: U.S. Government Printing
30 ///       Office, pp. 36-38.
31 #[derive(Clone, Copy, Debug)]
32 pub struct UnitCircle;
33 
34 impl UnitCircle {
35     /// Construct a new `UnitCircle` distribution.
36     #[inline]
new() -> UnitCircle37     pub fn new() -> UnitCircle {
38         UnitCircle
39     }
40 }
41 
42 impl Distribution<[f64; 2]> for UnitCircle {
43     #[inline]
sample<R: Rng + ?Sized>(&self, rng: &mut R) -> [f64; 2]44     fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> [f64; 2] {
45         let uniform = Uniform::new(-1., 1.);
46         let mut x1;
47         let mut x2;
48         let mut sum;
49         loop {
50             x1 = uniform.sample(rng);
51             x2 = uniform.sample(rng);
52             sum = x1*x1 + x2*x2;
53             if sum < 1. {
54                 break;
55             }
56         }
57         let diff = x1*x1 - x2*x2;
58         [diff / sum, 2.*x1*x2 / sum]
59     }
60 }
61 
62 #[cfg(test)]
63 mod tests {
64     use distributions::Distribution;
65     use super::UnitCircle;
66 
67     /// Assert that two numbers are almost equal to each other.
68     ///
69     /// On panic, this macro will print the values of the expressions with their
70     /// debug representations.
71     macro_rules! assert_almost_eq {
72         ($a:expr, $b:expr, $prec:expr) => (
73             let diff = ($a - $b).abs();
74             if diff > $prec {
75                 panic!(format!(
76                     "assertion failed: `abs(left - right) = {:.1e} < {:e}`, \
77                      (left: `{}`, right: `{}`)",
78                     diff, $prec, $a, $b));
79             }
80         );
81     }
82 
83     #[test]
norm()84     fn norm() {
85         let mut rng = ::test::rng(1);
86         let dist = UnitCircle::new();
87         for _ in 0..1000 {
88             let x = dist.sample(&mut rng);
89             assert_almost_eq!(x[0]*x[0] + x[1]*x[1], 1., 1e-15);
90         }
91     }
92 
93     #[test]
value_stability()94     fn value_stability() {
95         let mut rng = ::test::rng(2);
96         let dist = UnitCircle::new();
97         assert_eq!(dist.sample(&mut rng), [-0.8032118336637037, 0.5956935036263119]);
98         assert_eq!(dist.sample(&mut rng), [-0.4742919588505423, -0.880367615130018]);
99         assert_eq!(dist.sample(&mut rng), [0.9297328981467168, 0.368234623716601]);
100     }
101 }
102