GammaLib

Here some Python code implementing an exponential cutoff power law random number generator. n is the number of random samples to be drawn, alpha the exponent and Lambda the cutoff value.

        x = []
        y=[]
        for i in range(10*n):
            y.append(xmin - (1./Lambda)*log(1.-random()))
        while True:
            ytemp=[]
            for i in range(10*n):
                if random()<pow(y[i]/float(xmin),-alpha):ytemp.append(y[i])
            y = ytemp
            x = x+y
            q = len(x)-n
            if q==0: break;

            if (q>0):
                r = range(len(x))
                shuffle(r)

                xtemp = []
                for j in range(len(x)):
                    if j not in r[0:q]:
                        xtemp.append(x[j])
                x=xtemp
                break;

            if (q<0):
                y=[]
                for j in range(10*n):
                    y.append(xmin - (1./Lambda)*log(1.-random()))

Here another code:

/*
Simulate power law with cut-off x^(-alpha)*exp(-lambda*x)
To simulate power-law with cutoff , one can generate an
exponentially distributed random number using the formula
above     (as k>0 and integer, so k start at 1)
and then accept or reject it with probability p or 1 - p 
respectively (i.e accept U1 <p or reject U1>p, and
 U1 is a uniform [0,1] random variable),
where p = (x/x_min)^(-alpha)  and  x_min=1.

http://www.santafe.edu/~aaronc/powerlaws/
*/
double rand_gen::powerlaw_dist(double alpha, double lamda)
{
  double x;
  do {
    x = exponential_dist(lamda);
  } while (pow(x,-1*alpha) < uniform_dist(0.,1.));
  return (x);
}

/*
To simulate exponential distribution exp(-lambda*x), the inverse
 method is used.
The cumulative distribution function for the exponential
distribution is:1-exp(-lambda*x). The inversion function
is  -log(1-U)/lambda. The simplified form is -log(U)/lambda,
where U is a uniform [0,1] random variable.
 http://cg.scs.carleton.ca/~luc/rnbookindex.html
*/
double rand_gen::exponential_dist(double lambda)
{
  return(-1*log(uniform_dist(0.,1.))/lambda);
}