Why do these RNG's in C++ and R not produce similar results?

Question

Please excuse the disgustingly noobish nature of this post, but I have a question for those who program in C++ and R on their personal computer.

Question: Why are these random numbers produced from the two programs below not equal, and how do I resolve this issue?

1. Firstly, I suspect that I have misused the local function and the <<- operator in the R program.
2. Secondly, I suspect it may be a floating-accuracy issue. It's not immediately obvious to me how the two programs are different, so I don't know how to go about this problem.

I have tried casting all my calculations in C++ to double/float (even long double), and using fmod instead of the modulus operator %: different outputs again, but still not similar to the output in R. I do not know if it of any significant importance, but I want to add that I am compiling the C++ code using the G++ compiler.

Algorithm: The following algorithm can be used in any standard personal computer. It was proposed to use in parallel three word generators,

• m_k = 171 m_k-1 (mod 30269)
• m'_k = 172 m'_k-1 (mod 30307)
• m''_k = 172 m''_k-1 (mod 30323)

and to use as pseudorandom numbers the fractional parts

• g_k = {m_k / 30269 + m'_k / 30307 + m''_k / 30323}

I have used the initial values m₀ = 5, m'₀ = 11, and m''₀ = 17.

Programs: I have the following program in C++:

//: MC:Uniform.cpp
// Generate pseudo random numbers uniformly between 0 and 1
#include <iostream>
#include <math.h> // For using "fmod()"
using namespace std;

float uniform(){
    // A sequence of initial values
    static int x = 5;
    static int y = 11;
    static int z = 17;

    // Some integer arithmetic required
    x = 171 * (x % 177) - 2 * (x / 177);
    y = 172 * (x % 176) - 35 * (y / 176);
    z = 170 * (x % 178) - 63 * (z / 178);

    /* If both operands are nonnegative then the
    remainder is nonnegative; if not, the sign of
    the remainder is implementation-defined. */
    if(x < 0)
        x = x + 30269;
    if(y < 0)
        y = y + 30307;
    if(z < 0)
        z = z + 30323;

    return fmod(x / 30269. + y / 30307. + z / 30323., 1.);
}

int main(){
    // Print 5 random numbers
    for(int i = 0; i < 5; i++){
        cout << uniform() << ", ";
    }               
}///:~

The program exites with code and outputs the following:

0.686912, 0.329174, 0.689649, 0.753722, 0.209394,

I also have a program in R, that looks like the following:

## Generate pseudo random numbers uniformly between 0 and 1
uniform <- local({
    # A sequence of initial values
    x = 5
    y = 11
    z = 17

    # Use the <<- operator to make x, y and z local static 
    # variables in R.
    f <- function(){
        x <<- 171 * (x %% 177) - 2 * (x / 177)
        y <<- 172 * (y %% 176) - 35 * (y / 176)
        z <<- 170 * (z %% 178) - 63 * (z / 178)

        return((x / 30269. + y / 30307. + z / 30323.)%%1.)
    }
})

# Print 5 random numbers
for(i in 1:5){
    print(uniform())
}

This program exites with code as well and produces the output

[1] 0.1857093
[1] 0.7222047
[1] 0.05103441
[1] 0.7375034
[1] 0.2065817

Any suggestions are appreciated, thanks in advance.

Answer 1

You need a few more %/%'s (integer division) in your R code. Remember that numeric variables in R are floating-point, not integer, by default; so / will do ordinary division with a non-integral quotient. You've also left out the part where you deal with negative x/y/z.

f <- function(){
    x <<- 171 * (x %% 177) - 2 * (x %/% 177)
    y <<- 172 * (y %% 176) - 35 * (y %/% 176)
    z <<- 170 * (z %% 178) - 63 * (z %/% 178)

    if(x < 0)
        x <<- x + 30269;
    if(y < 0)
        y <<- y + 30307;
    if(z < 0)
        z <<- z + 30323;

    return((x / 30269. + y / 30307. + z / 30323.)%%1)
}

After making those changes, there doesn't seem to be anything seriously wrong with the result. A quick histogram of 100000 random draws looks very uniform, and there's no autocorrelation I can find. Still doesn't match your C++ result though....