George Marsaglia wrote a number of RNGs, either by himself, or with others. Here is one of them:
.
/***********************************************************************/
/* */
/* G. Marsaglia and W. W. Tsang, */
/* "The 64-bit Universal RNG". */
/* Generates double precision floating-point random numbers. */
/* */
/* It produces uniform [0,1) variates directly, as 64-bit */
/* floating-point numbers, without conversion from integers. */
/* The resultant numbers have a very long period (2**202 or 10**61) */
/* with excellent randomness. The generator produces the same random */
/* numbers in IEEE machines, regardless of computer language used. */
/* */ /***********************************************************************/

/* Translated to PL/I from Fortran, by R. A. Vowels, 23 November 2005. */

(subscriptrange, fixedoverflow, size):
main: proc options (main);
declare x float binary (53), i fixed binary (31);
declare U(97) static float binary (53) external;

/* George Marsaglia & Wai Wan Tsang's 64-bit pseudo-random number generator */ uni64: procedure returns (float binary (53)) options(reorder);
declare r float binary (53) static initial
((9007199254740881.0e0 / 9007199254740992.e0));
declare d float binary (53) static initial
((362436069876.0e0 / 9007199254740992.0e0));
declare c float binary (53) static initial (0),
i fixed binary (7) static initial (97),
j fixed binary (7) static initial (33),
x float binary (53);
declare U(97) float binary (53) external;

x=U(i)-U(j); if x<0 then x=x+1.0; U(i)=x;
i = i - 1; if i=0 then i=97; j = j - 1; if j=0 then j=97;
c=c-d; if c<0 then c=c+r;
x=x-c; if x<0 then return (x+1.); return (x);
end uni64;

/* A two-seed function for filling the static array U(97) one bit at a time. */ fillU: procedure (seed1, seed2) options(reorder);
declare (seed1, seed2) fixed binary (31);
declare (s, t) float binary (53),
(x, y) fixed binary (31),
(i, j) fixed binary (7);
declare U(97) float binary (53) external;

x=seed1; y=seed2;
do i=1 to 97;
s=0.0; t=0.5;
do j=1 to 53;
(nofixedoverflow, nosize):
x=REM(6969*x, 65543);
(nofixedoverflow, nosize):
y=REM(8888*y, 65579);
if IAND(IEOR(x,y), 32)>0 then s=s+t;
t=.5*t;
end;
U(i)=s;
end;
end fillU;

call fillU (123456789,987654321);

do i=1 to 10000000; x=uni64(); end;
put skip list (' x x*9007199254740992.0');
do i=1 to 5;
x = uni64();
put skip edit(x, x*9007199254740992.0E0) (f(20,17), f(18));
end;
end main;

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)