function [x b]=randpl(a,L,varargin)
% Random number with power-law distribution.
%
% x=randpl(a,L,dims); x is a list of dimensions dims drawn from the
% distribution (1+a)x^a/(L^(a+1)-1), a normalized power-law with min value
% 1 and max value L. abs(a) must be > 0.

% revision history:
% 02/23/98 Mark D. Shattuck <mds>   randpl.m

u=rand(varargin{:});
b=a+1;
x=((L^b-1)*u+1).^(1/b);