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



proc cmrg(nn);
   local a,b,m1,m2,x,y,z,k,i;
   a={0,63308,-183326}; m1=2^31-1;
   b={86098,0,-539608}; m2=2145483479;
   k=3;
/* initial random sequences */
   x=floor(rndu(k,1)*m1);                        /* integers between 0 and m1-1 */
   y=floor(rndu(k,1)*m2);                        /* integers between 0 and m2-1 */
   z=zeros(k,1);
   i=1;
   do while i<=k;
      z[i]=mod((x[i]-y[i]),m1);
      i=i+1;
   endo;
/* multiple recursive generator */
   x=x|zeros(nn-k,1);
   y=y|zeros(nn-k,1);
   z=z|zeros(nn-k,1);
   i=k+1;
   do while i<=nn;
      x[i]=mod((a'rev(x[i-k:i-1])),m1);           /* x[i]=(a1*x_1+a2*x_2+...+ak*x_k)%m1 */
      y[i]=mod((b'rev(y[i-k:i-1])),m2);           /* y[i]=(b1*y_1+b2*y_2+...+bk*y_k)%m1 */
      z[i]=mod((x[i]-y[i]),m1);
      i=i+1;
   endo;       
   z=z/m1;            /* If you want to avoid 0's, replace this line by z=(z+1)/(m1+1); */
   retp(z);
endp;

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