/* Trembling RAN1 U[0,1)     (RAN0 with Bays-Durham shuffling algorithm) */
/*                            Here, index j=(x+1)%k where k=8.           */
/* (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.               */
/* Notice: Initial seed of this version is based on rndu for simplicity. */
new; cls;
nn=10000;                       /* Up to 10000 in Light version of GAUSS */
x=ran1t(nn);
library pgraph;
graphset;
call hist(x,29);




proc ran1t(nn);
   local a,c,m,k,table,i,x,temp,j;
   a=16807; c=0; m=2^31-1; k=8;
   table=zeros(k+1,1);
/* seed */
   table[1]=floor(rndu(1,1)*(m-1))+1;      /* integer between 1 and m-1 */
/* initial table */
   i=1;
   do while i<=k;
      table[i+1]=(a*table[i]+c)%m;
      i=i+1;
   endo;        
   table=table[2:k+1];
/* Bays-Durham shuffling */
   x=zeros(nn,1);
   i=1;
   do while i<=nn-1;
      if i==1;
         x[i]=(a*table[k]+c)%m;
      else;
         x[i]=(a*x[i-1]+c)%m;
      endif;
      j=(x[i]+1)%k;                        /* functional form */
      temp=x[i];
      x[i]=table[j+1];
      table[j+1]=temp;
      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;