Welcome to Fractal Forums

//scale=2
//bailout=1000
sierpinski3(x,y,z){
   r=x*x+y*y+z*z;
   for(i=0;i<10 && r<bailout;i++){
      //Folding... These are some of the symmetry planes of the tetrahedron
      if(x+y<0){x1=-y;y=-x;x=x1;}
      if(x+z<0){x1=-z;z=-x;x=x1;}
      if(y+z<0){y1=-z;z=-y;y=y1;}
      
      //Stretche about the point [1,1,1]*(scale-1)/scale; The "(scale-1)/scale" is here in order to keep the size of the fractal constant wrt scale
      x=scale*x-(scale-1);//equivalent to: x=scale*(x-cx); where cx=(scale-1)/scale;
      y=scale*y-(scale-1);
      z=scale*z-(scale-1);
      r=x*x+y*y+z*z;
   }
   return (sqrt(r)-2)*scale^(-i);//the estimated distance
}

//scale=2
//bailout=1000
sierpinski3(x,y,z){
   r=x*x+y*y+z*z;
   for(i=0;i<10 && r<bailout;i++){
      rotate1(x,y,z);
      
      if(x+y<0){x1=-y;y=-x;x=x1;}
      if(x+z<0){x1=-z;z=-x;x=x1;}
      if(y+z<0){y1=-z;z=-y;y=y1;}
      
      rotate2(x,y,z);

      x=scale*x-(scale-1);
      y=scale*y-(scale-1);
      z=scale*z-(scale-1);
      r=x*x+y*y+z*z;
   }
   return (sqrt(r)-2)*scale^(-i);//the estimated distance
}

//scale=2
//bailout=1000
sierpinski3(x,y,z){
   r=x*x+y*y+z*z;
   for(i=0;i<10 && r<bailout;i++){
      rotate1(x,y,z);
      
      if(x+y<0){x1=-y;y=-x;x=x1;}
      if(x+z<0){x1=-z;z=-x;x=x1;}
      if(y+z<0){y1=-z;z=-y;y=y1;}
      
      rotate2(x,y,z);

      x=scale*x-CX*(scale-1);
      y=scale*y-CY*(scale-1);
      z=scale*z-CZ*(scale-1);
      r=x*x+y*y+z*z;
   }
   return (sqrt(r)-2)*scale^(-i);//the estimated distance
}

      if(x+y<0){x1=-y;y=-x;x=x1;}
      if(x+z<0){x1=-z;z=-x;x=x1;}
      if(y+z<0){y1=-z;z=-y;y=y1;}

      if(x-y<0){x1=y;y=x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}

      if(x-y<0){x1=y;y=x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}
      if(x+y<0){x1=-y;y=-x;x=x1;}
      if(x+z<0){x1=-z;z=-x;x=x1;}
      if(y+z<0){y1=-z;z=-y;y=y1;}

      x=abs(x);y=abs(y);z=abs(z);

      if(x-y<0){x1=y;y=x;x=x1;}
      if(x+y<0){x1=-y;y=-x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(x+z<0){x1=-z;z=-x;x=x1;}

      x=abs(x);y=abs(y);z=abs(z);
      if(x-y<0){x1=y;y=x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}

Menger3(x,y,z){
   r=x*x+y*y+z*z;
   for(i=0;i<MI && r<bailout;i++){
      rotate1(x,y,z);

      x=abs(x);y=abs(y);z=abs(z);
      if(x-y<0){x1=y;y=x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}

      rotate2(x,y,z);
    
      x=scale*x-CX*(scale-1);
      y=scale*y-CY*(scale-1);
      z=scale*z;
      if(z>0.5*CZ*(scale-1)) z-=CZ*(scale-1);
      
      r=x*x+y*y+z*z;
   }
   return (sqrt(x*x+y*y+z*z)-2)*scale^(-i);
}

Menger3(x,y,z){
   r=x*x+y*y+z*z;
   for(i=0;i<MI && r<bailout;i++){
      rotate1(x,y,z);

      x=abs(x);y=abs(y);z=abs(z);
      if(x-y<0){x1=y;y=x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}
      
      z-=0.5*CZ*(scale-1)/scale;
      z=-abs(-z);
      z+=0.5*CZ*(scale-1)/scale;

      rotate2(x,y,z);
    
      x=scale*x-CX*(scale-1);
      y=scale*y-CY*(scale-1);
      z=scale*z;
      
      r=x*x+y*y+z*z;
   }
   return sqrt(x*x+y*y+z*z)*scale^(-i);
}

sierpinski(x,y,z){
   scale=2;
   r=0;
   for(i=0;i<MI && r<1;i++){
      x1=rot[0][0]*x+rot[1][0]*y;
      y1=rot[0][1]*x+rot[1][1]*y;
      z1=z;
      x=x1;y=y1;z=z1;
      
      if(x-y<0){x1=y;y=x;x=x1;}
      if(x+y<0){x1=-y;y=-x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(x+z<0){x1=-z;z=-x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}
      if(y+z<0){y1=-z;z=-y;y=y1;}
      
      r=(x+y-z-1)*1/sqrt(3);
      
      x=scale*x-(scale-1);
      y=scale*y-(scale-1);
      z=scale*z-(scale-1);
   }
   return r*scale^(1-i);
}

menger(x,y,z){
   scale=3;
   r=0;
   for(i=0;i<MI && r<1;i++){
      x=abs(x);y=abs(y);z=abs(z);
      if(x-y<0){x1=y;y=x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}
      
      x1=x-1;y1=y-1;z1=z-1;
      r=max(x1,max(y1,z1));
      /*nice alternative that gives rounded (external) edges.
      It's not necessary to calc the sqrt() at eache iteration
      if(x1<0 && y1<0 && z1<0) r=max(x1,max(y1,z1));
      else r=sqrt(max(0,x1)^2+max(0,y1)^2+max(0,z1)^2);
      */
      
      x=scale*x-(scale-1);
      y=scale*y-(scale-1);
      z=scale*z;
      if(z>0.5*(scale-1)) z-=(scale-1);
   }
   r*scale^(1-i)
}

Mengercolor(x,y,z){
   r=x*x+y*y+z*z;
   c=0;
   for(i=0;i<MI;i++){//It is better to use a fixed iterations number to avoid shift when zooming.

      rotate1();

      c1=0;
      c1+=0.5*(1-sgn(x));c1+=0.5*(1-sgn(y));c1+=0.5*(1-sgn(z));
      x=abs(x);y=abs(y);z=abs(z);
      if(x-y<0){x1=y;y=x;x=x1;c1+=1;}
      if(x-z<0){x1=z;z=x;x=x1;c1+=1;}
      if(y-z<0){y1=z;z=y;y=y1;c1+=1;}

      rotate2();  
  
      x=scale*x-stc[0]*(scale-1);
      y=scale*y-stc[1]*(scale-1);
      z=scale*z;
      if(z>0.5*stc[2]*(scale-1)) {z-=stc[2]*(scale-1);c1+=1;}
      c+=c1*0.5;
      r=x*x+y*y+z*z;
   }
   return c;//this will be used as an index into a palette
}

 scale=2
       if x+y<0, x1=-y,y=-x,x=x1,endif
       if x+z<0, x1=-z,z=-x,x=x1,endif
       if x-y<0, x1=y,y=x,x=x1,endif
       if x-z<0, x1=z,z=x,x=x1,endif

      x=x*scale-1*(scale-1)
      y=y*scale-0*(scale-1)
      z=z*scale-0*(scale-1)

#define _IVNORM_ (0.5/_PHI_)
#define _PHI1_ (_PHI_*_IVNORM_)
#define _1PHI_ (_IVNORM_)
#define _PHI2_ (_PHI_*_PHI_*_IVNORM_)

#define _IKVNORM_ 1/sqrt((_PHI_*(1+_PHI_))^2+(_PHI_^2-1)^2+(1+_PHI_)^2)
#define _C1_ (_PHI_*(1+_PHI_)*_IKVNORM_)
#define _C2_ ((_PHI_*_PHI_-1)*_IKVNORM_)
#define _1C_ ((1+_PHI_)*_IKVNORM_)
dodecasierpinski3(x,y,z,scale,rot[3][3]){
   r=x*x+y*y+z*z;
   i=0;
   for(i=0;i<MaxIteration && r<Bailout; i++){
#ifdef PRE_ROTATE
      x1=rot[0][0]*x+rot[1][0]*y+rot[2][0]*z;
      y1=rot[0][1]*x+rot[1][1]*y+rot[2][1]*z;
      z1=rot[0][2]*x+rot[1][2]*y+rot[2][2]*z;
      x=x1;y=y1;z=z1;
#endif      
      x=abs(x);
      y=abs(y);
      z=abs(z);
      t=_PHI2_*x+_1PHI_*y-_PHI1_*z;
      if(t<0){x=x-2*t*_PHI2_;y=y-2*t*_1PHI_;z=z+2*t*_PHI1_;}

      t=-_PHI1_*x+_PHI2_*y+_1PHI_*z;
      if(t<0){x=x+2*t*_PHI1_;y=y-2*t*_PHI2_;z=z-2*t*_1PHI_;}

      t=_1PHI_*x-_PHI1_*y+_PHI2_*z;
      if(t<0){x=x-2*t*_1PHI_;y=y+2*t*_PHI1_;z=z-2*t*_PHI2_;}

      t=-_C1_*x+_C2_*y+_1C_*z;
      if(t<0){x=x+2*t*_C1_;y=y-2*t*_C2_;z=z-2*t*_1C_;}

      t=_1C_*x-_C1_*y+_C2_*z;
      if(t<0){x=x-2*t*_1C_;y=y+2*t*_C1_;z=z-2*t*_C2_;}
      
#ifdef POST_ROTATE
      x1=rot[0][0]*x+rot[1][0]*y+rot[2][0]*z;
      y1=rot[0][1]*x+rot[1][1]*y+rot[2][1]*z;
      z1=rot[0][2]*x+rot[1][2]*y+rot[2][2]*z;
      x=x1;y=y1;z=z1;
#endif      
      x=scale*x-(scale-1)*stc[0];//stc is the center of scaling
      y=scale*y-(scale-1)*stc[1];
      z=scale*z-(scale-1)*stc[2];
      r=x*x+y*y+z*z;
   }
   return (sqrt(r)-2)*scale^(-i);
}

procedure Rotate2(const Rx,Ry,Rz:Double; const x,y,z:Double; Var Nx,Ny,Nz:Double);
Var TempX  : Double;
    TempY  : Double;
    TempZ  : Double;
    SinX   : Double;
    SinY   : Double;
    SinZ   : Double;
    CosX   : Double;
    CosY   : Double;
    CosZ   : Double;
    XRadAng: Double;
    YRadAng: Double;
    ZRadAng: Double;
Begin

 XRadAng := Rx*PIDiv180;
 YRadAng := Ry*PIDiv180;
 ZRadAng := Rz*PIDiv180;

 SinX    := Sin(XRadAng);
 SinY    := Sin(YRadAng);
 SinZ    := Sin(ZRadAng);

 CosX    := Cos(XRadAng);
 CosY    := Cos(YRadAng);
 CosZ    := Cos(ZRadAng);

 Tempy   := y*CosY - z*SinY;
 Tempz   := y*SinY + z*CosY;
 Tempx   := x*CosX - Tempz*SinX;

 Nz      := x*SinX     + Tempz*CosX;
 Nx      := Tempx*CosZ - TempY*SinZ;
 Ny      := Tempx*SinZ + TempY*CosZ;

End;

procedure Rotate(const Rx,Ry,Rz:Double; const x,y,z,ox,oy,oz:Double; Var Nx,Ny,Nz:Double);
Begin
 Rotate2(Rx,Ry,Rz,x-ox,y-oy,z-oz,Nx,Ny,Nz);
 Nx := Nx+ox;
 Ny := Ny+oy;
 Nz := Nz+oz;
End;

// Return rotation matrix for rotating around vector v by angle
float3x3 rotationMatrix(float3 v, float angle)
{
	float c = cos(radians(angle));
	float s = sin(radians(angle));

	return float3x3(c + (1.0 - c) * v.x * v.x, (1.0 - c) * v.x * v.y - s * v.z, (1.0 - c) * v.x * v.z + s * v.y,
			    (1.0 - c) * v.x * v.y + s * v.z, c + (1.0 - c) * v.y * v.y, (1.0 - c) * v.y * v.z - s * v.x,
			    (1.0 - c) * v.x * v.z - s * v.y, (1.0 - c) * v.y * v.z + s * v.x, c + (1.0 - c) * v.z * v.z);
}

--

# 3x3 rotation matrix
float3x3 objRotation = rotationMatrix(float3(1, 0, 0), objectRotation.x) * rotationMatrix(float3(0, 1, 0), objectRotation.y) * rotationMatrix(float3(0, 0, 1), objectRotation.z);

      //salpha: rotation about Z axis angle
      //sbeta: rotation about Y axis angle
      //stheta: rotation about X axis angle
      s1=sin(sbeta);c1=cos(sbeta);
      s2=sin(stheta);c2=cos(stheta);
      s3=sin(salpha);c3=cos(salpha);
      rotmat[0][0]=c1*c3+s1*s2*s3; rotmat[0][1]=c2*s3; rotmat[0][2]=c1*s2*s3-c3*s1;
      rotmat[1][0]=c3*s1*s2-c1*s3; rotmat[1][1]=c2*c3; rotmat[1][2]=s1*s3+c1*c3*s2;
      rotmat[2][0]=c2*s1;          rotmat[2][1]=-s2;   rotmat[2][2]=c1*c2;

      x1=rot[0][0]*x+rot[1][0]*y+rot[2][0]*z;
      y1=rot[0][1]*x+rot[1][1]*y+rot[2][1]*z;
      z1=rot[0][2]*x+rot[1][2]*y+rot[2][2]*z;
      x=x1;y=y1;z=z1;

for each I
 float length = dot(point, Direction[I])
 if (length > Distance[I])
   point -= 2*Direction[I]*(length-Distance[I])  // fold
end
point.rotate(AngleVector)
point += Offset
point *= Scale

point.rotate(AngleVector[k1])
for each I
  point.rotate(AngleVector[i1])
  float length = dot(point, Direction[I])
  if (length > Distance[I])
    point -= 2*Direction[I]*(length-Distance[I])  // fold
  point.rotate(AngleVector[i2])
end
point.rotate(AngleVector[k2])
point += Offset
point *= Scale

point.rotate(AngleVector0)
for each I
  float length = dot(point, Direction[I])
  if (length > Distance[I])
    point -= 2*Direction[I]*(length-Distance[I])  // fold
  point.rotate(AngleVector[i])
end
point += Offset
point *= Scale

point = mRotation.RotateVector(point);
for(int i=0; i<IFSfoldingCount; i++)
{
   length = point.Dot(IFSDirection[i]);
   if(length > IFSDistance[i])
   {
	point -= IFSDirection[i]*2.0*(length - IFSDistance[i]);
   }
}
point += IFSOffset;
z *= scale;

Vector dir = randomDirectionVector()
float dot = point.dot(dir)
float dist = randomFloat0to1()
if (dot > dist)
 point -= 2*dir*(dot-dist)
elseif (dot < -dist)
 point -= 2*dir*(dot+dist)
point *= 1.15

#define _PHI_ (0.5*(1+sqrt(5)))
static stc[3]={_PHI_,1,0};
#define IN3 (1/sqrt(14+6*sqrt(5)))
static n3[3]={IN3*_PHI_,-IN3*(_PHI_^2),-IN3*(2*_PHI_+1)};

fractal(x,y,z,scale){
   
   for(i=0;i<4;i++){//Pre-fold
      y=abs(y);
      z=abs(z);
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t<0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}
   }
   r=x*x+y*y+z*z;
   for(i=0;i<MI && r<bailout;i++){//actual iterations
      y=abs(y);
      z=abs(z);
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t<0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}
      
      x=scale*x-stc[0]*(scale-1);
      y=scale*y-stc[1]*(scale-1);
      z=scale*z-stc[2]*(scale-1);
      
      r=x*x+y*y+z*z;
   }
   (sqrt(x*x+y*y+z*z)-2)*scale^(-i)
}

fractal(x,y,z,scale,rot[3][3]){
   
   for(i=0;i<5;i++){//5 pre-folds
      y=abs(y);
      z=abs(z);
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t<0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}
   }
   r=x*x+y*y+z*z;
   for(i=0;i<MI && r<bailout;i++){
#ifdef PRE_ROTATE
      x1=rot[0][0]*x+rot[1][0]*y+rot[2][0]*z;
      y1=rot[0][1]*x+rot[1][1]*y+rot[2][1]*z;
      z1=rot[0][2]*x+rot[1][2]*y+rot[2][2]*z;
      x=x1;y=y1;z=z1;
#endif
      y=abs(y);
      //z=abs(z);//I've removed this
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t<0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}
      
      x=scale*x-stc[0]*(scale-1);
      y=scale*y-stc[1]*(scale-1);
      z=scale*z-stc[2]*(scale-1);
      
      r=x*x+y*y+z*z;
   }
   (sqrt(x*x+y*y+z*z)-2)*scale^(-i)
}

fractal(x,y,z,scale,rot[3][3]){
   
   for(i=0;i<5;i++){//5 pre-folds
      y=abs(y);
      z=abs(z);x=abs(x);
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t<0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}
        
   }
   r=x*x+y*y+z*z;
   for(i=0;i<MI && r<bailout;i++){
#ifdef PRE_ROTATE
      x1=rot[0][0]*x+rot[1][0]*y+rot[2][0]*z;
      y1=rot[0][1]*x+rot[1][1]*y+rot[2][1]*z;
      z1=rot[0][2]*x+rot[1][2]*y+rot[2][2]*z;
      x=x1;y=y1;z=z1;
#endif
      y=abs(y);
      z=abs(z);x=abs(x);//added another fold
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t>0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}//changed the condition (t>0) instead of (t<0)
      
      x=scale*x-stc[0]*(scale-1);
      y=scale*y-stc[1]*(scale-1);
      z=scale*z-stc[2]*(scale-1);
      
      r=x*x+y*y+z*z;
   }
   (sqrt(x*x+y*y+z*z)-2)*scale^(-i)
}

      x=abs(x);y=abs(y);z=abs(z);
      if(x-y<0){x1=y;y=x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}

 		
int n = 0;
while (n < maxIterations) {
   p.y = abs(p.y); 
   p.z = abs(p.z);
   if (dot(p, n2)>0.0) { p *= n2Mat; }
   p.x = abs(p.x); 
   p.z = abs(p.z);
   if (dot(p, n1)>0.0) { p *= n1Mat; }        
   p = p*scale - offset*(scale-1.0);
   if (dot(p,p)> bailoutSquared) break;
   n++;
}
return (length(p) ) * pow(scale, -float(n));

phi = 1.61803399; // golden ratio.
scale = 2.0;
offset = normalize(float3(1.0,phi-1.0,0.0));
n1 = normalize(float3(-phi,phi-1.0,1.0));
n2 = normalize(float3(1.0,-phi,phi+1.0));
n1Mat = reflectionMatrix(n1);
n2Mat = reflectionMatrix(n2);
      
// Return reflection matrix for plane with normal 'n'
float3x3 reflectionMatrix(float3 n)
{
	return float3x3( 1.0 - 2.0*n.x*n.x,     - 2.0*n.y*n.x ,     - 2.0*n.z*n.x,
			      - 2.0*n.x*n.y, 1.0 - 2.0*n.y*n.y ,     - 2.0*n.z*n.y,
			      - 2.0*n.x*n.z,     - 2.0*n.y*n.z , 1.0 - 2.0*n.z*n.z );
}

      float t=dot(p,n);
      if (t>0.0) p-=2.0*t*n;

#define PHI (0.5*(1+sqrt(5))) // Phi is the golden ratio

The three points that are the vertices of a fundamental triangle.
static ftris[3]={PHI,1,0};// vertex
static ftria[3]={PHI,0,0};// middle of edge
static ftric[3]={1/3*(1+2*PHI),0,PHI/3};//center of an icosa triangle.

#define IN3 (1/sqrt(14+6*sqrt(5))) //This is the normalisation factor of the following vector
static n3[3]={IN3*PHI,-IN3*(PHI^2),-IN3*(2*PHI+1)};//this is the normal of one of the folding planes.
// The 2 others planes' normals are: {0,1,0} and {0,0,1}

//
fractal(x,y,z,scale,rot[3][3]){
   r=x*x+y*y+z*z;
   for(i=0;i<MI && r<bailout;i++){
// PRE_ROTATE
      x1=rot[0][0]*x+rot[1][0]*y+rot[2][0]*z;
      y1=rot[0][1]*x+rot[1][1]*y+rot[2][1]*z;
      z1=rot[0][2]*x+rot[1][2]*y+rot[2][2]*z;
      x=x1;y=y1;z=z1;

// Foldings. 10 in total but 7 of them are abs() ;o)
      x=abs(x);
      y=abs(y);
      z=abs(z);
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t<0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}
      y=abs(y);
      z=abs(z);
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t<0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}
      y=abs(y);
      z=abs(z);
      t=x*n3[0]+y*n3[1]+z*n3[2];
      if(t<0){x-=2*t*n3[0];y-=2*t*n3[1];z-=2*t*n3[2];}
      
// Stretching:
// for an icosahedron: stc[]=ftris[];
// for a dodecahedron: stc[]=ftric[];
// In general you can choose: stc[]=a*ftris[]+b*ftria[]+c*ftric[]; where: a+b+c=1;
      x=scale*x-stc[0]*(scale-1);
      y=scale*y-stc[1]*(scale-1);
      z=scale*z-stc[2]*(scale-1);
      
      r=x*x+y*y+z*z;
   }
   return (sqrt(x*x+y*y+z*z)-2)*(scale)^(-i);
}

// v is the input point
i=0;
r=dot(v,v);
while (i<max_iterations && r<bailout_squared) do{
      rotate1(v);
      fold(v);
      rotate2(v);
      v=scale*v-(scale-1)*offset;
      r=dot(v,v);
      i=i+1;
}
distance=sqrt(r)/scale^i;

// v is the input point
i=0;
r=dot(v,v);
while (i<max_iterations && r<bailout_squared) do{
      
      v.z=abs(v.z);//for example
      rotate2(v);//rotation
      v=v-offset*(scale-1)/scale;//translation, ((scale-1)/scale) is here just in order to keep the size of the fractal "constant" while changing the scale.
      v=scale*v;//scale
      rotate1(v);//rotation
      r=dot(v,v);
      i=i+1;
}
distance=sqrt(r)/scale^i;

 if (dot(p, n1)>0.0) { p *= n1Mat; }
 if (dot(p, n2)>0.0) { p *= n2Mat; }
 if (dot(p, n3)>0.0) { p *= n3Mat; }
 (repeated the necessary number of times)

 // The columns of the (skewed) transformation-matrix M should be the normals: n1,n2,n3
 p = M*abs(M_inverse*p);

Menger3(x,y,z){
   r=x*x+y*y+z*z;
   for(i=0;i<MI && r<bailout;i++){
      rotate1(x,y,z);

      x=abs(x);y=abs(y);z=abs(z);
      if(x-y<0){x1=y;y=x;x=x1;}
      if(x-z<0){x1=z;z=x;x=x1;}
      if(y-z<0){y1=z;z=y;y=y1;}
      
      z-=0.5*CZ*(scale-1)/scale;
      z=-abs(-z);
      z+=0.5*CZ*(scale-1)/scale;

      rotate2(x,y,z);
    
      x=scale*x-CX*(scale-1);
      y=scale*y-CY*(scale-1);
      z=scale*z;
      
      r=x*x+y*y+z*z;
   }
   return sqrt(x*x+y*y+z*z)*scale^(-i);
}

 // Prefolds.
p = abs(p);
t=dot(p,n1); if (t>0.0) { p-=2.0*t*n1; }
t=dot(p,n2); if (t>0.0) { p-=2.0*t*n2; }
  
while (n < maxIterations) {
   // folds
   p = abs(p);           
   t=dot(p,n1); if (t>0.0) { p-=2.0*t*n1; }
   t=dot(p,n2); if (t>0.0) { p-=2.0*t*n2; }
   
   // rotate1,scale,rotate2
   p4.xyz = p; p4.w = 1.0;
   p = (M*p4).xyz;
   
   n++;          
}

M = toMatrix4(rotationMatrix2) * translate(offset) * scale4x4(scale) * translate(-offset) * toMatrix4(rotationMatrix1);

for i=1 to maxiterations
    fold (whichever fold planes set)
    Multiply by An affine (or projective?) transform matrix. (in general a composition of translations, rotations and uniform scaling. You can obtain interresting results with other types of transformation but the Distance estimation must be adapted in this case)
    if radius>bailout then break
end

scale: 2
offset: -1 -1 -1
plane fold 1: 0.71 0.71 0
plane fold 2: 0.71 0 0.71
plane fold 3: 0 0.71 0.71
sphere fold: 0.25 1.00

scale: 2.5
offset: -2 -1 0
plane fold 1: 0.71 -0.71 0
plane fold 2: 0.71 0 -0.71
plane fold 3: 0 0.71 -0.71
sphere fold: 0.25 1.00