Thanks for the corrections.
My power 5 looks less mutant now. And some other weird effects (such as banding effects in flat surface regions) are gone. In addition to the errors you found, I was not adding a constant c vector in each iteration - instead I was adding the current x vector in some kind of recursive feedback loop. I am surprised it even gave sane images.
Here's the cleaned up code where I also renamed x to z.
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);