Calculating surface normals for 3D fractals

[Derakon]

I'm working on making my own raytracer for 3D fractals. I have a very basic system made (in about 120 lines of Python), that I think is rendering the eighth-order Mandelbulb...but it's a bit hard to tell because I don't have proper surface normals!



That's with the camera at <2, 2, 2>, looking at the origin. The shading that you can see is just varying with the distance of the surface from the camera.

What's the recommended way to generate surface normals for 3D fractals? The only method I can think of is to sample points around the ray intersection point to get the local slope, but that seems pretty hackish, not to mention inefficient. Any better ideas?