Intro
Producing simple analytic shapes with for a raymarching engine is very simple, as seen in the in this reference article in
this very webstite. Besides simple combination of shapes to construct more complex compound shapes, there is also the possibility to do this composition algorithmically. Probably
the most basic way of algorithmic composition, is the recursive introduction of regular smaller details. This naturally produces classic Cantor fractals. A good example of
that is the "untraceable" 1 kilobyte demo by TBC.
or even better, vec3 map( in vec3 p ) { float d = sdBox(p,vec3(1.0)); float s = 1.0; for( int m=0; m<3; m++ ) { vec3 a = mod( p*s, 2.0 )1.0; s *= 3.0; vec3 r = abs(1.0  3.0*abs(a)); float da = max(r.x,r.y); float db = max(r.y,r.z); float dc = max(r.z,r.x); float c = (min(da,min(db,dc))1.0)/s; d = max(d,c); } return vec3(d,1.0,1.0); } Last step is, to compute the material id based on the iteration count, and some fake occlusion too (see images on the right of the article to see how it looks like when all these elements are put together): vec3 map( in vec3 p ) { float d = sdBox(p,vec3(1.0)); vec3 res = vec3( d, 1.0, 0.0, 0.0 ); float s = 1.0; for( int m=0; m<3; m++ ) { vec3 a = mod( p*s, 2.0 )1.0; s *= 3.0; vec3 r = abs(1.0  3.0*abs(a)); float da = max(r.x,r.y); float db = max(r.y,r.z); float dc = max(r.z,r.x); float c = (min(da,min(db,dc))1.0)/s; if( c>d ) { d = c; res = vec3( d, 0.2*da*db*dc, (1.0+float(m))/4.0, 0.0 ); } } return res; } 

Some neat tricks can be applied during the iterative substraction of cubes, suh as translating and rotating the point p a little bit in each iteration. That produces less symmetrical patterns, as can be seen in the images to the right of this article or here below (which was rendered with a simple pathtracer):
And finally, this is the source code for a reference implementation, rendered in realtime (click the play button to see it mobe, and the title in the image to jump to the source code):