RE: Code to display LCG correlations in R

From: David Read

Dr. Suess,
 
Thanks for the suggestions, and here is some commented code you might find useful.
 
The first LCG is a bad one, but running it and the following plot alone should be 
a good visual representation of an LCG's period and the structure that comes about 
by the LCG algorithm.
 
The second LCG is a fantastic one with an immense period that was suggested by use 
by the U.S. Bureau of Standards in 1964. I believe the values it generates pass the 
DieHard suite, but don't quote me on that. At some point I want to learn how to run 
DieHard to see for myself, but I cannot get R to output the values in the proper format. 
(DieHard demands binary . . . )
 
If you take N down to say 500 the 2D graph looks completely random. At N=100000 and 
viewed in full screen, slightly tilted lines of structure appear. If the period were 
completed and R's silly circles could be turned into dots, then very nice stripes would 
appear with no breaks.
 
To see the 3D image, you will need to download the rgl package in R. I did this straight 
from cran using R itself. It was cake, the hard part was finding the rgl package. Then 
you have to load it into the session. Both are done with the pull down package menu from 
R, I can show you if this is new to you, but I doubt that.
 
I chose this 3D imaging package because it renders the image in OpenGL. This by default 
allows you to rotate and zoom in and out on your plot by clicking and dragging and 
rotating the mouse wheel. At N=100000 you can rotate the cubic cloud of values until 
lines appear. These are NOT rendering artifacts. Zoom in, rotate around, you can see that
there are no dots outside of the planes!! This is exactly what I anted to see, and your 
suggestion of using the current value on one axis, the one past value on another and 
the 2 past value on the third worked perfectly.
 
If you type in values for other LCG's you can see how the planes shift. Also, OpenGL can 
handle BILLIONS of points, so go head, crank N up and enjoy, but if your screen is small 
you'll see a lot of black. At N=1000000 the planes and other wonderful structures appear. 
Don't do this with R's standard 2D plot, it will crash R.
 
If I had to break codes, I would graph the encoded values, rotate the cube until planes 
appear, and note the angle I'm looking at them. I believe it could be shown that this 
angle where the planes appear is uniquely determined by one and only one set of M,C and D 
values. I like this approach, because you don't have to worry about the seed until you decode.
 
Anyway, that was a fun excursion into something more or less useless to my classes. But, 
since you answered my questions, this code may aid instruction.

 
Thanks again,
 
David Read



# Code to view correlation in pseudorandom numbers generated by a LCG
# David Read, California State University, Easy Bay, 3/30/09

# Note: This is a TERRIBLE LCG used only to show the correlation. 

N=26

Seed       = 44
Mvalue     = 26
Cvalue     = 27
Dvalue     = 3
v1         = numeric(N)
v1[1]      = Seed
for (i in 2:N)
	{Next = (Cvalue*Seed + Dvalue)%%Mvalue
       v1[i] = Next
       Seed = Next
      }

v2<-v1/Mvalue

plot(v2[1:(N-1)], v2[2:N], main="Correlation of Pseudorandom Numbers from a Linear Congruential Generator")

# Note: This LCG has passed, and still passes, many tests for "randomness" 
#       and is suggested for use by the U.S. Bureau of Standards in 1964. 

N=100000

Seed       = 762939453125                           	
Mvalue     = 247
Cvalue     = (27)+1
Dvalue     = 29741096258473
v1         = numeric(N)
v1[1]      = Seed
for (i in 2:N)
	{Next = (Cvalue*Seed + Dvalue)%%Mvalue
       v1[i] = Next
       Seed = Next
      }

v2<-v1/Mvalue

# View correlation in 2-dimensions. X-axis is v[i-1], Y-axis is v[i]. 
# Expand this image to full screen to see the "structure." 

plot(v2[1:(N-1)], v2[2:N], main="Correlation of Pseudorandom Numbers from a Linear Congruential Generator")

# View correlation in 3-dimensions. X-axis is v[i-2], Y-axis is v[i-1], 
# Z-axis is v[i]. 
# NOTE: You MUST install and load the R package rgl for this code to work. 
# The 3d image can be rotated by clicking and dragging. The mouse wheel zooms. 
# The "structure" can be seen by rotating the image until the empty planes can
# be seen. Again, full screen is best. 

plot3d(v2[1:(N-2)], v2[2:(N-1)], v2[3:N], main="Correlation of Pseudorandom Numbers from a Linear Congruential Generator")