blob: a5ec0e009ac929e32db83790405106a9a65be10f [file] [log] [blame]
 // Copyright 2018-2019 Developers of the Rand project. // // Licensed under the Apache License, Version 2.0 or the MIT license // , at your // option. This file may not be copied, modified, or distributed // except according to those terms. use num_traits::Float; use crate::{uniform::SampleUniform, Distribution, Uniform}; use rand::Rng; /// Samples uniformly from the surface of the unit sphere in three dimensions. /// /// Implemented via a method by Marsaglia[^1]. /// /// /// # Example /// /// ``` /// use rand_distr::{UnitSphere, Distribution}; /// /// let v: [f64; 3] = UnitSphere.sample(&mut rand::thread_rng()); /// println!("{:?} is from the unit sphere surface.", v) /// ``` /// /// [^1]: Marsaglia, George (1972). [*Choosing a Point from the Surface of a /// Sphere.*](https://doi.org/10.1214/aoms/1177692644) /// Ann. Math. Statist. 43, no. 2, 645--646. #[derive(Clone, Copy, Debug)] pub struct UnitSphere; impl Distribution<[F; 3]> for UnitSphere { #[inline] fn sample(&self, rng: &mut R) -> [F; 3] { let uniform = Uniform::new(F::from(-1.).unwrap(), F::from(1.).unwrap()); loop { let (x1, x2) = (uniform.sample(rng), uniform.sample(rng)); let sum = x1 * x1 + x2 * x2; if sum >= F::from(1.).unwrap() { continue; } let factor = F::from(2.).unwrap() * (F::one() - sum).sqrt(); return [x1 * factor, x2 * factor, F::from(1.).unwrap() - F::from(2.).unwrap() * sum]; } } } #[cfg(test)] mod tests { use super::UnitSphere; use crate::Distribution; /// Assert that two numbers are almost equal to each other. /// /// On panic, this macro will print the values of the expressions with their /// debug representations. macro_rules! assert_almost_eq { (\$a:expr, \$b:expr, \$prec:expr) => { let diff = (\$a - \$b).abs(); if diff > \$prec { panic!( "assertion failed: `abs(left - right) = {:.1e} < {:e}`, \ (left: `{}`, right: `{}`)", diff, \$prec, \$a, \$b ); } }; } #[test] fn norm() { let mut rng = crate::test::rng(1); for _ in 0..1000 { let x: [f64; 3] = UnitSphere.sample(&mut rng); assert_almost_eq!(x[0] * x[0] + x[1] * x[1] + x[2] * x[2], 1., 1e-15); } } }