function [X,U,S,V]=heavysing(n,a)
 B=sign(randn(n,n));         % Generates random signs 
 X=(rand(n,n).^(-1/a)-1).*B; % Generates heavy tails
 [U,S,V]=svd(X*n^(-1/a));    % Diagonalization, eigenvectors have unit 2-norm
 L=norm(V,Inf,'columns');    % Computes infinite norm of eigenvectors
 clf; plot(log(1+abs(S)),sqrt(n)*L,'ro'); hold on
 xlabel('log(1+|sigma|)'); ylabel('n^{1/2} |v|_{inf}|v|_2^{-1}')
 title(sprintf('SVD Xij = Rademacher x (Unif^{-1/alpha}-1) with alpha=%4.2f and n=%i',a,n))
 print(sprintf('heavysing_%4.2f_%i.png',a,n),'-dpng')
