More explanation of the code:

First of all a "func" is your own definition of a function call so this one is called by the main code as test3D7(a,b) where a and b are any complex variables though for this particular instance on entry the imaginary part of b is always zero.

The ampersands (&) in the function definition of entry value means that the complex a,b are passed by reference i.e. so within the func "z" actually refers directly to the passed variable "a" and w to b rather than a,b being copied into separate variables z,w - obviously pass by reference is faster and no return values are required at the end of the func.

Quaternions can be treated as a real plus a 3D vector i.e. (a,b,c,d) can be treated as (s,V) where s is just a real and V is the 3D vector made up of (b,c,d). Also for the purposes of this case it's generally OK to treat a 3D system as a restricted quaternion where the vector is just a 2D one, or indeed simply treated as a complex number i.e. (a,b,c) treated as (s,V) where V is made-up of (b,c). In addition in both cases the magnitude of V can be treated as the imaginary part of a complex number (s,m) where m is the vector's magnitude.

The code first gets the square of the full 3D magnitude - avoiding the unnecessary square root calculation at this point.

Then we only need to do further processing on non-zero values.

On the line declaring s there is also a re-assignment of w so that it contains the full complex vector made up of the imaginary part of z as the real part and what was the real part of w as the imaginary and s is assigned the magnitude of this (I've since realised that I should again be just calculating the square here and put the root calculation after the if test…)

if s is non-zero then we need to process the vector and first the real part of z is copied into r which I've now realised isn't strictly necessary at this point

Now we raise w to the power required by the angle scaling parameter and multiply by the magnitude to the power of the difference between the magnitude power parameter and the angle scaling parameter resulting in a new value for w such that the angle is scaled by bulba but the magnitude is raised to power bulbm. Normally bulbm should be 1 so that the magnitude of the vector is unchanged but the angle has been scaled as required.

If the value of bulbm is not 1 then we have to recalculate both m (the full 3D magnitude squared) and s (the magnitude of the 2d vector w), here for the new s raising a float to a power will be faster than using |w| whenever the power parameter is actually an integer (smart Frederik compiling).

In the next section the quaternionic calculations are done, for the power 2 version simply giving (x,y,z)^2 as (x^2-y^2-z^2, 2xy, 2xz) and for other powers using the classic angle scale and magnitude power method but obviously for quaternions this time.

Before the endfunc of course the first two parts are back in z and the third in the real part of w with the imaginary part of w being zero.

If the vector was zero then obviously we simply raise the real part to the power.

NB. edited to correct "bulbp" to "bulbm" and to mention that on exit of the func the imaginary part of w is zero....