/* MRG1 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=mrg1(nn);
library pgraph;
graphset;
call hist(x,29);



proc mrg1(nn);
   local a,m,k,i,x;
   a={107374182,0,0,0,104480}; m=2^31-1;
   k=rows(a);
/* initial random sequence: x[1],x[2],...,x[k] */
   x=floor(rndu(k,1)*m);                        /* integers between 0 and m-1 */
/* multiple recursive generator */
   x=x|zeros(nn-k,1);
   i=k+1;
   do while i<=nn;
      x[i]=(a'rev(x[i-k:i-1]))%m;             /* x[i]=(a1*x_1+a2*x_2+...+ak*x_k)%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;