%
%  macm 202 -- 31 mar 03 -- djm
%
%  shuffle.m:  produces a random permutation matrix
%              based on the GSR shuffling algorithm

function Pmatrix = shuffle(Ncards);

%  random binomial cut

cut = sum(floor(2*rand(Ncards,1)));

Ntop = cut;  Nbot = Ncards-cut;

%  build permutation matrix
%  - start with matrix of all zeros

Pmatrix = sparse(Ncards,Ncards);

%  use random choice of dropping from top or bot pile

for j= 1:Ncards
    if (rand < Ntop/(Ntop+Nbot))
        Pmatrix(j,cut   -Ntop+1) = 1;
        Ntop = Ntop-1;
    else
        Pmatrix(j,Ncards-Nbot+1) = 1;
        Nbot = Nbot-1;
    end
end