1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup arma_rng_cxx98
20 //! @{
21
22
23
24 class arma_rng_cxx98
25 {
26 public:
27
28 typedef unsigned int seed_type;
29
30 inline static void set_seed(const seed_type val);
31
32 arma_inline static int randi_val();
33 arma_inline static double randu_val();
34 inline static double randn_val();
35
36 template<typename eT>
37 inline static void randn_dual_val(eT& out1, eT& out2);
38
39 template<typename eT>
40 inline static void randi_fill(eT* mem, const uword N, const int a, const int b);
41
42 inline static int randi_max_val();
43 };
44
45
46
47 inline
48 void
set_seed(const arma_rng_cxx98::seed_type val)49 arma_rng_cxx98::set_seed(const arma_rng_cxx98::seed_type val)
50 {
51 std::srand(val);
52 }
53
54
55
56 arma_inline
57 int
randi_val()58 arma_rng_cxx98::randi_val()
59 {
60 #if (RAND_MAX == 32767)
61 {
62 // NOTE: this is a better-than-nothing solution
63 // NOTE: see also arma_rng_cxx98::randi_max_val()
64
65 u32 val1 = u32(std::rand());
66 u32 val2 = u32(std::rand());
67
68 val1 <<= 15;
69
70 return (val1 | val2);
71 }
72 #else
73 {
74 return std::rand();
75 }
76 #endif
77 }
78
79
80
81 arma_inline
82 double
randu_val()83 arma_rng_cxx98::randu_val()
84 {
85 return double( double(randi_val()) * ( double(1) / double(randi_max_val()) ) );
86 }
87
88
89
90 inline
91 double
randn_val()92 arma_rng_cxx98::randn_val()
93 {
94 // polar form of the Box-Muller transformation:
95 // http://en.wikipedia.org/wiki/Box-Muller_transformation
96 // http://en.wikipedia.org/wiki/Marsaglia_polar_method
97
98 double tmp1 = double(0);
99 double tmp2 = double(0);
100 double w = double(0);
101
102 do
103 {
104 tmp1 = double(2) * double(randi_val()) * (double(1) / double(randi_max_val())) - double(1);
105 tmp2 = double(2) * double(randi_val()) * (double(1) / double(randi_max_val())) - double(1);
106
107 w = tmp1*tmp1 + tmp2*tmp2;
108 }
109 while( w >= double(1) );
110
111 return double( tmp1 * std::sqrt( (double(-2) * std::log(w)) / w) );
112 }
113
114
115
116 template<typename eT>
117 inline
118 void
randn_dual_val(eT & out1,eT & out2)119 arma_rng_cxx98::randn_dual_val(eT& out1, eT& out2)
120 {
121 // make sure we are internally using at least floats
122 typedef typename promote_type<eT,float>::result eTp;
123
124 eTp tmp1 = eTp(0);
125 eTp tmp2 = eTp(0);
126 eTp w = eTp(0);
127
128 do
129 {
130 tmp1 = eTp(2) * eTp(randi_val()) * (eTp(1) / eTp(randi_max_val())) - eTp(1);
131 tmp2 = eTp(2) * eTp(randi_val()) * (eTp(1) / eTp(randi_max_val())) - eTp(1);
132
133 w = tmp1*tmp1 + tmp2*tmp2;
134 }
135 while( w >= eTp(1) );
136
137 const eTp k = std::sqrt( (eTp(-2) * std::log(w)) / w);
138
139 out1 = eT(tmp1*k);
140 out2 = eT(tmp2*k);
141 }
142
143
144
145 template<typename eT>
146 inline
147 void
randi_fill(eT * mem,const uword N,const int a,const int b)148 arma_rng_cxx98::randi_fill(eT* mem, const uword N, const int a, const int b)
149 {
150 if( (a == 0) && (b == RAND_MAX) )
151 {
152 for(uword i=0; i<N; ++i)
153 {
154 mem[i] = eT(std::rand());
155 }
156 }
157 else
158 {
159 const uword length = uword(b - a + 1);
160
161 const double scale = double(length) / double(randi_max_val());
162
163 for(uword i=0; i<N; ++i)
164 {
165 mem[i] = eT((std::min)( b, (int( double(randi_val()) * scale ) + a) ));
166 }
167 }
168 }
169
170
171
172 inline
173 int
randi_max_val()174 arma_rng_cxx98::randi_max_val()
175 {
176 #if (RAND_MAX == 32767)
177 return ( (32767 << 15) + 32767);
178 #else
179 return RAND_MAX;
180 #endif
181 }
182
183
184
185 //! @}
186