/* 2D Hammersley Sequence */
/* This is nothing but even intervals vs. van der Corput */
/* Notice: This version excludes the first zeros.        */
new; cls;
nn=1000;
x=Hammersley(nn);

library pgraph;
graphset;
_plctrl=-1;
xtics(0,1,0.1,0);
ytics(0,1,0.1,0);
xy(x[.,1],x[.,2]);





proc Hammersley(nn);
   local x1,bb,x,i,ff,n0,n1,r,n,U,index;
/* 1st dimension beginning from i/(nn+1) */
   x1=seqa(1,1,nn)/(nn+1);
/* Halton(2) or van der Corput */
   bb=2;
   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;
/* merge these */
   x=x1~x;
/* randomize the order as a pair */
   n=rows(x);
   U=x;
   x=zeros(n,2);
   i=1;
   do while i<=n-1;
      index=ceil((n-(i-1))*rndu(1,1));
      x[i,.]=U[index,.];
      if index==1;
         U=U[2:rows(U),.];
      elseif index==rows(U);
         U=U[1:rows(U)-1,.];
      else;
         U=U[1:(index-1) (index+1):rows(U),.];
      endif;
      i=i+1;
   endo;
   x[n,.]=U;
   retp(x);
endp;