/home/docs/checkouts/readthedocs.org/user_builds/advanced-micro-devices-rocrand/checkouts/develop/projects/rocrand/library/include/rocrand/rocrand_poisson.h Source File

• library
• include
• rocrand

rocrand_poisson.h

1 // Copyright (c) 2017-2025 Advanced Micro Devices, Inc. All rights reserved.

2 //

3 // Permission is hereby granted, free of charge, to any person obtaining a copy

4 // of this software and associated documentation files (the "Software"), to deal

5 // in the Software without restriction, including without limitation the rights

6 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

7 // copies of the Software, and to permit persons to whom the Software is

8 // furnished to do so, subject to the following conditions:

9 //

10 // The above copyright notice and this permission notice shall be included in

11 // all copies or substantial portions of the Software.

12 //

13 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

14 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

15 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

16 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

17 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

18 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN

19 // THE SOFTWARE.

20 

21 #ifndef ROCRAND_POISSON_H_

22 #define ROCRAND_POISSON_H_

23 

29 #include <math.h>

30 

31 #include "rocrand/rocrand_lfsr113.h"

32 #include "rocrand/rocrand_mrg31k3p.h"

33 #include "rocrand/rocrand_mrg32k3a.h"

34 #include "rocrand/rocrand_mtgp32.h"

35 #include "rocrand/rocrand_philox4x32_10.h"

36 #include "rocrand/rocrand_scrambled_sobol32.h"

37 #include "rocrand/rocrand_scrambled_sobol64.h"

38 #include "rocrand/rocrand_sobol32.h"

39 #include "rocrand/rocrand_sobol64.h"

40 #include "rocrand/rocrand_threefry2x32_20.h"

41 #include "rocrand/rocrand_threefry2x64_20.h"

42 #include "rocrand/rocrand_threefry4x32_20.h"

43 #include "rocrand/rocrand_threefry4x64_20.h"

44 #include "rocrand/rocrand_xorwow.h"

45 

46 #include "rocrand/rocrand_normal.h"

47 #include "rocrand/rocrand_uniform.h"

48 

49 #include <hip/hip_runtime.h>

50 

51 namespace rocrand_device {

52 namespace detail {

53 

54 constexpr double lambda_threshold_small = 64.0;

55 constexpr double lambda_threshold_huge = 4000.0;

56 

57 template<class State, typename Result_Type = unsigned int>

58 __forceinline__ __device__ __host__ Result_Type poisson_distribution_small(State& state,

59  double lambda)

60 {

61  // Knuth's method

62 

63  const double limit = exp(-lambda);

64  Result_Type k = 0;

65  double product = 1.0;

66 

67  do

68  {

69  k++;

70  product *= rocrand_uniform_double(state);

71  }

72  while (product > limit);

73 

74  return k - 1;

75 }

76 

77 __forceinline__ __device__ __host__ double lgamma_approx(const double x)

78 {

79  // Lanczos approximation (g = 7, n = 9)

80 

81  const double z = x - 1.0;

82 

83  const int g = 7;

84  const int n = 9;

85  const double coefs[n] = {

86  0.99999999999980993227684700473478,

87  676.520368121885098567009190444019,

88  -1259.13921672240287047156078755283,

89  771.3234287776530788486528258894,

90  -176.61502916214059906584551354,

91  12.507343278686904814458936853,

92  -0.13857109526572011689554707,

93  9.984369578019570859563e-6,

94  1.50563273514931155834e-7

95  };

96  double sum = 0.0;

97  #pragma unroll

98  for (int i = n - 1; i > 0; i--)

99  {

100  sum += coefs[i] / (z + i);

101  }

102  sum += coefs[0];

103 

104  const double log_sqrt_2_pi = 0.9189385332046727418;

105  const double e = 2.718281828459045090796;

106  return (log_sqrt_2_pi + log(sum) - g) + (z + 0.5) * log((z + g + 0.5) / e);

107 }

108 

109 template<class State, typename Result_Type = unsigned int>

110 __forceinline__ __device__ __host__ Result_Type poisson_distribution_large(State& state,

111  double lambda)

112 {

113  // Rejection method PA, A. C. Atkinson

114 

115  const double c = 0.767 - 3.36 / lambda;

116  const double beta = ROCRAND_PI_DOUBLE / sqrt(3.0 * lambda);

117  const double alpha = beta * lambda;

118  const double k = log(c) - lambda - log(beta);

119  const double log_lambda = log(lambda);

120 

121  while (true)

122  {

123  const double u = rocrand_uniform_double(state);

124  const double x = (alpha - log((1.0 - u) / u)) / beta;

125  const double n = floor(x + 0.5);

126  if (n < 0)

127  {

128  continue;

129  }

130  const double v = rocrand_uniform_double(state);

131  const double y = alpha - beta * x;

132  const double t = 1.0 + exp(y);

133  const double lhs = y + log(v / (t * t));

134  const double rhs = k + n * log_lambda - lgamma_approx(n + 1.0);

135  if (lhs <= rhs)

136  {

137  return static_cast<Result_Type>(n);

138  }

139  }

140 }

141 

142 template<class State, typename Result_Type = unsigned int>

143 __forceinline__ __device__ __host__ Result_Type poisson_distribution_huge(State& state,

144  double lambda)

145 {

146  // Approximate Poisson distribution with normal distribution

147 

148  const double n = rocrand_normal_double(state);

149  return static_cast<Result_Type>(round(sqrt(lambda) * n + lambda));

150 }

151 

152 template<class State, typename Result_Type = unsigned int>

153 __forceinline__ __device__ __host__ Result_Type poisson_distribution(State& state, double lambda)

154 {

155  if (lambda < lambda_threshold_small)

156  {

157  return poisson_distribution_small<State, Result_Type>(state, lambda);

158  }

159  else if (lambda <= lambda_threshold_huge)

160  {

161  return poisson_distribution_large<State, Result_Type>(state, lambda);

162  }

163  else

164  {

165  return poisson_distribution_huge<State, Result_Type>(state, lambda);

166  }

167 }

168 

169 template<class State, typename Result_Type = unsigned int>

170 __forceinline__ __device__ __host__ Result_Type poisson_distribution_itr(State& state,

171  double lambda)

172 {

173  // Algorithm ITR

174  // George S. Fishman

175  // Discrete-Event Simulation: Modeling, Programming, and Analysis

176  // p. 333

177  double L;

178  double x = 1.0;

179  double y = 1.0;

180  Result_Type k = 0;

181  int pow = 0;

182  // Algorithm ITR uses u from (0, 1) and uniform_double returns (0, 1]

183  // Change u to ensure that 1 is never generated,

184  // otherwise the inner loop never ends.

185  double u = rocrand_uniform_double(state) - ROCRAND_2POW32_INV_DOUBLE / 2.0;

186  double upow = pow + 500.0;

187  double ex = exp(-500.0);

188  do{

189  if (lambda > upow)

190  L = ex;

191  else

192  L = exp((double)(pow - lambda));

193 

194  x *= L;

195  y *= L;

196  pow += 500;

197  while (u > y)

198  {

199  k++;

200  x *= ((double)lambda / (double) k);

201  y += x;

202  }

203  } while((double)pow < lambda);

204  return k;

205 }

206 

207 template<class State, typename Result_Type = unsigned int>

208 __forceinline__ __device__ __host__ Result_Type poisson_distribution_inv(State& state,

209  double lambda)

210 {

211  if (lambda < 1000.0)

212  {

213  return poisson_distribution_itr<State, Result_Type>(state, lambda);

214  }

215  else

216  {

217  return poisson_distribution_huge<State, Result_Type>(state, lambda);

218  }

219 }

220 

221 } // end namespace detail

222 } // end namespace rocrand_device

223 

235 #ifndef ROCRAND_DETAIL_BM_NOT_IN_STATE

236 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_philox4x32_10* state,

237  double lambda)

238 {

239  return rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(

240  state,

241  lambda);

242 }

243 

255 __forceinline__ __device__ __host__

256 uint4 rocrand_poisson4(rocrand_state_philox4x32_10* state, double lambda)

257 {

258  return uint4{

259  rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(

260  state,

261  lambda),

262  rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(

263  state,

264  lambda),

265  rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(

266  state,

267  lambda),

268  rocrand_device::detail::poisson_distribution<rocrand_state_philox4x32_10*, unsigned int>(

269  state,

270  lambda)};

271 }

272 #endif // ROCRAND_DETAIL_BM_NOT_IN_STATE

273 

285 #ifndef ROCRAND_DETAIL_BM_NOT_IN_STATE

286 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_mrg31k3p* state,

287  double lambda)

288 {

289  return rocrand_device::detail::poisson_distribution<rocrand_state_mrg31k3p*, unsigned int>(

290  state,

291  lambda);

292 }

293 #endif // ROCRAND_DETAIL_BM_NOT_IN_STATE

294 

306 #ifndef ROCRAND_DETAIL_BM_NOT_IN_STATE

307 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_mrg32k3a* state,

308  double lambda)

309 {

310  return rocrand_device::detail::poisson_distribution<rocrand_state_mrg32k3a*, unsigned int>(

311  state,

312  lambda);

313 }

314 #endif // ROCRAND_DETAIL_BM_NOT_IN_STATE

315 

327 #ifndef ROCRAND_DETAIL_BM_NOT_IN_STATE

328 __forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_xorwow* state,

329  double lambda)

330 {

331  return rocrand_device::detail::poisson_distribution<rocrand_state_xorwow*, unsigned int>(

332  state,

333  lambda);

334 }

335 #endif // ROCRAND_DETAIL_BM_NOT_IN_STATE

336 

348 __forceinline__ __device__ __host__

349 unsigned int rocrand_poisson(rocrand_state_mtgp32* state, double lambda)

350 {

351  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_mtgp32*, unsigned int>(

352  state,

353  lambda);

354 }

355 

367 __forceinline__ __device__ __host__

368 unsigned int rocrand_poisson(rocrand_state_sobol32* state, double lambda)

369 {

370  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_sobol32*, unsigned int>(

371  state,

372  lambda);

373 }

374 

386 __forceinline__ __device__ __host__

387 unsigned int rocrand_poisson(rocrand_state_scrambled_sobol32* state, double lambda)

388 {

389  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_scrambled_sobol32*,

390  unsigned int>(state, lambda);

391 }

392 

404 __forceinline__ __device__ __host__

405 unsigned int rocrand_poisson(rocrand_state_sobol64* state, double lambda)

406 {

407  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_sobol64*, unsigned int>(

408  state,

409  lambda);

410 }

411 

423 __forceinline__ __device__ __host__

424 unsigned int rocrand_poisson(rocrand_state_scrambled_sobol64* state, double lambda)

425 {

426  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_scrambled_sobol64*,

427  unsigned int>(state, lambda);

428 }

429 

441 __forceinline__ __device__ __host__

442 unsigned int rocrand_poisson(rocrand_state_lfsr113* state, double lambda)

443 {

444  return rocrand_device::detail::poisson_distribution_inv<rocrand_state_lfsr113*, unsigned int>(

445  state,

446  lambda);

447 }

448 

460 __forceinline__ __device__ __host__

461 unsigned int rocrand_poisson(rocrand_state_threefry2x32_20* state, double lambda)

462 {

463  return rocrand_device::detail::poisson_distribution_inv(state, lambda);

464 }

465 

477 __forceinline__ __device__ __host__

478 unsigned int rocrand_poisson(rocrand_state_threefry2x64_20* state, double lambda)

479 {

480  return rocrand_device::detail::poisson_distribution_inv(state, lambda);

481 }

482 

494 __forceinline__ __device__ __host__

495 unsigned int rocrand_poisson(rocrand_state_threefry4x32_20* state, double lambda)

496 {

497  return rocrand_device::detail::poisson_distribution_inv(state, lambda);

498 }

499 

511 __forceinline__ __device__ __host__

512 unsigned int rocrand_poisson(rocrand_state_threefry4x64_20* state, double lambda)

513 {

514  return rocrand_device::detail::poisson_distribution_inv(state, lambda);

515 }

516  // end of group rocranddevice

518 

519 #endif // ROCRAND_POISSON_H_

rocrand_poisson

__forceinline__ __device__ __host__ unsigned int rocrand_poisson(rocrand_state_philox4x32_10 *state, double lambda)

Returns a Poisson-distributed unsigned int using Philox generator.

Definition: rocrand_poisson.h:236

rocrand_poisson4

__forceinline__ __device__ __host__ uint4 rocrand_poisson4(rocrand_state_philox4x32_10 *state, double lambda)

Returns four Poisson-distributed unsigned int values using Philox generator.

Definition: rocrand_poisson.h:256

rocrand_uniform_double

__forceinline__ __device__ __host__ double rocrand_uniform_double(rocrand_state_philox4x32_10 *state)

Returns a uniformly distributed random double value from (0; 1] range.

Definition: rocrand_uniform.h:299

rocrand_normal_double

__forceinline__ __device__ __host__ double rocrand_normal_double(rocrand_state_philox4x32_10 *state)

Returns a normally distributed double value.

Definition: rocrand_normal.h:442