## oggiemc 7

Hi everyone,

Ive been asked to do an assignment to test the randomness of a linear congruential generator..We've been given the code implemented in c..the problem is i dont really see how this code relates to the theoretical LCG equation i.e X[n+1] = (aX[n] + c) mod m..So i would be very appreciative if someone could explain the rand32() function in this code and how it implements the theoretical equation above..

many thanks

``````#include<stdio.h>
#include<stdint.h>

uint16_t rand16();
uint32_t rand32();

uint16_t n, userChosenValue = 0; //Change this for different seed values.
uint16_t numberOfValues = 10; //Change this depending on your needs.
uint32_t seed;
uint32_t mlcg,p,q;
uint64_t tmpseed;

int main(){

/* Calculate and print a series of 16 bit random numbers */
seed = (uint32_t)(userChosenValue + 1);

printf("16 Bit:\n\n");
for (n=0;n<numberOfValues;n++){
printf("%.4x\n",rand16());
}

/* Calculate and print a series of 32 bit random numbers */
seed = (uint32_t)(userChosenValue + 1);

printf("\n\n32 Bit:\n\n");
for (n=0;n<numberOfValues;n++){
printf("%.8x\t\t\n",rand32());
}

return 0;
}

/* Return the next 32 bit random number */
uint32_t rand32() {

tmpseed =  (uint64_t)33614U * (uint64_t)seed;
q = tmpseed;  /* low */
q = q >> 1;
p = tmpseed >> 32 ;           /* hi */
mlcg = p + q;
if (mlcg & 0x80000000) {
mlcg = mlcg & 0x7FFFFFFF;
mlcg++;
}
seed = mlcg;

return mlcg;
}

/* Return low 16 bits of next 32 bit random number */
uint16_t rand16() {
return (uint16_t)rand32();
}``````
commented: interesting +7

## jephthah 1,888

interesting question.

i believe this is a proper LCG. At least i can't say that it isn't. I don't fully understand how, but it looks like he's using shifts to perform modulus operation. see here.

im interested to know the details if you find out.

## oggiemc 7

Hi jephthah,

Got this explanation:

The below statement is aX[n] + c

q = q >> 1; <==== This is a part of aX[n]
p = tmpseed >> 32 ; /* hi */
mlcg = p + q;

And mod operator can be done like this.
x%y = if(x>y) {
z=x&(y-1);
z=z+1;
}

Above one is true only if y is 2^n-1

i dont really understand it..does it make any sense yo you?

## jephthah 1,888

well, this is some kind of voodoo guru stuff. pretty clever use of shifts and AND'ing to perform a modulus.

But i think i'm seeing the basics: this is the LCG as it is implemented by Apple's Carbon C-language API, according to the wiki link i posted above

given the LCG equation, `X[n+1] = (a*X[n] + c) mod m` ,

the Carbon LCG implements it with a = 16807, c = 0 (a special case), and m = 2^31 - 1. The choice of these numbers involve modular arithmetic discussions on Mersenne Primes and primitive roots modulo, and i have to admit at that point it gets beyond me pretty quickly.

but notice that in your program, as it stands, your 'tempseed' value = 33614 * 1, and then this is then shifted right 1 bit, which is the same as dividing by 2 ... so now you have 16807 (the 'q' variable in the program).

this next part is a little confusing:

to teh above 'q' value it adds the 'p' value, which is the upper 32 bits from the previous 64-bit seed. this initially is less than 32 bits, so that added value = 0, but as this iterative process develops, the seed varies across most of the range of a 64-bit value so the corresponding upper 32 bits (the 'p' variable) is a meaningful value.

now this part follows the explanation you received:

this summed value has its modulo 'm' taken. rather than actually perform a modulo operation, which is costly in cpu cycles, it does the tricky business with the AND and add one, which is fast. because the 'm' is a Mersenne prime, you can use that method that you described:

if x > y, and y = 2^n - 1, then z = x % y can be performed:

z = x & y
z++

if you want a simple example:

``````x = 45
y = 31   (31 = 2^5 - 1)

101101
&  011111
-----------
=  001101

+       1
=  001110 = 14

45 % 31 = 14``````

so the same applies in your case where y = 2^31 - 1 ... Hope this helped. i'm still a little shaky on this, so if anyone else has anything to add, please do.

Maybe Narue will stop by, i think she knows a good bit about PRNGs.

.

commented: Nice attept to elucidate it -- I couldn't begin to tell you if it's correct +4

## oggiemc 7

thanks jephthah, sorry only replying to you now..your explanation makes sense although the p and q addition is a little bit confusing.. can you explain "this initially is less than 32 bits, so that added value = 0"

## invisal 381

I think jephthah got it right and it is LCG in my opinion because:

``````tmpseed =  (uint64_t)33614U * (uint64_t)seed;	// 33614 * X[n]
q = tmpseed;					// 33614 * X[n] mod 2^32
q = q >> 1;					// 16807 * X[n] mod 2^32
p = tmpseed >> 32;				// (33614*X[n]) / (2^32) = (16807*X[n]) / (2^31)
mlcg = p + q;
// ((16807 * X[n]) mod 2^31) + (16807*X[n]) / (2^31)
// 16807 * X[n] + (16807*X[n]) / (2^31)  mod 2^31
// X[n] * ((16807 * 2^31 + 16807) / (2^31)) mod 2^31``````

so it has form of aX[n] + c mod m where
a = (16807 * 2^31 + 16807) / (2^31)
c = 0
m = 2^31

It is similar to Park–Miller random number generator

My opinion might be wrong.

## jephthah 1,888

the Apple Carbon C API generator, which this appears to be, is a Park-Miller PRNG.

note, the mod 'm' is `2^31 - 1 = 0x7FFFFFFF` ... the modulus operation is being performed at the end where in teh case where mlcg >= 2^31, the mlcg is ANDed with 0x7FFFFFFF and then mlcg++ ...

that bit of trickery IS a modulus operation for Mersenne primes, the set of which 2^31 - 1 is a member

to be honest though, i'm still not certain how the iterative procedure where the top 32bits are added to the lower 32 bits (mlcg = p + q) all fits together.

## invisal 381

To be honest, `mlcg = mlcg & 0x7FFFFFFF;` is the same as `mlcg = mlcg % 0x80000000` So, m cannot be `2^31 - 1`

## jephthah 1,888

sorry, that's incorrect.

In 1988, Park and Miller [1] suggested RNG with particular parameters n=2^31−1 = 2,147,483,647 (a Mersenne prime M31) and g=16,807 (a primitive root modulo M31), now known as MINSTD. Despite that MINSTD was later criticized by Marsaglia and Sullivan, [2] it is still in use today (in particular, in CarbonLib).

for another thing, you're forgetting the +1 increment which is required when using this method to perform the modulus using a Marsenne prime.

.

## invisal 381

(X[n] mod (2^31) )+1 is not the same as X[n] mod (2^31-1)

For example:
89 mod 31 = 27
(89 mod 32) + 1 = 26

## invisal 381

Look like I found the trick of the code.

If you want to find modulus of B-1, you can add the result of `A mod B` with A/B,
so `A mod B-1 = (A mod B)+ (int)(A/B)`

``````tmpseed = (uint64_t)33614U * (uint64_t)seed;
q = tmpseed; /* low */
q = q >> 1;
p = tmpseed >> 32 ; /* hi */
mlcg = p + q;``````

So the purpose of this code is `16807 * X[n] mod (2^31-1)` . The purpose of `p+q` is to make conversion from `(2^31) to (2^31)-1`

``````mlcg = mlcg & 0x7FFFFFFF;
mlcg++;``````

Since `p + q` cannot reach 0xF000000, even if `p + q` is greater than 08F000000, the result of `mlcg / 0x80000000` will always be 1. So you only need to +1 to make conversion from 2^31 to 2^31-1

sorry, that's incorrect.

It seem that my answer was not totally incorrect.

so it has form of aX[n] + c mod m where
a = (16807 * 2^31 + 16807) / (2^31)
c = 0
m = 2^31

If I make conversion from 2^31 to 2^31-1, Then I will have:
m = 2^31-1
c = 0
a = (16807 * 2^31 + 16807) / (2^31) - (16807)/(2^31)
a = (16807 * 2^31 - 16807) / (2^31)
a = (16807 * 2^31) / (2^31)
a = 16807