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 surface of the unit sphere in three dimensions.
13 ///
14 /// Implemented via a method by Marsaglia[^1].
15 ///
16 ///
17 /// # Example
18 ///
19 /// ```
20 /// use rand::distributions::{UnitSphereSurface, Distribution};
21 ///
22 /// let sphere = UnitSphereSurface::new();
23 /// let v = sphere.sample(&mut rand::thread_rng());
24 /// println!("{:?} is from the unit sphere surface.", v)
25 /// ```
26 ///
27 /// [^1]: Marsaglia, George (1972). [*Choosing a Point from the Surface of a
28 ///       Sphere.*](https://doi.org/10.1214/aoms/1177692644)
29 ///       Ann. Math. Statist. 43, no. 2, 645--646.
30 #[derive(Clone, Copy, Debug)]
31 pub struct UnitSphereSurface;
32 
33 impl UnitSphereSurface {
34     /// Construct a new `UnitSphereSurface` distribution.
35     #[inline]
new() -> UnitSphereSurface36     pub fn new() -> UnitSphereSurface {
37         UnitSphereSurface
38     }
39 }
40 
41 impl Distribution<[f64; 3]> for UnitSphereSurface {
42     #[inline]
sample<R: Rng + ?Sized>(&self, rng: &mut R) -> [f64; 3]43     fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> [f64; 3] {
44         let uniform = Uniform::new(-1., 1.);
45         loop {
46             let (x1, x2) = (uniform.sample(rng), uniform.sample(rng));
47             let sum = x1*x1 + x2*x2;
48             if sum >= 1. {
49                 continue;
50             }
51             let factor = 2. * (1.0_f64 - sum).sqrt();
52             return [x1 * factor, x2 * factor, 1. - 2.*sum];
53         }
54     }
55 }
56 
57 #[cfg(test)]
58 mod tests {
59     use distributions::Distribution;
60     use super::UnitSphereSurface;
61 
62     /// Assert that two numbers are almost equal to each other.
63     ///
64     /// On panic, this macro will print the values of the expressions with their
65     /// debug representations.
66     macro_rules! assert_almost_eq {
67         ($a:expr, $b:expr, $prec:expr) => (
68             let diff = ($a - $b).abs();
69             if diff > $prec {
70                 panic!(format!(
71                     "assertion failed: `abs(left - right) = {:.1e} < {:e}`, \
72                      (left: `{}`, right: `{}`)",
73                     diff, $prec, $a, $b));
74             }
75         );
76     }
77 
78     #[test]
norm()79     fn norm() {
80         let mut rng = ::test::rng(1);
81         let dist = UnitSphereSurface::new();
82         for _ in 0..1000 {
83             let x = dist.sample(&mut rng);
84             assert_almost_eq!(x[0]*x[0] + x[1]*x[1] + x[2]*x[2], 1., 1e-15);
85         }
86     }
87 
88     #[test]
value_stability()89     fn value_stability() {
90         let mut rng = ::test::rng(2);
91         let dist = UnitSphereSurface::new();
92         assert_eq!(dist.sample(&mut rng),
93                    [-0.24950027180862533, -0.7552572587896719, 0.6060825747478084]);
94         assert_eq!(dist.sample(&mut rng),
95                    [0.47604534507233487, -0.797200864987207, -0.3712837328763685]);
96         assert_eq!(dist.sample(&mut rng),
97                    [0.9795722330927367, 0.18692349236651176, 0.07414747571708524]);
98     }
99 }
100