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



proc ranlux(nn);
   local j,k,m,i,x,c;
   j=10; k=24; m=2^24;
/* initial random sequence: x[1],x[2],...,x[k] */
   x=floor(rndu(k,1)*m);                        /* integers between 0 and m-1 */
/* Lagged Fibonacci with Subtract-With-Borrow */
   x=x|zeros(nn-k,1);
   c=0;
   i=k+1;
   do while i<=nn;
      x[i]=mod(x[i-j]-x[i-k]-c,m);
      if x[i-j]-x[i-k]-c<0;
         c=1;
      else;
         c=0;
      endif;
      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;