/* RAN3 U[0,1) */
/* Notice: Initial seeds of this version, which are all kept and used for the results, */
/*         are based on rndu for simplicity.                                           */
new; cls;
nn=10000;                      /* Up to 10000 in Light version of GAUSS */
x=ran3(nn);
library pgraph;
graphset;
call hist(x,29);


proc ran3(nn);
   local j,k,m,i,x;
   j=33; k=55; m=10^9;
/* initial random sequence: x[1],x[2],...,x[k] */
   x=floor(rndu(k,1)*m);                        /* integers between 0 and m-1 */
/* Lagged Fibonacci */
   x=x|zeros(nn-k,1);
   i=k+1;
   do while i<=nn;
      x[i]=mod(x[i-k]-x[i-j],m);
      i=i+1;
   endo;       
   x=x/m;         /* If you want to avoid 0's, replace this line by x=(x+1)/(m+1); */
   retp(x);
endp;

proc mod(x,m);
   retp( (x.>=0).*(x%m)+(x.<0).*((x+m.*ceil(abs(x)/m))%m) );
endp;