Welcome to Fractal Forums

vec4 mul(vec4 a, vec4 b) {
	return vec4( a.x*b.x, a.x*b.y + a.y*b.x, a.x*b.z+a.z*b.x,  a.x*b.w+a.y*b.z+a.z*b.y+a.w*b.x);
}

float DE(vec3 pos) {
	vec4 p = vec4(pos,C4);
	vec4 dp = vec4(1.0,0.0,0.0,0.0);
	for (int i = 0; i < Iterations; i++) {
		dp = 2.0* mul(dp,p) + vec4(1.0,0.0,0.0,0.0);
		p = mul(p,p) + vec4(pos,C4);
		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);
}

This formula was inspired by:
http://www.fractalforums.com/general-discussion-b77/hyper-dual-numbers/
It is mathematically meaningless - but can give cool fractals

 // original formula:
 // a.x*b.x, a.x*b.y + a.y*b.x, a.x*b.z+a.z*b.x,  a.x*b.w+a.y*b.z+a.z*b.y+a.w*b.x
 // was too unregular - especially x component so I fixed
 
 a=x, b=y, c=z, d=w,  pow = IntPower-1;
 do { // do allows you to choose a higher pow
 newa = xa - yb - ( zc + wd ) - ( zd + wc )*zwMul; // instead of xa
 // if you remove yb and ( zc + wd ) you get weird things...
 newb = xb + ya;
 newc = xc + za + ( yd + wb )*ywMul;
 newd = xd + wa + ( yc + zb )*yzMul;
 pow--, a=newa, ... d=newd;
 } while pow>=0;
 x=a+Cx, y=b+Cy, z=c+Cz, w=d+Cw