linear congruential generator code

[oggiemc]

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(); 
}

[jephthah]

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]

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
In your case , m=7FFFFFFF

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

[jephthah]

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.

.

[oggiemc]

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"

Thanks for your input

[invisal]

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]

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]

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

[jephthah]

sorry, that's incorrect.

from your own link that you posted:

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]

(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]

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