Welcome to Fractal Forums

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

i*i=-1
(i*i)*j=-j   // multiply by j on the right
i*(i*j)=-j   // multiplication is associative 
i*k=-j       // use the fact that j*i = i*j = k (commutativity)

j*i=k
j*j*i=j*k    // multiply by j on the left
-i*i=j*k     // use j*j = -i
j*k=1        // use i*i = -1

j*i=k
j*i*j*i = k*k  // multiply by j*i on left side, k on right side
j*k*i = k*k    // i*j = k
-i*i = k*k     
k*k = 1

nx=sqr(sx)-sqr(sy)+2*sz*sk;
ny=2*sx*sy-sqr(sz)+sqr(sk);
nz=2*sx*sz-2*sy*sk;
nk=2*sx*sk+2*sy*sz;
				
if (juliaMode) {
		sx=nx+cr;
		sy=ny+ci;
		sz=nz+cj;
		sk=nk+ck;
} else {		
		sx=nx+(pixelr);
		sy=ny+(pixeli);
		sz=nz+(pixelj);
		sk=nk+(pixelk);
}

#include "DE-Raytracer.frag"
#group Hypercomplex Thing

// Number of fractal iterations.
uniform int Iterations;  slider[0,16,100]
// Breakout distance
uniform float Threshold; slider[0,10,100]
// Mandel or Julia
uniform bool JuliaMode; checkbox[false]

// Quaterion Constant (first three components - only used for Julia)
uniform vec3 C123; slider[(-1,-1,-1),(0.18,0.88,0.24),(1,1,1)]
// Quaterion Constant (last component)
uniform float C4; slider[-1,0.16,1]

vec4 c = vec4(C123,C4); // We don't support 4-component sliders yet...

void init() {}

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

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* mul(p,dp);
		p = mul(p,p) + (JuliaMode ? c : 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);
}