Pictures of Julia and Mandelbrot Sets/Computer programs
We will show two ways of making a program that can draw a Mandelbrot set where the iterations are towards ∞: a program where all is done from the ground, but which does nothing more than drawing a single picture, and a program where you can do all you desire (zoom, alter colouring and render as file), but where these facilities are handed over to another program, namely Ultra Fractal.
A single picture[edit]
Let us make a program that is as simple as possible: it draws a single picture, and does nothing more (when it has done its work it closes). It was originally written in assembly language, but here we write it in pseudocode. It consists of two parts: the drawing procedure and the colour scale.
We draw a section of the Mandelbrot set for the family , for p = 0.009. We let the size of the window be 800x600 pixels, and we imagine that the section has its lower left corner in the point with coordinates (ax, ay) and that is has width h. We have , so that the finite critical point we need, is solution to the equation . We choose the real solution cri = .
The drawing procedure
 p = 0.009
 cri =
 r = 1.0E200 (square of the bailout radius)
 u = log(log(r))
 v = 1/log(2) (2 is the degree of the rational function)
 ax = 0.7028233689263 (xcoordinate of lower left corner)
 ay = 0.1142331238418 (ycoordinate of lower left corner)
 h = 0.0092906805859 (width of the section)
 g = h / 800
 m = 850 (maximum iteration number)
 thick = h * 0.0005 (thickness of the boundary)
 dens = 24.9 (density of the colours)
 disp = 432.0 (displacement of the colour scale)
for i = 0 to 799 do
 begin
 cx = ax + i * g (xcoordinate of pixel (i, j))
 for j = 0 to 599 do
 begin
 cy = ay + j * g (ycoordinate of pixel (i, j))
 x = cri
 y = 0
 xd = 0
 yd = 0
 f = 0
 n = 0
 while (n < m) and (f < r) do
 begin
 n = n + 1
 x1 = x * x  y * y
 y1 = 2 * x * y
 f = x * x + y * y
 x2 = 2 * x  p * x1 / (f * f)
 y2 = 2 * y + p * y1 / (f * f)
 temp = xd
 xd = x2 * xd  y2 * yd + 1
 yd = x2 * yd + y2 * temp
 fd = xd * xd + yd * yd
 x = x1 + p * x / f + cx
 y = y1  p * y / f + cy
 f = x * x + y * y
 end
 begin
 if (n = m) or (log(f) * sqrt(f) < thick * sqrt(fd)) then
 setpixel(i, 599  j, 0)
 else
 begin
 s = n  v * (log(log(f))  u)
 n = round(dens * s + disp) mod 720
 col = paletteRGB(col1[n], col2[n], col3[n])
 setpixel(i, 599  j, col)
 end
 begin
 end
 begin
 end
Explanation
We have set , , and , and we have used that .
The successive calculation of the derivative z' is , where and . The next point in the iteration is . The distance function is
 log(√f)*√f/√fd (= log(f)*√f/(2*√fd))
for the last zvalue, here and . The number in the interval [0, 1[ to be subtracted from the iteration number (in order to get the real iteration number), is , for the last z, here and .
paletteRGB(r, g, b) is the integer indexing the colour with RGBvalues (r, g, b), 0 gives black. col, col2 and col3 are the arrays from 0 to 719 of integers from 0 to 255 constructed in the next section.
The colour scale
A colour is immediately (in the computer) given by its composition of the three primary colours red, green and blue, and their shares are measured in whole numbers between 0 and 255. This triple of numbers is the set of RGB values of the colour. The colours correspond accordingly to the entire points in a cube with side length 256, and we construct our cyclic colour scala by going regularly along an ellipse in this cube. An ellipse in the space is determined by:

 a centre with coordinates (a, b, c)
 a major axe rmaj and a minor axe rmin
 two angles angle g and h determining the direction of the plane of the ellipse
 an angle q determining the direction of the ellipse in this plane
 u = cos(g * pi / 180) * cos(h * pi / 180) {(u,v,w) is a unit vector orthogonal to the plane}
 v = cos(g * pi / 180) * sin(h * pi / 180)
 w = sin(g * pi / 180)
 x = a {(x,y,z) is a vector in the plane}
 y = b
 z = a * u / w + b * v / w
 e = sqrt(sqr(x) + sqr(y) + sqr(z))
 x1 = x / e {the unit vector corresponding to (x,y,z)}
 y1 = y / e
 z1 = z / e
 x2 = v * z1  w * y1 {the unit vector in the plane orthogonal to (x1,y1,z1)}
 y2 = w * x1  u * z1
 z2 = u * y1  v * x1
 e = cos(q * pi / 180)
 f = sin(q * pi / 180)
 x1 = e * x1 + f * x2 {the unit vector in the plane in the direction of the angle q}
 y1 = e * y1 + f * y2
 z1 = e * z1 + f * z2
 x2 = v * z1  w * y1 {the unit vector in the plane orthogonal to this direction}
 y2 = w * x1  u * z1
 z2 = u * y1  v * x1
 for i = 0 to 719 do
 begin
 e = rmaj * cos(i * pi / 360)
 f = rmin * sin(i * pi / 360)
 col1[i] = round(a + e * x1 + f * x2) {the three coordinates of the point on the ellipse}
 col2[i] = round(b + e * y1 + f * y2)
 col3[i] = round(c + e * z1 + f * z2)
 end
 begin
In the picture we have: (a, b, c) = (176, 176, 176), rmaj = rmin = 88, g = 146, h = 32, q = 0.
Ultra Fractal[edit]
Ultra Fractal is a program for drawing fractals where the user to a great extent can write his own formula programs and where the main program performs all that is common for all the fractal procedures, that is, the procedure of going from pixel to pixel, zooming, colouring, production of a large picture as a file. Your effort is reduced to a minimum, and apart from pictures that cannot be drawn regularly from pixel to pixel  attractors and landscapes, for instance  you can do all you desire: noncomplex functions, quaternions, ... You can operate directly with complex numbers.
First you should download a free trialversion of the program and learn how it works. Then you should see the article Make attractive pictures of Mandelbrot and Julia sets with Ultra Fractal, and download the (free) programs from this site to be run with Ultra Fractal.
For the drawing two programs have to be combined: a formula program and a colouring program. In the formula program is a section called loop, and when this loop has done its work, the number of iterations and the last value of the complex number z are transferred to the colouring program, which calculates the colour from these two numbers. However, Ultra fractal is designed to "all" types of fractals, and unfortunately the Julia and Mandelbrot sets do not quite fit into this form: for iterations towards cycles, the iteration has to be continued in order to find the order and the attraction of the cycle, and it is not the last value of z in the sequence we do use for the colour. Therefore we replace this loop by a fictive loop (which only runs one time) and perform all the operations in the section init. This implies that we cannot use the default colouring program of Ultra Fractal, we must write a new colouring program.
Another thing you must be aware of, is that in Ultra Fractal the norm z of the complex number z is the square of its length: z = .
We will show a program that draws the Mandelbrot set for the family , where is a polynomial whose first two terms are zero, so that it begins with or a larger power of z, because we then can take 0 as the finite critical point.
The two programs can be copied and inserted in an empty formula and colouring document, respectively. We have inserted the polynomial of degree 4, where p is a real parameter. The picture shows a section of the Mandelbrot set for p = 0.26.
The formula program
 Mandelbrot {
 global:
 float p = 0.26 ;parameter
 float deg = 4 ;degree of the polynomial
 float v = 1 / log(deg)
 float g = 10 * log(10)
 float r = exp(g) ;square of the radius of the bailout circle
 u = log(g)
 float tb = sqr(@thick / (1000 * #magn)) ;for the thickness of the boundary
 float h = 1 / (1500 * #magn * @width) ;a very small real number
 init:
 complex z = 0 ;critical point
 complex zd = 0 ;the sequence of the derivatives
 complex z1 = 0
 float w = 0
 int n = 0
 while n < #maxit && z < r
 n = n + 1
 z1 = z^2/2 + p*z^4
 zd = ((z+h)^2/2 + p*(z+h)^4  z1) * zd / h + 1
 z = z1 + #pixel
 endwhile
 if n == #maxit  sqr(log(z)) * z < tb * zd
 w = 1
 else
 w = n  v * (log(log(z))  u)
 endif
 ;begin fictive loop
 z = 0
 n = 0
 loop:
 n = n + 1
 z = z + #pixel
 if n == 1
 z = w
 endif
 bailout:
 n < 1
 ;end fictive loop
 default:
 title = "Mandelbrot"
 maxiter = 100
 param thick
 caption = "boundary"
 default = 1.0
 endparam
 param width
 caption = "width"
 default = 640
 endparam
 }
Explanation
We have set the square of the radius r of the bailout circle to , and set the thickness of the boundary to "thick/(1000 * #magn)", so that it becomes thinner when we zoom in.
We have calculated the derivative by for a small number h, namely "h = 1/(1500*#magn*width)", where "width" is the width of the picture. The default value is the width of the window in pixels (here set to 640), but if you draw a large picture, you should enter the width of this.
The successive calculation of the derivative is "zd = (f(z+h)  f(z)) * zd / h + 1". The next iteration is "z = f(z) + #pixel". The square of the distance estimation is "sqr(log(z)) * z / zd". The number to be subtracted from the iteration number is "v * (log(log(z))  u)", where "v = 1/log(deg)" and "u = log(log(r))".
The final complex number z (to be transferred to the colouring program), which here actually is real, is set to 1, when the point belongs to the Mandelbrot set or to the boundary, otherwise it is set to the real iteration number. It is transferred to the colouring program, which we have called Gradient, as #z, and is here set equal to the real number s. When s is negative the colour is set to "#solid", otherwise we multiply s by a number "dens" determining the density, and add to this number a number "disp" determining the displacement of the scale, and as we want this displacement to be a per cent, we divide the result by 100: "u = (dens * s + disp) / 100". In Ultra Fractal the colours of the cyclic colour scales are indexed by the real numbers in the interval [0, 1[, and we let this number, "#index", be the nonintegral part of the number u: "#index = u  trunc(u)".
The colouring program
 Gradient {
 final:
 float s = real(#z)
 float u = 0
 if s < 0
 #solid = true
 else
 u = (@dens * s + @disp) / 100
 #index = u  trunc(u)
 endif
 default:
 title = "Gradient"
 param disp
 caption = "displace"
 default = 0
 endparam
 param dens
 caption = "density"
 default = 1.0
 endparam
 }