If *y* is a quaternion, how many solutions does *x*² = *y* have? That is, does every quaternion have a square root, and if so, how many square roots does it have?

A quaternion can have more than two roots. There is a example right in the definition of quaternions:

*i*² = *j*² = *k*² = -1

and so -1 has at least three square roots. It follows that every negative real number has at least three square roots.

## Negative reals

In fact, every square root of -1 corresponds to a square root of any other negative number. To see this, let *p* be a positive real number and let *r* = √*p* > 0. If *x* is a square root of -1, then *rx* is a square root of –*p*. Conversely, if *x* is a square root of –*p*, then *x*/*r* is a square root of -1.

So all negative real numbers have the same number of square roots in the quaternions, and that number is at least 3. In fact, it turns out that all negative real numbers have infinitely many quaterion roots. A calculation shows that

(*a* + *bi* + *cj* + *dk*)² = -1

if and only if *b*²* + c*²* + d*² = 1. So there are as many square roots of -1 as there are points on a sphere.

## Not negative reals

If a quaternion is not a negative real number, then the story is much more familiar: the equation *x*² = *y* has one solution if *y* = 0 and two solutions otherwise.

To summarize, *x*² = *y* has

- one solution if
*y*= 0, - infinitely many solutions if
*y*is negative and real, and - two solutions otherwise.

It’s easy to show that 0 has only one root: |*x*²| = |*x*|², and so if *x**² *= 0, then |*x*| = 0 and so *x* = 0. Here the norm of a quaternion is defined by

|*a* + *bi* + *cj* + *dk| = *√(*a²* + *b²* + *c²* + *d²*).

If *p* > 0, the solutions to *x*² = –*p* are *bi* + *cj* + *dk* where *b*²* + c*²* + d*² = *p*.

What’s left is to show that if *y* is not zero and not a negative real, then *x*² = *y* has two solutions.

## Polar form and calculating roots

Quaternions have a polar form analogous to the polar form for complex numbers. That is, we can write a quaternion *y* in the form

*y* = *r*(cos θ + *u* sin θ)

where *r* is real and positive, θ is real, and *u* is a unit quaternion with zero real part. We can find *r*, θ, and *u* as follows. Write *y* as

*a* + *bi* + *cj* + *dk* = a + *q*

so that *a* is the real part and *q* is the pure quaternion part. Then we can find *r*, θ, and *u* by

*r* = |*y*|

cos θ = a/|*y*|

*u* = *q* / |*q*|

Then the two roots are

*x* = ±√*r *(cos θ/2 + *u* sin θ/2).

## Python code

First we’ll implement the algorithm above, then show that it produces the same results as the Python package `quaternionic`

.

First, our code to compute *x* satisfying *x*² = *y*.

import numpy as np np.random.seed(20210106) def norm(q): return sum(t*t for t in q)**0.5 y = np.random.random(size=4) r = norm(y) theta = np.arccos(y[0]/r) u = y + 0 # make a copy u[0] = 0 u /= norm(u) x = np.sin(theta/2)*u x[0] = np.cos(theta/2) x *= r**0.5

And now we check our results using `quaternionic`

.

import quaternionic qx = quaternionic.array(x) qy = quaternionic.array(y) print(f"y = {y}") print(f"x*x = {qx*qx}") # should equal y print(f"x = {x}") # our root print(f"x = {np.sqrt(qy)}") # root via package

This produces

y = [0.61615367 0.07612092 0.09606777 0.11150865] x*x = [0.61615367 0.07612092 0.09606777 0.11150865] x = [0.79189641 0.04806243 0.06065678 0.07040609] x = [0.79189641 0.04806243 0.06065678 0.07040609]

## More quaternion posts

- Finnegans Wake and quaternions
- How far is
*xy*from*yx*on average for quaternions? - Probability of commuting

For an exploration of the roots of general polynomial equations over polynomials, see Equations in Quaternions by Ivan Niven. The American Mathematical Monthly, Dec., 1941, Vol. 48, pp. 654-661.