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



proc mrg2(nn);
   local a,m,k,i,x;
   a=zeros(17,1);
   a[1]=2029258094; a[2]=250807368; a[3]=792743881; a[4]=659053828; a[5]=1554859284; a[6]=1198681177;
   a[7]=88961812; a[8]=1337122565; a[9]=109683638; a[10]=698262476; a[11]=442136888; a[12]=498352877;
   a[13]=184405303; a[14]=354728717; a[15]=949102210; a[16]=584030810; a[17]=1775088280;
   m=2^31-151;
   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;