← Back to branch summary

~ubuntu-branches/ubuntu/utopic/libthrust/utopic

« back to all changes in this revision

Viewing changes to random/detail/normal_distribution_base.h

  • Committer: Package Import Robot
  • Author(s): Andreas Beckmann
  • Date: 2013-07-10 12:57:33 UTC
  • mfrom: (1.1.4)
  • Revision ID: package-import@ubuntu.com-20130710125733-my19jic71sqsabaj
Tags: 1.7.0-1
http://bugs.debian.org/715362
* New upstream release.  (Closes: #715362)
* Update watch file.

Show diffs side-by-side

added added

removed removed

Lines of Context:
random/detail/normal_distribution_base.h
 
1
/*
 
2
 *  Copyright 2008-2012 NVIDIA Corporation
 
3
 *
 
4
 *  Licensed under the Apache License, Version 2.0 (the "License");
 
5
 *  you may not use this file except in compliance with the License.
 
6
 *  You may obtain a copy of the License at
 
7
 *
 
8
 *      http://www.apache.org/licenses/LICENSE-2.0
 
9
 *
 
10
 *  Unless required by applicable law or agreed to in writing, software
 
11
 *  distributed under the License is distributed on an "AS IS" BASIS,
 
12
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 
13
 *  See the License for the specific language governing permissions and
 
14
 *  limitations under the License.
 
15
 */
 
16
 
 
17
/*
 
18
 * Copyright Jens Maurer 2000-2001
 
19
 * Distributed under the Boost Software License, Version 1.0. (See
 
20
 * accompanying file LICENSE_1_0.txt or copy at
 
21
 * http://www.boost.org/LICENSE_1_0.txt)
 
22
 */
 
23
 
 
24
#pragma once
 
25
 
 
26
#include <thrust/detail/config.h>
 
27
#include <thrust/pair.h>
 
28
#include <thrust/random/uniform_real_distribution.h>
 
29
#include <limits>
 
30
#include <cmath>
 
31
 
 
32
namespace thrust
 
33
{
 
34
namespace random
 
35
{
 
36
namespace detail
 
37
{
 
38
 
 
39
// this version samples the normal distribution directly 
 
40
// and uses the non-standard math function erfcinv
 
41
template<typename RealType>
 
42
  class normal_distribution_nvcc
 
43
{
 
44
  protected:
 
45
    template<typename UniformRandomNumberGenerator>
 
46
    __host__ __device__
 
47
    RealType sample(UniformRandomNumberGenerator &urng, const RealType mean, const RealType stddev)
 
48
    {
 
49
      typedef typename UniformRandomNumberGenerator::result_type uint_type;
 
50
      const uint_type urng_range = UniformRandomNumberGenerator::max - UniformRandomNumberGenerator::min;
 
51
 
 
52
      // Constants for conversion
 
53
      const RealType S1 = static_cast<RealType>(1) / urng_range;
 
54
      const RealType S2 = S1 / 2;
 
55
 
 
56
      RealType S3 = static_cast<RealType>(-1.4142135623730950488016887242097); // -sqrt(2)
 
57
      
 
58
      // Get the integer value
 
59
      uint_type u = urng() - UniformRandomNumberGenerator::min;
 
60
 
 
61
      // Ensure the conversion to float will give a value in the range [0,0.5)
 
62
      if(u > (urng_range / 2))
 
63
      {
 
64
        u = urng_range - u;
 
65
        S3 = -S3;
 
66
      }
 
67
 
 
68
      // Convert to floating point in [0,0.5)
 
69
      RealType p = u*S1 + S2;
 
70
 
 
71
      // Apply inverse error function
 
72
      return mean + stddev * S3 * erfcinv(2 * p);
 
73
    }
 
74
 
 
75
    // no-op
 
76
    __host__ __device__
 
77
    void reset() {}
 
78
};
 
79
 
 
80
// this version samples the normal distribution using 
 
81
// Marsaglia's "polar method"
 
82
template<typename RealType>
 
83
  class normal_distribution_portable
 
84
{
 
85
  protected:
 
86
    normal_distribution_portable()
 
87
      : m_valid(false)
 
88
    {}
 
89
 
 
90
    normal_distribution_portable(const normal_distribution_portable &other)
 
91
      : m_valid(other.m_valid)
 
92
    {}
 
93
 
 
94
    void reset()
 
95
    {
 
96
      m_valid = false;
 
97
    }
 
98
 
 
99
    // note that we promise to call this member function with the same mean and stddev
 
100
    template<typename UniformRandomNumberGenerator>
 
101
    __host__ __device__
 
102
    RealType sample(UniformRandomNumberGenerator &urng, const RealType mean, const RealType stddev)
 
103
    {
 
104
      // implementation from Boost
 
105
      // allow for Koenig lookup
 
106
      using std::sqrt; using std::log; using std::sin; using std::cos;
 
107
 
 
108
      if(!m_valid)
 
109
      {
 
110
        uniform_real_distribution<RealType> u01;
 
111
        m_r1 = u01(urng);
 
112
        m_r2 = u01(urng);
 
113
        m_cached_rho = sqrt(-RealType(2) * log(RealType(1)-m_r2));
 
114
 
 
115
        m_valid = true;
 
116
      }
 
117
      else
 
118
      {
 
119
        m_valid = false;
 
120
      }
 
121
 
 
122
      const RealType pi = RealType(3.14159265358979323846);
 
123
 
 
124
      RealType result = m_cached_rho * (m_valid ?
 
125
                          cos(RealType(2)*pi*m_r1) :
 
126
                          sin(RealType(2)*pi*m_r1));
 
127
 
 
128
      return result;
 
129
    }
 
130
 
 
131
  private:
 
132
    RealType m_r1, m_r2, m_cached_rho;
 
133
    bool m_valid;
 
134
};
 
135
 
 
136
template<typename RealType>
 
137
  struct normal_distribution_base
 
138
{
 
139
#if THRUST_DEVICE_COMPILER == THRUST_DEVICE_COMPILER_NVCC
 
140
  typedef normal_distribution_nvcc<RealType> type;
 
141
#else
 
142
  typedef normal_distribution_portable<RealType> type;
 
143
#endif
 
144
};
 
145
 
 
146
} // end detail
 
147
} // end random
 
148
} // end thrust
 
149