This article is from the sci.fractals FAQ, by Michael C. Taylor and Jean-Pierre Louvet with numerous contributions by others.
The Julia set can be computed by iteration similar to the
Mandelbrot computation. The only difference is that the c value is
fixed and the initial z value varies.
Alternatively, points on the boundary of the Julia set can be computed
quickly by using inverse iterations. This technique is particularly
useful when the Julia set is a Cantor Set. In inverse iteration, the
equation z1 = z0^2 + c is reversed to give an equation for z0: z0 =
±sqrt(z1 - c). By applying this equation repeatedly, the resulting
points quickly converge to the Julia set boundary. (At each step,
either the positive or negative root is randomly selected.) This is a
nonlinear iterated function system.
z = 1 (or any value)
if (random number < .5) then
z = sqrt(z - c)
z = -sqrt(z - c)