/* Sobol Sequence     */
/* Notice: Very Slow. */
new; cls;
nn=100;
dim=3;
x=Sobol(nn,dim);

library pgraph;
graphset;
pqgwin auto;
_plctrl=-1;
call hist(x[.,1],19);
call hist(x[.,2],19);
call hist(x[.,3],19);
xyz(x[.,1],x[.,2],x[.,3]);



proc Sobol(nn,dim);
   local x,lastq,Sseed,maxdim,seed,v,p,m,i,j,k,l,jj,newv,inc,rec,tempseed;
   if dim<2 or dim>40;
      errorlog "ERROR: Dimension number must be an integer where 2 <= dim <= 40.";
      retp(".");
   endif;
   x=zeros(nn,dim);
   lastq=zeros(1,dim);
   Sseed=0;
   maxdim=40;
   seed=2;  /* exclude seed=1 corresponding to zero vector at the beginning */
   do while seed<=nn+1;
      let v[40,30]=
          1  0  0  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  0  0  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  0  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  5  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1  1  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3  7  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3  9  9  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7 13  3  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  5 11 27  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

          1  3  5  1 15  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  7  3 29  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7  7 21  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  1  9 23 37  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3  5 19 33  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3 13 11  7  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  7 13 25  5  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  5 11  7 11  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  1  3 13 39  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1 15 17 63 13  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

          1  1  5  5  1 27 33  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3  3 25 17 115 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3 15 29 15 41  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1  7  3 23 79  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7  9 31 29 17  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  5 13 11  3 29  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1  9  5 21 119 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3  1 23 13 75  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3 11 27 31 73  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  7  7 19 25 105 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

          1  3  5  5 21  9  7  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  1 15  5 49 59  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  1  1  1 33 65  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  5 15 17 19 21  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  7 11 13 29  3  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7  5  7 11 113 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  5  3 15 19 61  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1  1  9 27 89  7  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3  7 31 15 45 23  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3  9  9 25 107 39 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
         ;
      p={1,   3,   7,  11,  13,  19,  25,  37,  59,  47,
        61,  55,  41,  67,  97,  91, 109, 103, 115, 131,
        193, 137, 145, 143, 241, 157, 185, 167, 229, 171,
        213, 191, 253, 203, 211, 239, 247, 285, 369, 299};
      i=1;
      v[i,1:30]=ones(1,30);
      i=2;
      do while i<=maxdim;
         m=sumc(p[i].>=(2^seqa(1,1,9)));
         j=p[i];
         inc=zeros(m,1);
         k=m;
         do while k>=1;
            jj=floor(j/2);
            inc[k]=(j/=2*jj);
            j=jj;
            k=k-1;
         endo;
         j=m+1;
         do while j<=30;
            newv=v[i,j-m];
            l=1;
            k=1;
            do while k<=m;
               l=2*l;
               if inc[k];
                  newv=eor(newv,l*v[i,j-k]);
               endif;
               k=k+1;
            endo;
            v[i,j]=newv;
            j=j+1;
         endo;
         i=i+1;
      endo;
      l=1;
      j=30-1;
      do while j>=1;
         l=2*l;
         v[1:dim,j]=v[1:dim,j]*l;
         j=j-1;
      endo;
      rec=1/(2*l);
      if seed==Sseed+1;
         l=low0(seed);
      elseif seed<=Sseed;
         Sseed=0;
         l=1;
         lastq[1:dim]=0;
         tempseed=Sseed;
         do while tempseed<=seed-1;
            l=low0(tempseed);
            i=1;
            do while i<=dim;
               lastq[i]=eor(lastq[i],v[i,l]);
               i=i+1;
            endo;
            tempseed=tempseed+1;
         endo;
         l=low0(seed);
      elseif Sseed+1<seed;
         tempseed=Sseed+1;
         do while tempseed<=seed-1;
            l=low0(tempseed);
            i=1;
            do while i<=dim;
               lastq[i]=eor(lastq[i],v[i,l]);
               i=i+1;
            endo;
            tempseed=tempseed+1;
         endo;
         l=low0(seed);
      endif;
      i=1;
      do while i<=dim;
         x[seed-1,i]=lastq[i]*rec;
         lastq[i]=eor(lastq[i],v[i,l]);
         i=i+1;
      endo;
      Sseed=seed;
      seed=seed+1;
   endo;
   retp(x);
endp;


proc low0(x);
   local n,z,y;
   if x==0;
      retp(0);
   else;
      n=maxc(log(x)/log(2))+1;
      z=reshape(2,n,rows(x));
      y=rev(recserrc(x,z));
      y=0|y;
      y=rev(y);
      y=seqa(1,1,rows(y)).*(y.==0);
      y=delif(y,y.==0);
      retp(y[1]);
   endif;
endp;

proc eor(a,b);
   local a2,b2,dim,y,x;
   a2=dectobin(a);
   b2=dectobin(b);
   dim=maxc(rows(a2)|rows(b2));
   if rows(a2)>rows(b2);
      b2=zeros(rows(a2)-rows(b2),1)|b2;
   elseif rows(a2)<rows(b2);
      a2=zeros(rows(b2)-rows(a2),1)|a2;
   endif;
   y=a2+b2;
   y=y.*(y./=2);
   x=bintodec(y);
   retp(x);
endp;

proc dectobin(x);           /* dec to bin */
   local n,z,y;
   if x==0;
      y=0;
   else;
      n=maxc(log(x)/log(2))+1;
      z=reshape(2,n,rows(x));
      y=rev(recserrc(x,z)); 
   endif;
   retp(y);
endp;

proc bintodec(y);           /* bin to dec */
   if rows(y)==1;
      retp(y);
   else;
      retp(polyeval(2,y));
   endif;
endp;






/* Another version of Sobol */
/* Notice: Faster but you need to set coefficient number vector. */

co={1,5,16,17};      /* Integers between 1 and 21 given coefficient matrix c. */
x=Sobolc(1000,co);

library pgraph;
graphset;
pqgwin auto;
_plctrl=-1;
title("another version of Sobol");
xyz(x[.,1],x[.,2],x[.,4]);




proc Sobolc(nn,no);
   local dim,m,c,vv,v,j,i,dm,r,a,bn,first,last,x;
   dim=rows(no);
   m=29;
   c={0 0 0 0 0 0 0 1 1,
      0 0 0 0 0 1 0 1 1,
      0 0 0 0 0 1 1 0 1,
      0 0 0 0 1 0 0 1 1,
      0 0 0 0 1 1 0 0 1,
      0 0 0 0 1 1 1 1 1,
      0 0 0 1 0 0 1 0 1,
      0 0 0 1 0 1 0 0 1,
      0 0 0 1 0 1 1 1 1,
      0 0 0 1 1 0 1 1 1,
      0 0 0 1 1 1 0 1 1,
      0 0 0 1 1 1 1 0 1,
      0 0 1 0 0 0 0 1 1,
      0 0 1 0 0 1 0 0 1,
      0 0 1 0 1 0 1 1 1,
      0 0 1 0 1 1 0 1 1,
      0 0 1 1 0 0 0 0 1,
      0 0 1 1 0 0 1 1 1,
      0 0 1 1 0 1 1 0 1,
      0 0 1 1 1 0 0 1 1,
      0 0 1 1 1 0 1 0 1};
   vv=zeros(dim,m);
   j=seqa(1,1,m);
   v=2^(m-j).*(2^j-1);
   dm=1;
   do while dm<=dim;
      a=find(c[no[dm],.]);
      first=a[1];
      r=cols(c)-first;
      vv[dm,.]=v';
      j=r+1;
      do while j<=m;
         vv[dm,j]=vv[dm,j-r]/2^r;
         i=1;
         do while i<=r-1;
            vv[dm,j]=eor(vv[dm,j],c[no[dm],first+i]*vv[dm,j-1]);
            i=i+1;
         endo;
         vv[dm,j]=eor(vv[dm,j],vv[dm,j-r]);
         j=j+1;
      endo;
      dm=dm+1;
   endo;
   x=zeros(dim,nn+1);
   x[.,1]=zeros(dim,1);
   i=2;
   do while i<=nn+1;
      bn=dectobin(i-1);
      a=find(bn);
      last=rows(bn)-a[rows(a)]+1;
      j=1;
      do while j<=dim;
         x[j,i]=eor(x[j,i-1],vv[j,last]);
         j=j+1;
      endo;
      i=i+1;
   endo;
   x=x/2^m;
   x=x[.,2:nn]';
   retp(x);
endp;
 
proc eor(a,b);
   local a2,b2,dim,y,x;
   a2=dectobin(a);
   b2=dectobin(b);
   dim=maxc(rows(a2)|rows(b2));
   if rows(a2)>rows(b2);
      b2=zeros(rows(a2)-rows(b2),1)|b2;
   elseif rows(a2)<rows(b2);
      a2=zeros(rows(b2)-rows(a2),1)|a2;
   endif;
   y=a2+b2;
   y=y.*(y./=2);
   x=bintodec(y);
   retp(x);
endp;

proc dectobin(x);           /* dec to bin */
   local n,z,y;
   if x==0;
      y=0;
   else;
      n=maxc(log(x)/log(2))+1;
      z=reshape(2,n,rows(x));
      y=rev(recserrc(x,z)); 
   endif;
   retp(y);
endp;

proc bintodec(y);           /* bin to dec */
   if rows(y)==1;
      retp(y);
   else;
      retp(polyeval(2,y));
   endif;
endp;

proc find(x);
   local i,y;
   x=vec(x');
   i=seqa(1,1,rows(x));
   y=(x./=0).*i;
   y=delif(y,y.==0);
   retp(y);
endp;