/* Fluctuating RAN0 U[0,1)                                 */
/* (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=ran0f(nb,nn);    
library pgraph;
graphset;
call hist(x,29);



proc ran0f(nb,nn);
   local a,c,m,x,i;
   a=16807; c=0; m=2^31-1;
   x=zeros(nn,1);
/* seed */
   x[1]=floor(rndu(1,1)*(m-1))+1;         /* integer between 1 and m-1 */
/* LCG */
   i=1;
   do while i<=nn-1;
      x[i+1]=(a*x[i]+c)%m;
      i=i+1;
   endo;
/* add fluctuating factor */
   x=x-nb+floor((2*nb+1)*rndu(nn,1));
   x=mod(x,m);
   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;