/* (One-dimensional) Halton Sequence with arbitrary base                           */
/* This is also called as van der Corput in most books, but I follow Paul Wilmott. */
new; cls;
bb=2;             /* Any positive integer. If base b is 2, that is the same as van der Corput. */
nn=10000;
x=Haltonb(bb,nn);

library pgraph;
graphset;
call hist(x,29);





proc Haltonb(bb,nn);
   local x,i,ff,n0,n1,r;
   if bb<=1 or (bb-floor(bb))/=0;
      errorlog "ERROR: Base parameter bb should be a positive integer except one.";
      retp(".");
   endif;
   x=zeros(nn,1);
   i=1;
   do while i<=nn;
      n0=i;
      ff=1/bb;
      do while n0>0;
         n1=floor(n0/bb);
         r=n0-n1*bb;
         x[i]=x[i]+ff*r;
         ff=ff/bb;
         n0=n1;
      endo;
      i=i+1;
   endo;
   retp(x);
endp;