Stable distribution sampler in Julia -- forced to have scale of 1 instead of root 2.

stable.jl

# Utilities needed in rand()

# (exp(x) - 1) / x
function d2(z::Real)
	const p1 = 0.840066852536483239e3
	const p2 = 0.200011141589964569e2
	const q1 = 0.168013370507926648e4
	const q2 = 0.18013370407390023e3
	const q3 = 1.0

	if abs(z) > 0.1
		return (exp(z) - 1) / z
	else
		zz = z * z
		pv = p1 + zz * p2
		return 2 * pv / (q1 + zz * (q2 + zz * q3) - z * pv)
	end
end

# tan(x) / x
function tan2(xarg::Real)
	const p0 = 0.129221035e3
	const p1 = -0.887662377e1
	const p2 = 0.528644456e-1
	const q0 = 0.164529332e3
	const q1 = -0.451320561e2
	const q2 = 1.0
	const piby4 = π / 4

	x = abs(xarg)
	if x > piby4
		return tan(xarg) / xarg
	else
		x /= piby4
		xx = x * x
		return (p0 + xx * (p1 + xx * p2)) /
		         (piby4 * (q0 + xx * (q1 + xx * q2)))
	end
end

immutable Stable <: ContinuousUnivariateDistribution
	alpha::Float64 # Characteristic exponent
	beta::Float64 # Skewness in standard parameterization
	bprime::Float64 # Skewness in revised parameterization
	gamma::Float64 # Scale
	delta::Float64 # Location
	function Stable(a::Real, b::Real, g::Real, d::Real)
		0 < a <= 2 || throw(ArgumentError("alpha must be in (0, 2]"))
		abs(b) <= 1 || throw(ArgumentError("beta must be in [-1, 1]"))
		g > 0 || throw(ArgumentError("gamma must be non-negative"))
		phi0 = pi / 2 * b * ((1 - abs(1 - a)) / a) # Always zero for a = 2
		bprime = a == 1.0 ? b : -tan(pi / 2 * (1 - a)) * tan(a * phi0)
		new(float64(a), float64(b), float64(bprime), float64(g), float64(d))
	end
end

# 2-parameter stable distribution
Stable(a::Real, b::Real) = Stable(a, b, 1.0, 0.0)

# Algorithm from Mallows
function Distributions.rand(d::Stable)
	const piby2 = π / 2

	α, β, γ, δ = d.alpha, d.bprime, d.gamma, d.delta

	u = rand()
	w = rand(Exponential())

	ε = 1 - α

	# Compute some tangents
	phiby2 = piby2 * (u - 0.5)
	a = phiby2 * tan2(phiby2)
	bb = tan2(ε * phiby2)
	b = ε * phiby2 * bb
	if ε > -0.99
		τ = β / (tan2(ε * piby2) * piby2)
	else
		τ = β * piby2 * ε * (1 - ε) * tan2((1 - ε) * piby2)
	end

	# Compute some necessary subexpressions
	da = a^2
	db = b^2
	a2 = 1.0 - da
	a2p = 1.0 + da
	b2 = 1.0 - db
	b2p = 1.0 + db

	# compute coefficient
	z = a2p * (b2 + 2 * phiby2 * bb * τ) / (w * a2 * b2p)

	# compute the exponential-type expression
	logz = log(z)
	d = d2(ε * logz / (1 - ε)) * (logz / (1 - ε))

	# compute "standard" stable with std = sqrt(2)
	s = (1 + ε * d) * 2 *
	      ((a - b) * (1 + a * b) - phiby2 * τ * bb * (b * a2 - 2 * a))
	s /= (a2 * b2p)
	s += τ * d

	# shift scale and location
	s = s / sqrt(2) * γ + δ

	return s
end