function [X,V,S]=heavy(n,a)
 B=sign(randn(n,n));         % Generates random signs 
 X=(rand(n,n).^(-1/a)-1).*B; % Generates heavy tails
 X=tril(X,-1)+tril(X)';      % Symmetrization
 [V,S]=eig(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+|lambda|)'); ylabel('n^{1/2} |v|_{inf}|v|_2^{-1}')
 title(sprintf('Xij = Xji = Rademacher x (Unif^{-1/alpha}-1) with alpha=%4.2f and n=%i',a,n))
 print(sprintf('heavy_%4.2f_%i.png',a,n),'-dpng')
