Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Examples of generators of pseudorandom numbers

Linear congruent generator

Pseudorandom numbers are generated by

(Subscript[x, i+1])=(a (Subscript[x, i])+c) mod m

According to the choice of a, c, m, and (Subscript[x, 0]) we get pseudorandom number generators of different quality

myRandom[n_,dim_,par_,seed_]:=Module[{out=ConstantArray[0,{n,dim}]},
xRandom=seed;
Do[
Do[
xRandom=Mod[par[[1]] xRandom+par[[2]],par[[3]]];
out[[i,j]]=xRandom,
{j,1,dim}],
{i,1,n}
];
Return[out]
];

Generator RANDU

Expected value, variance and autocorrelation test:

n=10^5;
points=N[Flatten[myRandom[n,1,{(2^16)+3,0,(2^31)},100]/(2^31)]];
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and variance should be: ", N[1/2],", ",N[1/12]];
Print["Expected value and variance obtained:  ", expectedValue,", ",variance];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

But points in 3D are in 15 planes:

n=10^4;
points=myRandom[n,3,{(2^16)+3,0,(2^31)},1]/(2^31);
ListPointPlot3D[points,PlotStyle->PointSize[Small],PlotRange->{{0,1},{0,1},{0,1}}]

Generator from Apple CarbonLib

Expected value, variance and autocorrelation test:

n=10^5;
m=(2^31)-1;
points=N[Flatten[myRandom[n,1,{16807,0,m},100]/m]];
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and variance should be: ", N[1/2],", ",N[1/12]];
Print["Expected value and variance obtained:  ", expectedValue,", ",variance];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

ListPointPlot3D[myRandom[10000,3,{16807,0,m},1]/m,PlotStyle->PointSize[Small],PlotRange->{{0,1},{0,1},{0,1}}]

Mathematica's LCG generator

Expected value, variance and autocorrelation test:

n=10^5;
points=RandomReal[{0,1},n];
expectedValue=Total[points]/n ;
variance=Total[points^2]/n-expectedValue^2 ;
Print["Expected value and variance should be: ", N[1/2],", ",N[1/12]];
Print["Expected value and variance obtained:  ", expectedValue,", ",variance];
kmax=10; corr=ConstantArray[1.0,kmax+1];
Do[
corr[[k+1]]={k,(points[[1;;n-k]].points[[k+1;;n]]/(n-k)-expectedValue^2)/variance},
{k,0,kmax}
];
ListPlot[corr,PlotRange->All]

SeedRandom[1,Method->"Congruential"];
ListPointPlot3D[RandomReal[{0,1},{10000,3}],PlotStyle->PointSize[Small],PlotRange->{{0,1},{0,1},{0,1}}]