searle_gfaure.m

来自「very good matlab of Monte carlo」· M 代码 · 共 32 行

function Pn = searle_gfaure(n,d)
%Calculates the first n elements of the d-dimensional Generalised Faure sequence
%Method taken from Tezuka and Tokuyama (See paper by L'Ecuyer and Lemieux:Recent advances in Randomised... )

tic;
prime = primes(367);           %Determiming the smallest prime >=d
sprime = find(prime>=d);
p = prime(sprime(1)); 
maxb = ceil(log(n)/log(p));    %Upper bound for length of sequence

Pas = zeros(maxb,maxb);        %Calculating Pascal matrix
Pas(:,1) = 1;
for i = 2:maxb
   for j = 2:i
      Pas(i,j) = Pas(i-1,j) + Pas(i-1,j-1);
   end;
end;
Pas = mod(Pas,p);              %Ensuring elements are in the required field

bob = -[1:maxb];
rad_in = p.^bob;

y = [zeros(length(dec2bas([0:n-1],p)),maxb-size(dec2bas([0:n-1],p),2)),dec2bas([0:n-1],p)];
y=fliplr(y);
Pn(1,:) = rad_in*y';  %Since C^1 = A_1(P')^0 is an identity matrix
for l = 2:d
   y = mod(y*Pas*Pas',p);
   Pn(l,:) = rad_in*y';
end;
toc;