% sampleRandomSings(n,k) calls testRandomSings(n) k times and find the
% frequency with which LU decomposition generates a zero vector in the
% last row and the frequency with which the rank is measured at n-1.
 
function [p,q] = sampleRandomSings(n,k)
 sz = 0;
 sr = 0;
 for i=1:k
    [z,r] = testRandomSing(n);
    sz = sz+z;
    sr = sr+r;
 end
 p = sz/k;
 q = sr/k;
end

% Construct a random singular n x n matrix by making the nth column a random
% weighted sum of the others. Check whether row echelon decomposition
% gives an exact zero vector in the final row. Test whether rank(u) gives
% the right rank (n-1)
 
function [t,r] = testRandomSing(n)
   a = rand(n,n-1);
   b = rand(n-1,1);
   a(:,n) = a*b;
   [l,u] = lu(a); % u is the row-echelon decomposition of a;
   t=(sum(u(n,:)==0)==n); % check that the last row of u is the zero vector
   r=(rank(u)==n-1);
end
