END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

Start by learning the math. Check the mandelbulb wiki. Instead of complex numbers, triplex numbers are used and
Z=Z^2+C becomes:(in glsl)

float r=length(Z);   //Z is now a 3d vector
float zo = asin(Z.z / r) * p;  //p is power, usually 8 instead of 2
float zi = atan(Z.y, Z.x) * p;
Z=pow(r,p)*vec3(cos(zo)*cos(zi),cos(zo)*sin(zi),sin(zo))+C; //there are variations to this swapping sin/cos

But now you have a bigger problem - how to render a 3d object to a 2d screen.
Learn about ray-marching with distance estimation. A distance estimation formula tells you how far you can safely move the ray until you need to check again. A mandelbulb DE formula looks like this:
float DE(vec3 Z){
  vec3 C=Z;
  float dr = 1.0; //this is a running derivative
  float r = length(Z);
  for (int n = 0; n < iterations && r<2.0; n++) {//about 7 iterations will do for starters
    float zo = asin(Z.z / r) * p;
    float zi = atan(Z.y, Z.x) * p;
    float zr = pow(r, p-1.0);
    dr = dr * zr * p + 1.0; //do you see what is meant by a running derivative here?
    Z=zr*r*vec3(cos(zo)*cos(zi),cos(zo)*sin(zi),sin(zo))+C;
    r = length(Z);
  }
  return 0.5*r*log(r)/dr;
}
TL;DR: Check out Boxplorer code for an OpenGL sample.
 Oh wait a MandelBOT!?! Sorry can't help

Very good. When the return value dist=DE(ray) is less than the size of a pixel at that distance than you quit and say I hit the fractal. Like this
rayDir is the normalized direction of the ray march
ray is the point we are testing, it starts at the camera lens
eps is an estimate of the pixel size
maxRayLen is the bound and marching steps can vary, about 100 to start
for(int i=0;i<marchingSteps && dist>eps && rayLength<maxRayLen;i++){
  rayLen+=dist=DE(ray+rayLen*rayDir);
}
This is actually quite fast! Think about a simple DE like a sphere of radius 1...
float DE(Z){return length(Z)-1.0;}
If you are aiming right at it you hit in 1 step. If you miss it completely it takes a few more and you reach maxRayLen.
The interesting part is near the edges. This will take many steps and may seem inefficient but the number of steps gives you an approximate occlusion of light so you could even just color the fractal with...
color=diffuse*(1.0-float(i)/float(marchingSteps));

But thats your next step learning lighting calculations.

Yes there are programs that use monte carlo methods etc. If you are using the GPU each pixel is oblivious of its neighbors.
The method I've described runs near real time depending on your hardware. Many thanks to the folks here who have contributed to it!