First one was very discontinuous. Images at the bottom.... messed around a bit in between:
UPDATE: I also tried your old {x,y,z,w}2 = {x2-y2-z2-w2, 2(xy+zw), 2(xz+yw), 2(xw+yz)} set with the signs... nothing interesting: fractured all up. Then I did something a bit freaky:
New code:
nx=(sqr(sx)-r1); /nx is new x sx is starting x ... same for all variables
ny=-2*abs(sx*sy+sz*sk); // abs is absolute value
nz=2*abs(sx*sz+sy*sk); //
nk=-2*abs(sx*sk+sy*sz); //
Guess what I got? Dun dun dun dun.... a square 3d BS variant:
Notice the similarity to my old fractal (1st image is square formula, 2nd image is my old one, formula follows):

new x= x^2 -y^2 -z^2 + x pixel value
new y= - abs [ 2* sqrt(x^2+y^2) * z] + absolute value of y pixel value
new z= - abs [(x^2+y^2-z^2) * (2*x*y) * (x^2+y^2)^-1] + absolute value of z pixel value
Your original code:
nx=sqr(sx)-sqr(sy) - sqr(sz) - sqr(sk);
if (sy<0) { sy2=-1;} else {sy2=1;}
if (sz<0) {sz2=-1;} else {sz2=1;}
if (sk<0) {sk2=-1;} else {sk2=1;}
ny=2*sz2*sk2*sx*sy;
nz=2*sy2*sk2*sx*sz;
nk=2*sy2*sz2*sx*sk;
then pixel or julia values...


