// Version 1 

function [a] = binaire(p, n)
  x = n;
  i = 1;
  while x>0,
    a(i) = pmodulo( x , p );
    x = (x - a(i)) / p;
    i = i + 1;
  end;
endfunction

function [phi] = VanDerCorput(p, n )
  phi(1)=0;
  phi(2) = 1/p;
  for( i = 2 : n)
    a = binaire(p,i);
    x(1:length(a)) = 1/p;
    x = cumprod( x );
    phi(i+1) = sum( a .* x )
  end,
endfunction

// Version 2 

fin=-36;

function [nouvel_etat]= vdc(p,ancien_etat)
  nouvel_etat=ancien_etat;
  i=1;  
  while (nouvel_etat(i) == p-1) do
	nouvel_etat(i) = 0;
	i = i + 1;
  end;
  if(nouvel_etat(i) == fin) then
	nouvel_etat(i) = 1;
	nouvel_etat(i+1) = fin;
  else;
	 nouvel_etat(i) = nouvel_etat(i) + 1;
  end;
endfunction;

function [x] = Convert(etat)
  x = 0;
  i = 1;
  while (~(etat(i) == fin)) do
	x = x + etat(i) * PPuissance(i);
	i=i+1;
  end
endfunction


PPuissance=zeros(1,1000);

function [etat] = initialize_state(p,taille)
  etat(1:taille-1) = round((p-1)*rand(1,taille-1,"unif"));
  etat(taille)=fin;
endfunction

function [Vec] = vector_vdc(p,n)
  etat=initialize_state(p,100);
  PPuissance(1)=1/p;
  for i=2:1000
    PPuissance(i) = PPuissance(i - 1) / p;
  end;
  Vec=zeros(1,n);
  for i=1:n
	etat=vdc(p,etat);
	Vec(i)=Convert(etat);
  end
endfunction

p=2;
N=1000;
VDC=vector_vdc(p,N);
Rand=rand(1,N,"unif");

// pour voir la tranformation de Kakutani
X=VDC;
N=prod(size(X));
Y=VDC(2:N);
X=VDC(1:N-1);

xbasc();
isoview(0,1,0,1);
plot2d(X,Y,0)

for i=1:100
    conv_vdc(i) = abs(mean(VDC(1:i*100))-1/2);
    conv_ran(i)	= abs(mean(Rand(1:i*100))-1/2);
end
plot2d([conv_vdc,conv_ran]);

n=1000000;
n=1;
omega=2*%pi*n;
X=sin(omega*VDC);
Y=sin(omega*Rand);
for i=1:100
    conv_vdc(i) = mean(X(1:i*100));
    conv_ran(i)	= mean(Y(1:i*100));
end
plot2d([conv_vdc,conv_ran]);