DE Formula for Quaternion + folding

[Patryk Kizny]

I played a bit with the original quaternion adding iteration rotation, box fold and sphere fold.
The question I have is regarding the accuracy of DE formula in this case.

If the 'original' quaternion looked as follows:

Code:

float DE(vec3 pos) {
	vec4 p = vec4(pos, 0.0);
	vec4 dp = vec4(1.0, 0.0,0.0,0.0);
	
for (int i = 0; i < Iterations; i++) {
		dp = 2.0* vec4(p.x*dp.x-dot(p.yzw, dp.yzw), p.x*dp.yzw+dp.x*p.yzw+cross(p.yzw, dp.yzw));
		p = vec4(p.x*p.x-dot(p.yzw, p.yzw), vec3(2.0*p.x*p.yzw)) + C;

		float p2 = dot(p,p);
		orbitTrap = min(orbitTrap, abs(vec4(p.xyz,p2)));
		if (p2 > Threshold) break;
	}
	float r = length(p);
	return  0.5 * r * log(r) / length(dp);
}

And roughly we wanted to tweak as follows:

Code:

float DE(vec3 pos) {
	vec4 p = vec4(pos, 0.0);
	vec4 dp = vec4(1.0, 0.0,0.0,0.0);
	
for (int i = 0; i < Iterations; i++) {
		dp = 2.0* vec4(p.x*dp.x-dot(p.yzw, dp.yzw), p.x*dp.yzw+dp.x*p.yzw+cross(p.yzw, dp.yzw));
		p = vec4(p.x*p.x-dot(p.yzw, p.yzw), vec3(2.0*p.x*p.yzw)) + C;

====> roatate(p);
====> boxfold(p);
====> spherefold(p);

		float p2 = dot(p,p);
		orbitTrap = min(orbitTrap, abs(vec4(p.xyz,p2)));
		if (p2 > Threshold) break;
	}
	float r = length(p);
	return  0.5 * r * log(r) / length(dp);
}

So we get something like

Code:

float DEQT(vec3 pos) {
	
	vec4 p = vec4(pos, 0.0);
	vec4 dp = vec4(1.0, 0.0,0.0,0.0);
	
	for (int i = 0; i <QTIterations; i++){
		
		// Quaternion iterations
		 dp = 2.0* vec4(p.x*dp.x-dot(p.yzw, dp.yzw), p.x*dp.yzw+dp.x*p.yzw+cross(p.yzw, dp.yzw));
		 p = vec4(p.x*p.x-dot(p.yzw, p.yzw), vec3(2.0*p.x*p.yzw)) + QTC;
		
		//rotation
		p.xyz *= rot;

		//  Box Folding
		if (foldBoxIts > 0){
			if ( i >= foldBoxStartIt && i < foldBoxStartIt + foldBoxIts){				
				//  Box Folding (v2) - separate axes + offset
				p.x = clamp(p.x, -foldingLimitX, foldingLimitX) * foldingValueX - p.x + Trans.x;
				p.y = clamp(p.y, -foldingLimitY, foldingLimitY) * foldingValueY - p.y + Trans.y;
				p.z = clamp(p.z, -foldingLimitZ, foldingLimitZ) * foldingValueZ - p.z + Trans.z;					
			}
		}
		// Sphere Folding
		if(foldSphereIts > 0){
			// Sphere folding & scale
			if ( i >= foldSphereStartIt && i < foldSphereStartIt + foldSphereIts){
				float r2 = dot(p.xyz, p.xyz);
				p *= clamp(max(MinRad2/r2, MinRad2), 0.0, 1.0);
				p = p*scale + vec4(Julia,0);
			}
		}

		float p2 = dot(p,p);
		orbitTrap = min(orbitTrap, abs(vec4(p.xyz,p2)));
		if (p2 > QTThreshold) break;
		
	}
	
	//QT
	float r = length(p);
	return  0.5 * r * log(r) / length(dp);
	
}

So we got p, we got dp then we introduce a few other transformations on p while leaving dp untouched for the final DE calculation.
Is the DE formula still valid if we do so? I have a feeling this is not exactly accurate and that might be a reason for much noise I am getting out of these images.

My math skills are really below my aspirations, so wrapping my head around all this gives me a hard time.
Just started thinking about it after having read syntopia's blogposts on DE fractals.