/* Trembling R250 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;
nn=10000;                      /* Up to 10000 in Light version of GAUSS */
x=r250t(nn);
library pgraph;
graphset;
call hist(x,29);



proc r250t(nn);
   local a,c,m,q,r,x,i,high,low,t;
   a=16807; c=0; m=2^31-1;
   q=127773; r=2836;
   x=zeros(nn,1);
/* seed */
   x[1]=floor(rndu(1,1)*(m-1))+1;                    /* integer between 1 and m-1 */
/* R250 */
   i=1;
   do while i<=nn-1;
      high=x[i]/q;
      low=x[i]%q;
      t=a*low-r*high;
      if t>=0;
         x[i+1]=t;
      else;
         x[i+1]=t+m;
      endif;
      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;