1 /*
2 * Copyright (c) 2015, The University of Oxford
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 * 1. Redistributions of source code must retain the above copyright notice,
8 * this list of conditions and the following disclaimer.
9 * 2. Redistributions in binary form must reproduce the above copyright notice,
10 * this list of conditions and the following disclaimer in the documentation
11 * and/or other materials provided with the distribution.
12 * 3. Neither the name of the University of Oxford nor the names of its
13 * contributors may be used to endorse or promote products derived from this
14 * software without specific prior written permission.
15 *
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
20 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
26 * POSSIBILITY OF SUCH DAMAGE.
27 */
28
29 #ifndef OSKAR_PRIVATE_RANDOM_HELPERS_H_
30 #define OSKAR_PRIVATE_RANDOM_HELPERS_H_
31
32 #include <oskar_global.h>
33 #include <math/oskar_cmath.h>
34 #include <math/private_random_generators.h>
35 #ifdef __CUDACC__
36 #include <cuda_runtime.h>
37 #endif
38
39 #ifdef __cplusplus
40 extern "C" {
41 #endif
42
43 OSKAR_INLINE
oskar_int_to_range_0_to_1_f(const unsigned long in)44 float oskar_int_to_range_0_to_1_f(const unsigned long in)
45 {
46 const float factor = 1.0f / (1.0f + 0xFFFFFFFFuL);
47 const float half_factor = 0.5f * factor;
48 return (in * factor) + half_factor;
49 }
50
51 OSKAR_INLINE
oskar_int_to_range_0_to_1_d(const unsigned long in)52 double oskar_int_to_range_0_to_1_d(const unsigned long in)
53 {
54 const double factor = 1.0 / (1.0 + 0xFFFFFFFFuL);
55 const double half_factor = 0.5 * factor;
56 return (in * factor) + half_factor;
57 }
58
59 OSKAR_INLINE
oskar_int_to_range_minus_1_to_1_f(const unsigned long in)60 float oskar_int_to_range_minus_1_to_1_f(const unsigned long in)
61 {
62 const float factor = 1.0f / (1.0f + 0x7FFFFFFFuL);
63 const float half_factor = 0.5f * factor;
64 return (((long)in) * factor) + half_factor;
65 }
66
67 OSKAR_INLINE
oskar_int_to_range_minus_1_to_1_d(const unsigned long in)68 double oskar_int_to_range_minus_1_to_1_d(const unsigned long in)
69 {
70 const double factor = 1.0 / (1.0 + 0x7FFFFFFFuL);
71 const double half_factor = 0.5 * factor;
72 return (((long)in) * factor) + half_factor;
73 }
74
75 OSKAR_INLINE
oskar_box_muller_f(unsigned long u0,unsigned long u1,float * f0,float * f1)76 void oskar_box_muller_f(unsigned long u0, unsigned long u1,
77 float* f0, float* f1)
78 {
79 float r;
80 #ifdef __CUDACC__
81 sincospif(oskar_int_to_range_minus_1_to_1_f(u0), f0, f1);
82 #else
83 float t = (float) M_PI;
84 t *= oskar_int_to_range_minus_1_to_1_f(u0);
85 *f0 = sinf(t);
86 *f1 = cosf(t);
87 #endif
88 r = sqrtf(-2.0f * logf(oskar_int_to_range_0_to_1_f(u1)));
89 *f0 *= r;
90 *f1 *= r;
91 }
92
93 OSKAR_INLINE
oskar_box_muller_d(unsigned long u0,unsigned long u1,double * f0,double * f1)94 void oskar_box_muller_d(unsigned long u0, unsigned long u1,
95 double* f0, double* f1)
96 {
97 double r;
98 #ifdef __CUDACC__
99 sincospi(oskar_int_to_range_minus_1_to_1_d(u0), f0, f1);
100 #else
101 double t = M_PI;
102 t *= oskar_int_to_range_minus_1_to_1_d(u0);
103 *f0 = sin(t);
104 *f1 = cos(t);
105 #endif
106 r = sqrt(-2.0 * log(oskar_int_to_range_0_to_1_d(u1)));
107 *f0 *= r;
108 *f1 *= r;
109 }
110
111 #ifdef __cplusplus
112 }
113 #endif
114
115 #endif /* OSKAR_PRIVATE_RANDOM_HELPERS_H_ */
116