/* Fluctuating RAN2 U[0,1) (LCG(2^31-85,40014) with Bays-Durham + LCG(2^31-249,40692))%2^31-85 */
/*                     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;
nb=100;        /* +-nb of band width (if nb=0, original) */
nn=10000;                       /* Up to 10000 in Light version of GAUSS */
x=ran2f(nb,nn);
library pgraph;
graphset;
call hist(x,29);




proc ran2f(nb,nn);
   local a1,c1,m1,a2,c2,m2,k,table,i,x,temp,j,y,z;
   a1=40014; c1=0; m1=2^31-85; k=8;
   a2=40692; c2=0; m2=2^31-249;
/* x: LCG(2^31-85,40014) with Bays/Durham */
   table=zeros(k+1,1);
   /* seed */
   table[1]=floor(rndu(1,1)*(m1-1))+1;      /* integer between 1 and m1-1 */
   /* initial table */
   i=1;
   do while i<=k;
      table[i+1]=(a1*table[i]+c1)%m1;
      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]=(a1*table[k]+c1)%m1;
      else;
         x[i]=(a1*x[i-1]+c1)%m1;
      endif;
      j=(x[i]+1)%k;                        /* functional form */
      temp=x[i];
      x[i]=table[j+1];
      table[j+1]=temp;
      i=i+1;
   endo;
/* y: LCG(2^31-249,40692) */
   y=zeros(nn,1);
   /* seed */
   y[1]=floor(rndu(1,1)*(m2-1))+1;         /* integer between 1 and m2-1 */
   /* LCG */
   i=1;
   do while i<=nn-1;
      y[i+1]=(a2*y[i]+c2)%m2;
      i=i+1;
   endo;
/* combine (x+y)%(2^31-85) */
   z=(x+y)%m1;
/* add fluctuating factor */
   z=z-nb+floor((2*nb+1)*rndu(nn,1));
   z=mod(z,m1);
   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;