I like the way you are constructing the fractal directly - it is easier to understand - but typically since the code becomes simpler we would write the map loop like this:
float fd=(length(p)-1.0)/bs; //divide the result by how much we have scaled up p
p=abs(p);
if(p.x<p.y)p.xy=p.yx;//sorting the longest dim into x
if(p.y<p.z)p.yz=p.zy;
if(p.x<p.y)p.xy=p.yx;
p.xy-=1.0; //move xy back towards origin 0.92 makes them just touch
p*=3.3333; //scale up - makes the result smaller
bs*=3.3333; //bs is now keeping track of the running scale
if (fd < d) {
d = fd;
r = float(i)/3.0;
}
edit: corrected the offset - my bad