% This program generates triplets (U1, U2, U3) of iid Uniform[0,1] variates
% according to the linear congruential generator RANDU and the pseudorandom
% number generator that MATLAB uses.


Z0 = 1; %  Seed for RANDU
m = 2^(31);
a = 2^(16)+3;

% First we generate numbers in the interval [0,1]

Z = Z0;
U = Z0/m;
for i = 1:300
    newZ = mod(a*Z(i),m);
    newU = newZ/m;
    Z = [Z; newZ];
    U = [U; newU];
end
figure;
plot(U,0,'xr')
title('RANDU random number generator')
xlabel('Random numbers (Ui)')

matlabU = rand(1,300);
figure;
plot(matlabU,0,'xb')
title('MATLAB random number generator')
xlabel('Random numbers (Ui)')

% Next we generate 3-dim vectors in the cube [0,1]x[0,1]x[0,1]

for i = 300:2100
    newZ = mod(a*Z(i),m);
    newU = newZ/m;
    Z = [Z; newZ];
    U = [U; newU];
end
U1 = [];
U2 = [];
U3 = [];
for j = 1:floor(length(U)/3)
    U1 = [U1; U(3*(j-1)+1)];
    U2 = [U2; U(3*(j-1)+2)];
    U3 = [U3; U(3*(j-1)+3)];
end
figure;
scatter3(U1, U2, U3,'r','filled')
title('RANDU random number generator')
xlabel('Ui')
ylabel('U(i+1)')
zlabel('U(i+2)')

matlabU = rand(1,2100);
matlabU1 = [];
matlabU2 = [];
matlabU3 = [];
for j = 1:floor(length(U)/3)
    matlabU1 = [matlabU1; matlabU(3*(j-1)+1)];
    matlabU2 = [matlabU2; matlabU(3*(j-1)+2)];
    matlabU3 = [matlabU3; matlabU(3*(j-1)+3)];
end
figure;
scatter3(matlabU1, matlabU2, matlabU3,'b','filled')
title('MATLAB random number generator')
xlabel('Ui')
ylabel('U(i+1)')
zlabel('U(i+2)')
