1- // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
2- //
1+
32// Copyright (C) 2013-2014 Conrad Sanderson
43// Copyright (C) 2013-2014 NICTA (www.nicta.com.au)
5- //
4+ //
65// This Source Code Form is subject to the terms of the Mozilla Public
76// License, v. 2.0. If a copy of the MPL was not distributed with this
87// file, You can obtain one at http://mozilla.org/MPL/2.0/.
98//
10- // This file is based on Conrad's default generators and as such licensed under both
9+ // This file is based on Conrad's default generators and as such licensed under both
1110// the MPL 2.0 for his as well as the GNU GPL 2.0 or later for my modification to it.
1211
1312// Copyright (C) 2014 Dirk Eddelbuettel
2827// along with RcppArmadillo. If not, see <http://www.gnu.org/licenses/>.
2928
3029
31- // NB This files use R's uniform generator and can be compiled only when the R
30+ // NB This files use R's uniform generator and can be compiled only when the R
3231// headers are available as is the case for RcppArmadillo.
3332//
3433// Also note that you MUST set / reset the R RNG state. When using RcppArmadillo
4039
4140class arma_rng_alt {
4241public:
43-
42+
4443 typedef unsigned int seed_type;
45-
44+
4645 inline static void set_seed (const seed_type val);
47-
46+
4847 arma_inline static int randi_val ();
4948 arma_inline static double randu_val ();
5049 inline static double randn_val ();
51-
50+
5251 template <typename eT>
5352 inline static void randn_dual_val (eT& out1, eT& out2);
54-
53+
5554 template <typename eT>
5655 inline static void randi_fill (eT* mem, const uword N, const int a, const int b);
57-
56+
5857 inline static int randi_max_val ();
5958};
6059
@@ -72,7 +71,7 @@ inline void arma_rng_alt::set_seed(const arma_rng_alt::seed_type val) {
7271}
7372
7473arma_inline int arma_rng_alt::randi_val () {
75- return ::Rf_runif (0 , RAND_MAX); // std::rand();
74+ return static_cast < int >( ::Rf_runif (0 , RAND_MAX) ); // std::rand();
7675}
7776
7877arma_inline double arma_rng_alt::randu_val () {
@@ -84,43 +83,43 @@ inline double arma_rng_alt::randn_val() {
8483 // polar form of the Box-Muller transformation:
8584 // http://en.wikipedia.org/wiki/Box-Muller_transformation
8685 // http://en.wikipedia.org/wiki/Marsaglia_polar_method
87-
86+
8887 double tmp1;
8988 double tmp2;
9089 double w;
91-
90+
9291 do {
9392 tmp1 = double (2 ) * double (::Rf_runif (0 , 1 )) - double (1 );
9493 tmp2 = double (2 ) * double (::Rf_runif (0 , 1 )) - double (1 );
9594 // tmp1 = double(2) * double(std::rand()) * (double(1) / double(RAND_MAX)) - double(1);
9695 // tmp2 = double(2) * double(std::rand()) * (double(1) / double(RAND_MAX)) - double(1);
97-
96+
9897 w = tmp1*tmp1 + tmp2*tmp2;
9998 } while ( w >= double (1 ) );
100-
99+
101100 return double ( tmp1 * std::sqrt ( (double (-2 ) * std::log (w)) / w) );
102101}
103102
104- template <typename eT>
103+ template <typename eT>
105104inline void arma_rng_alt::randn_dual_val (eT& out1, eT& out2) {
106105 // make sure we are internally using at least floats
107106 typedef typename promote_type<eT,float >::result eTp;
108-
107+
109108 eTp tmp1;
110109 eTp tmp2;
111110 eTp w;
112-
111+
113112 do {
114113 tmp1 = eTp (2 ) * eTp (::Rf_runif (0 , RAND_MAX)) * (eTp (1 ) / eTp (RAND_MAX)) - eTp (1 );
115114 tmp2 = eTp (2 ) * eTp (::Rf_runif (0 , RAND_MAX)) * (eTp (1 ) / eTp (RAND_MAX)) - eTp (1 );
116115 // tmp1 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
117116 // tmp2 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
118-
117+
119118 w = tmp1*tmp1 + tmp2*tmp2;
120119 } while ( w >= eTp (1 ) );
121-
120+
122121 const eTp k = std::sqrt ( (eTp (-2 ) * std::log (w)) / w);
123-
122+
124123 out1 = eT (tmp1*k);
125124 out2 = eT (tmp2*k);
126125}
@@ -131,14 +130,14 @@ template<typename eT>
131130inline void arma_rng_alt::randi_fill (eT* mem, const uword N, const int a, const int b) {
132131 if ( (a == 0 ) && (b == RAND_MAX) ) {
133132 for (uword i=0 ; i<N; ++i) {
134- mem[i] = ::Rf_runif (0 , RAND_MAX);
133+ mem[i] = static_cast <eT>( ::Rf_runif (0 , RAND_MAX)); // std::rand( );
135134 // mem[i] = std::rand();
136135 }
137136 } else {
138137 const uword length = b - a + 1 ;
139-
138+
140139 const double scale = double (length) / double (RAND_MAX);
141-
140+
142141 for (uword i=0 ; i<N; ++i) {
143142 mem[i] = (std::min)( b, (int ( double (::Rf_runif (0 ,RAND_MAX)) * scale ) + a) );
144143 // mem[i] = (std::min)( b, (int( double(std::rand()) * scale ) + a) );
0 commit comments