Imaginary numbers come straight from math class in school. Imaginary numbers are usually taught first in a highschool Algebra II class, but anyone who has had any algebra at all should be able to understand this document. So here goes.
     
Lets start with the basics. When you square a number, you multiply that number by itself. The square of 2 is 4 then, since 2*2=4. When you take the square root of a number, you do the reverse. If the square of 2 is 4, the square root of 4 is 2. Ok, so you knew that already.
     
When you multiply a number by itself, the sign of the result is always positive. This is because whenever two numbers with like signs (that is two positive numbers or two negative numbers) are multiplied, the result is always positive.
     
The question I'm leading to then is what happens if you need to take the square root of a negative number. We call the square of a negative number an imaginary number. An imaginary number is denoted by an i.
sqrt(-1) = i
sqrt(-4) = sqrt(-1) * sqrt(4) = i * 2 = 2i
     
Imaginary numbers cannot be evaulated, as they have i as a product, and i is undefined. This is why they are called imaginary numbers. They don't really exist on our number line. Mathematicians, especially those who work in the fields of chaos and fractals, are familiar enough with imaginary numbers to no longer consider them imaginary just because their value is undefined. While they can't be represented on our number line, they can be represented on a two-dimensional number plane: the complex plane.
     
The complex plane is very similar to the Cartesian plane used to graph y in terms of x for a function. The Cartesian plane uses one coordinate for input (usually x), and one coordinate for output (usually y). The complex plane, however, uses both coordinates just to describe one input number, and outputs a color or height value at each coordinate to generate a fractal image or 3d graph. Let's take a look.

     
To code the Julia and Mandelbrot Sets, there are only three operations that we need to perform with complex numbers: add, square and find the absolute value. Complex numbers add just like polynomials with two terms. They square like polynomials (a^2 + b^2 + 2ab), but you wind up subtracting the square of b from a, since an imaginary number squared is a negative.
Addition:
  (a + bi) + (c + di) = (a + c) + (b + d)i
Computing the square:
  (a + bi)^2 = (a^2 - b^2) + 2abi
The absolute value of a number is defined as the distance from zero. This can be evaluated using the Pythagorean theorem. If you drew an imaginary vertical line between the point you want to evaluate to the x-axis and another from the point to the origin, you can use the triangle that is formed to determine the distance of the point from the origin, or from zero.

     
You may be a bit confused as to how to get your computer to recognize a number as imaginary. Certain Fortran versions can internally recognize a complex number data type. If you're using a complex number friendly language, then this part might not apply to you.
     
The trick to programming complex math is to use two variables for each number. For each number, have a real part variable, and an imaginary part variable. So for 3 - 2i, 3 would be stored in the real part, and 2 in the imaginary part. In Pascal or C, you can make a complex number structure, with a real and imaginary float. In C++, you can make a complex number class, and if you code it right, you could have the class operate using the standard math operators, as if there was an internal complex number data type. In assembly, you're on your own. You can do structures with most assemblers, and TASM even has some rudimentary pseudo object oriented features, or so I'm told.
     
Here are the basics of the complex number math that we'll need again, this time in pseudocode.
Variables:
  xr = real part of complex number x
  xi = imaginary part of complex number x
  yr = real part of complex number y
  yi = imaginary part of complex number y
  z  = floating point number (not complex)
Addition: x = x + y
  xr = xr + yr
  xi = xi + yi
Computing the square: x = y^2
  xr = (yr * yr) - (yi * yi)
  xi = 2 * yr * yi
Finding the absolute value: z = abs(x)
  z = sqrt(xr * xr + xi * xi)
     
To graph a function on a Cartesian plane, you walk the x-axis graphing the result of a funciton performed on the x-position to the y-axis. In a fractal, you walk both axes, and graph the result of a function as a color to each pixel. A fractal is just a function done for every point on a complex plane, with each resulting number interpreted as a color.
     
Most fractals use the number of iterations of an operation needed for a certain condition to be met as the color number or color 0 if that condition is never met. This is how Mandelbrot and Julia fractals work.
     
The Julia Set is graphed by calculating the complex coordinate of each pixel (Z), and by squaring that number and adding a second constant complex number (C) each iteration. If that number never goes to infinity, then that number is in the set, and the pixel is colored with color 0 (traditionally black). If not, the color corresponding to the number of iterations it took to determine that the number is going to infinity is plotted.
     
If a complex number's absolute value is greater than or equal to two, it can be assumed that the number will go to infinity if we continued squaring it and adding to it. I do not know if this can be proved. If you know, then please post an explanation or a link on the messageboard and I'll put it on the site.
     
Here's the pseudocode:
c = magic number
repeat (num_pixels) times {
  z = calculate coordinate of current pixel
  iterations = 0
  repeat (maximum_iterations) times {
    iterations = iterations + 1
    z = (z * z) + c
    if (abs(z) >= 2) then goto nextpixel
  }
  nextpixel:
  if (iterations = maximum_iterations) then iterations = 0
  colorpixel(iterations)
}
     
There are an infinite number of Julia Sets, as the constant C can be any complex number. There are also an infinite number of optimizations that can be done to that pseudocode, but that was written for readibility.
     
The Mandelbrot Set, modeled after the Julia Set, uses the same basic process for each pixel, starting with Z equal to zero instead of starting out with the pixel's complex number value, and adding the pixel's number value (C) each iteration after squaring what you've got (Z) intead of adding a constant. Here's some more pseudocode.
repeat (num_pixels) times {
  c = calculate coordinate of current pixel
  z = 0
  iterations = 0
  repeat (maximum_iterations) times {
    iterations = iterations + 1
    z = (z * z) + c
    if (abs(z) >= 2) then goto nextpixel
  }
  nextpixel:
  if (iterations = maximum_iterations) then iterations = 0
  colorpixel(iterations)
}
     
Just as with the Julia Set, you can create different fractals by using different C values, in the Mandelbrot Set, you can use different starting Z values. You can also you compute Z to a power other than 2 before adding C. Z^3+C is a cubic Mandelbrot, etc.
     
Fractals are fairly simple, and with a decent palette they can be very impressive. I've made a Mandelbrot Set grapher with several palettes which is available with source (dos assembly) and I've made a much simpler Mandelbrot Set plotter in BASIC for the TI-82 graphic calculator (roughly 3 hours to render each frame) which will be available once I ask Andrew Logan to transfer it to computer w/ his link and email it to me.