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



proc ranmar(nn);
   local i97,j97,x,ix,ij,kl,i,j,k,l,u,ii,jj,s,t,m,c,cd,cm,rval;
/* seeds */
   ij=floor(rndu(1,1)*(31328+1));      /* integer between 0 and 31328 */
   kl=floor(rndu(1,1)*(30081+1));      /* integer between 0 and 30081 */
/* initial u's */
   i97=97;
   j97=33;
   i=ij/177%177+2;
   j=ij%177+2;
   k=kl/169%178+1;
   l=kl%169;
   u=zeros(97,1);
   ii=1;
   do while ii<=97;
      s=0;
      t=0.5;
      jj=1;
      do while jj<=24;
         m=(i*j%179)*k%179;
         i=j;
         j=k;
         k=m;
         l=(53*l+1)%169;
         if (l*m%64)>=32;
            s=s+t;
         endif;
         t=0.5*t;
         jj=jj+1;
      endo;
      u[ii]=s;
      ii=ii+1;
   endo;
/* main */
   x=zeros(nn,1);
   ix=1;
   do while ix<=nn;
      c=362436/16777216;
      cd=7654321/16777216;
      cm=16777213/16777216;
      rval=u[i97]-u[j97];
      if rval<0;
         rval=rval+1;
      endif;
      u[i97]=rval;
      i97=i97-1;
      if i97==0;
         i97=97;
      endif;
      j97=j97-1;
      if j97==0;
         j97=97;
      endif;
      c=c-cd;
      if c<0;
         c=c+cm;
      endif;
      rval=rval-c;
      if rval<0;
         rval=rval+1;
      endif;
      x[ix]=rval;
      ix=ix+1;
   endo;
   retp(x);
endp;