Quaternion square root Background

Quaternion is a number system that extends complex numbers. A quaternion has the following form

$$ a + bi + cj + dk $$

where \$ a,b,c,d \$ are real numbers and \$ i,j,k \$ are three fundamental quaternion units. The units have the following properties:

$$ i^2 = j^2 = k^2 = -1 $$
$$ ij = k, jk = i, ki = j $$
$$ ji = -k, kj = -i, ik = -j $$

Note that quaternion multiplication is not commutative.

Task

Given a non-real quaternion, compute at least one of its square roots.

How?

According to this Math.SE answer, we can express any non-real quaternion in the following form:

$$ q = a + b\vec{u} $$

where \$ a,b\$ are real numbers and \$ \vec{u} \$ is the imaginary unit vector in the form \$ xi + yj + zk \$ with \$ x^2 + y^2 + z^2 = 1 \$. Any such \$ \vec{u} \$ has the property \$ \vec{u}^2 = -1 \$, so it can be viewed as the imaginary unit.

Then the square of \$ q \$ looks like this:

$$ q^2 = (a^2 - b^2) + 2ab\vec{u} $$

Inversely, given a quaternion \$ q' = x + y\vec{u} \$, we can find the square root of \$ q' \$ by solving the following equations

$$ x = a^2 - b^2, y = 2ab $$

which is identical to the process of finding the square root of a complex number.

Note that a negative real number has infinitely many quaternion square roots, but a non-real quaternion has only two square roots.

Input and output

Input is a non-real quaternion. You can take it as four real (floating-point) numbers, in any order and structure of your choice. Non-real means that at least one of \$ b,c,d \$ is non-zero.

Output is one or two quaternions which, when squared, are equal to the input.

Test cases

   Input (a, b, c, d)  =&gt;  Output (a, b, c, d) rounded to 6 digits

 0.0,  1.0,  0.0,  0.0 =&gt;  0.707107,  0.707107,  0.000000,  0.000000
 1.0,  1.0,  0.0,  0.0 =&gt;  1.098684,  0.455090,  0.000000,  0.000000
 1.0, -1.0,  1.0,  0.0 =&gt;  1.168771, -0.427800,  0.427800,  0.000000
 2.0,  0.0, -2.0, -1.0 =&gt;  1.581139,  0.000000, -0.632456, -0.316228
 1.0,  1.0,  1.0,  1.0 =&gt;  1.224745,  0.408248,  0.408248,  0.408248
 0.1,  0.2,  0.3,  0.4 =&gt;  0.569088,  0.175720,  0.263580,  0.351439
99.0,  0.0,  0.0,  0.1 =&gt;  9.949876,  0.000000,  0.000000,  0.005025


Generated using this Python script. Only one of the two correct answers is specified for each test case; the other is all four values negated.

Scoring &amp; winning criterion

Standard code-golf rules apply. The shortest program or function in bytes in each language wins.
APL (NARS), 2 bytes

√

NARS has built-in support for quaternions. ¯\_(⍨)_/¯