vSMC
vSMC: Scalable Monte Carlo
gammak1.h
Go to the documentation of this file.
1 //============================================================================
2 // vSMC/include/vsmc/rng/gammak1.h
3 //----------------------------------------------------------------------------
4 // vSMC: Scalable Monte Carlo
5 //----------------------------------------------------------------------------
6 // Copyright (c) 2013,2014, Yan Zhou
7 // All rights reserved.
8 //
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are met:
11 //
12 // Redistributions of source code must retain the above copyright notice,
13 // this list of conditions and the following disclaimer.
14 //
15 // Redistributions in binary form must reproduce the above copyright notice,
16 // this list of conditions and the following disclaimer in the documentation
17 // and/or other materials provided with the distribution.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 // POSSIBILITY OF SUCH DAMAGE.
30 //============================================================================
31 
32 #ifndef VSMC_RNG_GAMMAK1_H
33 #define VSMC_RNG_GAMMAK1_H
34 
81 
82 #include <vsmc/rng/internal/common.h>
83 #include <vsmc/rng/urng.h>
84 #include <vsmc/rng/normal01.h>
85 #include <vsmc/rng/u01.h>
86 
87 #if VSMC_HAS_OPENCL_DOUBLE
88 
89 #define GAMMAK1_2x32 gammak1_2x32_53
90 #define GAMMAK1_2x32_INIT gammak1_2x32_53_init
91 #define GAMMAK1_2x32_RAND gammak1_2x32_53_rand
92 
93 #define GAMMAK1_2x64 gammak1_2x64_53
94 #define GAMMAK1_2x64_INIT gammak1_2x64_53_init
95 #define GAMMAK1_2x64_RAND gammak1_2x64_53_rand
96 
97 #define GAMMAK1_4x32 gammak1_4x32_53
98 #define GAMMAK1_4x32_INIT gammak1_4x32_53_init
99 #define GAMMAK1_4x32_RAND gammak1_4x32_53_rand
100 
101 #define GAMMAK1_4x64 gammak1_4x64_53
102 #define GAMMAK1_4x64_INIT gammak1_4x64_53_init
103 #define GAMMAK1_4x64_RAND gammak1_4x64_53_rand
104 
105 #else // VSMC_HAS_OPENCL_DOUBLE
106 
107 #define GAMMAK1_2x32 gammak1_2x32_24
108 #define GAMMAK1_2x32_INIT gammak1_2x32_24_init
109 #define GAMMAK1_2x32_RAND gammak1_2x32_24_rand
110 
111 #define GAMMAK1_4x32 gammak1_4x32_24
112 #define GAMMAK1_4x32_INIT gammak1_4x32_24_init
113 #define GAMMAK1_4x32_RAND gammak1_4x32_24_rand
114 
115 #endif // VSMC_HAS_OPENCL_DOUBLE
116 
117 #define VSMC_DEFINE_RNG_GAMMAK1(N, W, F, FT) \
118  typedef struct { \
119  FT c_shape, c_s2, c_s, c_d, c_b, c_si, c_c, c_q0; \
120  normal01_##N##x##W##_##F rnorm; \
121  } gammak1_##N##x##W##_##F;
122 
123 #define VSMC_DEFINE_RNG_GAMMAK1_INIT(N, W, F, FT) \
124  VSMC_STATIC_INLINE void gammak1_##N##x##W##_##F##_init ( \
125  gammak1_##N##x##W##_##F *rgamma, \
126  cburng##N##x##W##_rng_t *rng, FT shape) \
127  { \
128  rgamma->c_shape = shape; \
129  normal01_##N##x##W##_##F##_init(&(rgamma->rnorm), rng); \
130  \
131  if (shape <= 1) { \
132  rgamma->c_s2 = 0; \
133  rgamma->c_s = 0; \
134  rgamma->c_d = 0; \
135  rgamma->c_b = 0; \
136  rgamma->c_si = 0; \
137  rgamma->c_c = 0; \
138  rgamma->c_q0 = 0; \
139  return; \
140  } \
141  \
142  FT c_sqrt32 = 5.6568542494923802; \
143  FT c_s2 = shape - 0.5f; \
144  FT c_s = sqrt(c_s2); \
145  FT c_d = c_sqrt32 - c_s * 12; \
146  \
147  FT c_b, c_si, c_c; \
148  FT r = 1 / shape; \
149  FT q1 = 0.04166669; \
150  FT q2 = 0.02083148; \
151  FT q3 = 0.00801191; \
152  FT q4 = 0.00144121; \
153  FT q5 = -7.388e-5; \
154  FT q6 = 2.4511e-4; \
155  FT q7 = 2.424e-4; \
156  FT c_q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r + \
157  q2) * r + q1) * r; \
158  if (shape <= 3.686f) { \
159  c_b = 0.463f + c_s + 0.178f * c_s2; \
160  c_si = 1.235f; \
161  c_c = 0.195f / c_s - 0.079f + 0.16f * c_s; \
162  } else if (shape <= 13.022f) { \
163  c_b = 1.654f + 0.0076f * c_s2; \
164  c_si = 1.68f / c_s + 0.275f; \
165  c_c = 0.062f / c_s + 0.024f; \
166  } else { \
167  c_b = 1.77f; \
168  c_si = 0.75f; \
169  c_c = 0.1515f / c_s; \
170  } \
171  \
172  rgamma->c_s2 = c_s2; \
173  rgamma->c_s = c_s; \
174  rgamma->c_d = c_d; \
175  rgamma->c_b = c_b; \
176  rgamma->c_si = c_si; \
177  rgamma->c_c = c_c; \
178  rgamma->c_q0 = c_q0; \
179  }
180 
181 #define VSMC_DEFINE_RNG_GAMMAK1_RAND(N, W, F, FT) \
182  VSMC_STATIC_INLINE FT gammak1_##N##x##W##_##F##_rand ( \
183  gammak1_##N##x##W##_##F *rgamma, cburng##N##x##W##_rng_t *rng) \
184  { \
185  const FT c_shape = rgamma->c_shape; \
186  const FT c_s2 = rgamma->c_s2; \
187  const FT c_s = rgamma->c_s; \
188  const FT c_d = rgamma->c_d; \
189  const FT c_b = rgamma->c_b; \
190  const FT c_si = rgamma->c_si; \
191  const FT c_c = rgamma->c_c; \
192  const FT c_q0 = rgamma->c_q0; \
193  \
194  if (c_shape < 0) \
195  return -1; \
196  \
197  if (c_shape == 0) \
198  return 0; \
199  \
200  if (c_shape == 1) { \
201  FT u = u01_open_open_##W##_##F(cburng##N##x##W##_rand(rng)); \
202  return -log(u); \
203  } \
204  \
205  if (c_shape < 1) { \
206  const FT c_exp_m1 = 0.3678794411714423; \
207  FT b = 1 + c_exp_m1 * c_shape; \
208  while (true) { \
209  FT u = u01_open_open_##W##_##F(cburng##N##x##W##_rand(rng)); \
210  FT v = u01_open_open_##W##_##F(cburng##N##x##W##_rand(rng)); \
211  FT p = b * u; \
212  FT e = -log(v); \
213  if (p >= 1) { \
214  FT x = -log((b - p) / c_shape); \
215  if (e >= (1 - c_shape) * log(x)) \
216  return x * x; \
217  } else { \
218  FT x = exp(log(p) / c_shape); \
219  if (e >= x) \
220  return x * x; \
221  } \
222  } \
223  } \
224  \
225  FT t = normal01_##N##x##W##_##F##_rand(&(rgamma->rnorm), rng); \
226  FT x = c_s + t * 0.5f; \
227  if (t >= 0) \
228  return x * x; \
229  \
230  FT u = u01_open_open_##W##_##F(cburng##N##x##W##_rand(rng)); \
231  if (c_d * u <= t * t * t) \
232  return x * x; \
233  \
234  const FT a1 = 0.3333333; \
235  const FT a2 = -0.250003; \
236  const FT a3 = 0.2000062; \
237  const FT a4 = -0.1662921; \
238  const FT a5 = 0.1423657; \
239  const FT a6 = -0.1367177; \
240  const FT a7 = 0.1233795; \
241  \
242  if (x > 0) { \
243  FT v = t / (c_s + c_s); \
244  FT q; \
245  if (fabs(v) <= 0.25f) { \
246  q = c_q0 + 0.5f * t * t * \
247  ((((((a7 * v + a6) * v + a5) * v + a4) * \
248  v + a3) * v + a2) * v + a1) * v; \
249  } else { \
250  q = c_q0 - c_s * t + 0.25f * t * t + \
251  (c_s2 + c_s2) * log(1 + v); \
252  } \
253  \
254  if (log(1 - u) <= q) \
255  return x * x; \
256  } \
257  \
258  while (true) { \
259  FT e = u01_open_open_##W##_##F(cburng##N##x##W##_rand(rng)); \
260  FT u = u01_open_open_##W##_##F(cburng##N##x##W##_rand(rng)); \
261  e = -log(e); \
262  u = u + u - 1; \
263  FT t = u < 0 ? c_b - c_si * e : c_b + c_si * e; \
264  \
265  const FT c_tau = -0.71874483771719; \
266  if (t <= c_tau) \
267  continue; \
268  \
269  FT v = t / (c_s + c_s); \
270  FT q; \
271  if (fabs(v) <= 0.25f) { \
272  q = c_q0 + 0.5f * t * t * \
273  ((((((a7 * v + a6) * v + a5) * v + a4) * \
274  v + a3) * v + a2) * v + a1) * v; \
275  } else { \
276  q = c_q0 - c_s * t + 0.25f * t * t + \
277  (c_s2 + c_s2) * log(1 + v); \
278  } \
279  if (q <= 0) \
280  continue; \
281  if (c_c * fabs(u) > expm1(q) * exp(e - 0.5f * t * t)) \
282  continue; \
283  \
284  FT x = c_s + 0.5f * t; \
285  return x * x; \
286  } \
287  }
288 
290 VSMC_DEFINE_RNG_GAMMAK1(2, 32, 24, float)
292 VSMC_DEFINE_RNG_GAMMAK1(4, 32, 24, float)
293 
295 VSMC_DEFINE_RNG_GAMMAK1_INIT(2, 32, 24, float)
297 VSMC_DEFINE_RNG_GAMMAK1_INIT(4, 32, 24, float)
298 
300 VSMC_DEFINE_RNG_GAMMAK1_RAND(2, 32, 24, float)
302 VSMC_DEFINE_RNG_GAMMAK1_RAND(4, 32, 24, float)
303 
304 #if VSMC_HAS_OPENCL_DOUBLE
305 
307 VSMC_DEFINE_RNG_GAMMAK1(2, 32, 53, double)
309 VSMC_DEFINE_RNG_GAMMAK1(4, 32, 53, double)
310 
312 VSMC_DEFINE_RNG_GAMMAK1_INIT(2, 32, 53, double)
314 VSMC_DEFINE_RNG_GAMMAK1_INIT(4, 32, 53, double)
315 
317 VSMC_DEFINE_RNG_GAMMAK1_RAND(2, 32, 53, double)
319 VSMC_DEFINE_RNG_GAMMAK1_RAND(4, 32, 53, double)
320 
322 VSMC_DEFINE_RNG_GAMMAK1(2, 64, 53, double)
324 VSMC_DEFINE_RNG_GAMMAK1(4, 64, 53, double)
325 
327 VSMC_DEFINE_RNG_GAMMAK1_INIT(2, 64, 53, double)
329 VSMC_DEFINE_RNG_GAMMAK1_INIT(4, 64, 53, double)
330 
332 VSMC_DEFINE_RNG_GAMMAK1_RAND(2, 64, 53, double)
334 VSMC_DEFINE_RNG_GAMMAK1_RAND(4, 64, 53, double)
335 
336 #endif // VSMC_HAS_OPENCL_DOUBLE
337 
338 #endif // VSMC_RNG_GAMMAK1_H
VSMC_DEFINE_RNG_GAMMAK1_INIT
#define VSMC_DEFINE_RNG_GAMMAK1_INIT(N, W, F, FT)
Definition: gammak1.h:123
VSMC_DEFINE_RNG_GAMMAK1
#define VSMC_DEFINE_RNG_GAMMAK1(N, W, F, FT)
Definition: gammak1.h:117
VSMC_DEFINE_RNG_GAMMAK1_RAND
#define VSMC_DEFINE_RNG_GAMMAK1_RAND(N, W, F, FT)
Definition: gammak1.h:181
Generated by   doxygen 1.8.9.1