diff options
Diffstat (limited to 'boost/compute/random/normal_distribution.hpp')
| -rw-r--r-- | boost/compute/random/normal_distribution.hpp | 24 |
1 files changed, 20 insertions, 4 deletions
diff --git a/boost/compute/random/normal_distribution.hpp b/boost/compute/random/normal_distribution.hpp index d025faeb2e..4693e4fffe 100644 --- a/boost/compute/random/normal_distribution.hpp +++ b/boost/compute/random/normal_distribution.hpp @@ -13,6 +13,9 @@ #include <limits> +#include <boost/assert.hpp> +#include <boost/type_traits.hpp> + #include <boost/compute/command_queue.hpp> #include <boost/compute/function.hpp> #include <boost/compute/types/fundamental.hpp> @@ -90,11 +93,19 @@ public: BOOST_COMPUTE_FUNCTION(RealType2, box_muller, (const uint2_ x), { - const RealType x1 = x.x / (RealType) (UINT_MAX - 1); - const RealType x2 = x.y / (RealType) (UINT_MAX - 1); + const RealType one = 1; + const RealType two = 2; + + // Use nextafter to push values down into [0,1) range; without this, floating point rounding can + // lead to have x1 = 1, but that would lead to taking the log of 0, which would result in negative + // infinities; by pushing the values off 1 towards 0, we ensure this won't happen. + const RealType x1 = nextafter(x.x / (RealType) UINT_MAX, (RealType) 0); + const RealType x2 = x.y / (RealType) UINT_MAX; - const RealType z1 = sqrt(-2.f * log2(x1)) * cos(2.f * M_PI_F * x2); - const RealType z2 = sqrt(-2.f * log2(x1)) * sin(2.f * M_PI_F * x2); + const RealType rho = sqrt(-two * log(one-x1)); + + const RealType z1 = rho * cos(two * M_PI_F * x2); + const RealType z2 = rho * sin(two * M_PI_F * x2); return (RealType2)(MEAN, MEAN) + (RealType2)(z1, z2) * (RealType2)(STDDEV, STDDEV); }); @@ -116,6 +127,11 @@ public: private: RealType m_mean; RealType m_stddev; + + BOOST_STATIC_ASSERT_MSG( + boost::is_floating_point<RealType>::value, + "Template argument must be a floating point type" + ); }; } // end compute namespace |