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

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.

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?

Edited 6 Years Ago by oggiemc: n/a

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.


.

Edited 6 Years Ago by jephthah: n/a

Comments
Nice attept to elucidate it -- I couldn't begin to tell you if it's correct

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

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.

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.

Edited 6 Years Ago by jephthah: n/a

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

Edited 6 Years Ago by invisal: n/a

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.


.

Edited 6 Years Ago by jephthah: n/a

(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

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

Edited 6 Years Ago by invisal: n/a

This article has been dead for over six months. Start a new discussion instead.