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