I have a strange problem when calling runif in an Rcpp routine

I have a routine called test written in Rcpp which has a single argument, a function called g. In the routine I call another Rcpp routine called ru that returns a uniform [0,1]. (The same issue happens with other R random variable routines such as rbinom). I have a small loop that calls ru a few times and prints the result.

If that is all everything works fine. However, if I also make a call to g, the random numbers returned by ru are all the same! Even so g and ru have nothing in common.

Actually this happens only if g is also a Rcpp routine. If it is a R routine everything works as expected. In the example below both fC and fR simple return 0.0.

Here are the routines:

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector fC() { return Rcpp::wrap(0.0); }

Rcpp::NumericVector ru() { return Rcpp::runif(1); }                  

// [[Rcpp::export]]
void test(Rcpp::Function g) { 
      
  Rcpp::Rcout << "first ru" << "\n";
  for(int j=0; j<10; ++j) {
    Rcpp::NumericVector z = g() ;
    Rcpp::Rcout << ru()  <<" ";
  }
  
  Rcpp::Rcout<<"\n\n second ru \n";
  for(int j=0; j<10; ++j) Rcpp::Rcout<< ru() <<" ";
  Rcpp::Rcout<<"\n";
}

/***R
# and the R routine is
fR = function() {0.0}

# Calling the R function works just fine 
set.seed(1)
test(fR)

# but not when the Rcpp function is called
set.seed(1)
test(fC)
*/

This produces

> test(fR)
first ru
0.265509 0.372124 0.572853 0.908208 0.201682 0.89839 0.944675 0.660798 0.629114 0.0617863 

 second ru 
0.205975 0.176557 0.687023 0.384104 0.769841 0.497699 0.717619 0.991906 0.380035 0.777445 

> set.seed(1)

> test(fC)
first ru
0.265509 0.265509 0.265509 0.265509 0.265509 0.265509 0.265509 0.265509 0.265509 0.265509 

 second ru 
0.372124 0.572853 0.908208 0.201682 0.89839 0.944675 0.660798 0.629114 0.0617863 0.205975 

In the second call, the first loop with the call to g always returns the same number. So somehow a call to g seems to set a seed for runif if g is an Rcpp function, even though g does nothing.

Adding the GetRNGstate(); and PutRNGstate(); statements as described in the document does indeed to the trick! Thanks user20650.

The crucial part of the toy example now looks like this:

  for(int j=0;j<10;++j) {
     NumericVector z=g();
     GetRNGstate();
     Rcout<<ru()<<" ";
     PutRNGstate();
  }