/* Fluctuating RANECU U[0,1)                               */
/* (C) Copyright 2005 Yosuke Amijima. All Rights Reserved. */
/* Notice: Initial seeds of this version are 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=ranecuf(nb,nn);
library pgraph;
graphset;
call hist(x,29);




proc ranecuf(nb,nn);
   local a1,a2,c1,c2,m1,m2,m,x1,x2,x,i;
   a1=40014; c1=0; m1=2^31-85;
   a2=40692; c2=0; m2=2^31-249;
   m=2^31-85;
   x1=zeros(nn,1);
   x2=zeros(nn,1);
   x=zeros(nn,1);
/* seed */
   x1[1]=floor(rndu(1,1)*(m1-1))+1;                    /* integer between 1 and m1-1 */
   x2[1]=floor(rndu(1,1)*(m2-1))+1;                    /* integer between 1 and m2-1 */
/* Combine two LCG's and then take mod by m */
   x[1]=(x1[1]+x2[1])%m;
   i=1;
   do while i<=nn-1;
      x1[i+1]=(a1*x1[i]+c1)%m1;
      x2[i+1]=(a2*x2[i]+c2)%m2;
      x[i+1]=(x1[i+1]+x2[i+1])%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;