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

   float R  = sqrtf(x.x*x.x + x.y*x.y + x.z*x.z);
    float th = atan2(x.y, x.x);
    float ph = atan2(x.z, sqrtf(x.x*x.x + x.y*x.y));

    float3 dz;
    const float c = 1.0f;
    dz.x = 1.0f;
    dz.y = 0.0f;
    dz.z = 0.0f;
    float ph_dz = 0.0f;
    float th_dz = 0.0f;
    float R_dz  = 1.0f;

    const float sq_threshold = 4.0f;   // divergence threshold

    // Iterate to compute the distance estimator.
    int i = m_max_iterations;

    const float p = power;             // power of fractal
    while( i-- )
    {
      // derivative
      dz.x = p*powf(R, p-1.0f) * R_dz*cosf(ph_dz+(p-1.0f)*ph)*cos(th_dz+(p-1.0f)*th)+c;
      dz.y = p*powf(R, p-1.0f) * R_dz*cosf(ph_dz+(p-1.0f)*ph)*sin(th_dz+(p-1.0f)*th);
      dz.z = p*powf(R, p-1.0f) * R_dz*sinf(ph_dz+(p-1.0f)*ph                       );
      
      R_dz  = sqrtf(dz.x*dz.x + dz.y*dz.y + dz.z*dz.z);
      th_dz = atan2(dz.y, dz.x);
      ph_dz = acos(dz.z/R_dz);
      
      // new x,y,z coordinate
      x.x = powf(R, p) * cosf(p*ph)*cos(p*th) + x.x;
      x.y = powf(R, p) * cosf(p*ph)*sin(p*th) + x.y;
      x.z = powf(R, p) * sinf(p*ph)           + x.z;

      R  = sqrtf(x.x*x.x + x.y*x.y + x.z*x.z);
      th = atan2(x.y, x.x);
      ph = atan2(x.z, sqrtf(x.x*x.x + x.y*x.y));

      // Stop when we know the point diverges.
      if( R > sq_threshold )
        break;
    }
    return  0.5 * R * logf(R)/R_dz;

   float R  = sqrtf(z.x*z.x + z.y*z.y + z.z*z.z);
    float th = atan2(z.y, z.x);
    float ph = asinf(z.z / R);

    float3 dz;
    const float3 c = z;
    float ph_dz = 0.0f;
    float th_dz = 0.0f;
    float R_dz  = 1.0f;

    const float sq_threshold = 4.0f;   // divergence threshold

    // Iterate to compute the distance estimator.
    int i = m_max_iterations;

    const float p = power;             // power of fractal
    while( i-- )
    {
      // derivative
      dz.x = p*powf(R, p-1.0f) * R_dz*cosf(ph_dz+(p-1.0f)*ph)*cos(th_dz+(p-1.0f)*th)+1.0f;
      dz.y = p*powf(R, p-1.0f) * R_dz*cosf(ph_dz+(p-1.0f)*ph)*sin(th_dz+(p-1.0f)*th);
      dz.z = p*powf(R, p-1.0f) * R_dz*sinf(ph_dz+(p-1.0f)*ph                       );
      
      R_dz  = sqrtf(dz.x*dz.x + dz.y*dz.y + dz.z*dz.z);
      th_dz = atan2(dz.y, dz.x);
      ph_dz = asinf(dz.z/R_dz);
      
      // new x,y,z coordinate
      z.x = powf(R, p) * cosf(p*ph)*cos(p*th) + c.x;
      z.y = powf(R, p) * cosf(p*ph)*sin(p*th) + c.y;
      z.z = powf(R, p) * sinf(p*ph)           + c.z;

      R  = sqrtf(z.x*z.x + z.y*z.y + z.z*z.z);
      th = atan2(z.y, z.x);
      ph = asinf(z.z / R);

      // Stop when we know the point diverges.
      if( R > sq_threshold )
        break;
    }
    return  0.5 * R * logf(R)/R_dz;



And finally after optimization (sincosf etc...) the computation part looks like this:

      // derivative dz iteration
      float pP_ = p*powf(R, p-1.0f);
      float c1, s1;   sincosf(ph_dz+(p-1.0f)*ph, &s1, &c1);
      float s1_, c1_; sincosf(th_dz+(p-1.0f)*th, &s1_, &c1_);
      dz.x = pP_ * R_dz*c1*c1_+1.0f;
      dz.y = pP_ * R_dz*c1*s1_;
      dz.z = pP_ * R_dz*s1;
      
      // polar coordinates of derivative dz
      R_dz  = sqrtf(dz.x*dz.x + dz.y*dz.y + dz.z*dz.z);
      th_dz = atan2(dz.y, dz.x);
      ph_dz = asinf(dz.z /R_dz);
      
      // z iteration
      float P_ = powf(R, p);
      float s2,c2;    sincosf(p*ph, &s2, &c2);
      float s2_, c2_; sincosf(p*th, &s2_, &c2_);
      z.x = P_ * c2*c2_ + c.x;
      z.y = P_ * c2*s2_ + c.y;
      z.z = P_ * s2    + c.z;

      // polar coordinates of z
      R  = sqrtf(z.x*z.x + z.y*z.y + z.z*z.z);
      th = atan2(z.y, z.x);
      ph = asinf(z.z / R);