function [y] = randompowerlawints(n,alpha)

rand('state',sum(100 * clock))
x = rand(n,1);
if (alpha <= 1) 
        disp('warning: this function valid for alpha > 1. Returning zero vector');
        y = x*0;
end
disp(alpha)
y = (2^(alpha-1)*(1-x)).^(-1/(alpha-1));
y = round(y);