/* Trembling RAN3 U[0,1)                                   */
/* (C) Copyright 2005 Yosuke Amijima. All Rights Reserved. */
/* 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=ran3t(nn);
library pgraph;
graphset;
call hist(x,29);


proc ran3t(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;
/* add trembling factor */
   x=x-0.5+rndu(nn,1);    
   x=x/m;
   x=x.*(x.>=0 .and x.<1)+(1+x).*(x.<0)+(x-1).*(x.>=1);
   retp(x);
endp;

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