# 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.

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)  =>  Output (a, b, c, d) rounded to 6 digits

0.0,  1.0,  0.0,  0.0 =>  0.707107,  0.707107,  0.000000,  0.000000
1.0,  1.0,  0.0,  0.0 =>  1.098684,  0.455090,  0.000000,  0.000000
1.0, -1.0,  1.0,  0.0 =>  1.168771, -0.427800,  0.427800,  0.000000
2.0,  0.0, -2.0, -1.0 =>  1.581139,  0.000000, -0.632456, -0.316228
1.0,  1.0,  1.0,  1.0 =>  1.224745,  0.408248,  0.408248,  0.408248
0.1,  0.2,  0.3,  0.4 =>  0.569088,  0.175720,  0.263580,  0.351439
99.0,  0.0,  0.0,  0.1 =>  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 & winning criterion

Standard rules apply. The shortest program or function in bytes in each language wins.

• Can we take the quaternion as a, (b, c, d)? – nwellnhof Nov 16 '18 at 13:19
• @nwellnhof Sure. Even something like a,[b,[c,[d]]] is fine, if you can somehow save bytes with it :) – Bubbler Nov 16 '18 at 13:50

## APL (NARS), 2 bytes

√

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

• I can't help it: you should include " ¯_(ツ)_/¯ " In your answer – Barranka Nov 16 '18 at 5:22
• You dropped this \ – Andrew Nov 16 '18 at 8:33
• @Barranka Done. – Adám Nov 16 '18 at 9:08
• @Andrew blame it on the Android app... Thank you for picking it up :) – Barranka Nov 16 '18 at 14:31
• It'd be better if it's ¯\_(⍨)√¯ – Zacharý Nov 16 '18 at 21:03

# Python 2, 72 bytes

def f(a,b,c,d):s=((a+(a*a+b*b+c*c+d*d)**.5)*2)**.5;print s/2,b/s,c/s,d/s


Try it online!

More or less a raw formula. I thought I could use list comprehensions to loop over b,c,d, but this seems to be longer. Python is really hurt here by a lack of vector operations, in particular scaling and norm.

Python 3, 77 bytes

def f(a,*l):r=a+sum(x*x for x in[a,*l])**.5;return[x/(r*2)**.5for x in[r,*l]]


Try it online!

Solving the quadratic directly was also shorter than using Python's complex-number square root to solve it like in the problem statement.

• "Input is a non-real quaternion. You can take it as four real (floating-point) numbers, in any order and structure of your choice." So you can consider it to be a pandas series or numpy array. Series have scaling with simple multiplication, and there are various ways to get norm, such as (s*s).sum()**.5. – Acccumulation Nov 16 '18 at 20:04

# Wolfram Language (Mathematica), 19 bytes

Sqrt
<<Quaternions


Try it online!

Mathematica has Quaternion built-in too, but is more verbose.

Although built-ins look cool, do upvote solutions that don't use built-ins too! I don't want votes on questions reaching HNQ be skewed.

# JavaScript (ES7), 55 53 bytes

Based on the direct formula used by xnor.

Takes input as an array.

q=>q.map(v=>1/q?v/2/q:q=((v+Math.hypot(...q))/2)**.5)


Try it online!

### How?

Given an array $$\q=[a,b,c,d]\$$, this computes:

$$x=\sqrt{\frac{a+\sqrt{a^2+b^2+c^2+d^2}}{2}}$$

And returns:

$$\left[x,\frac{b}{2x},\frac{c}{2x},\frac{d}{2x}\right]$$

q =>                            // q[] = input array
q.map(v =>                    // for each value v in q[]:
1 / q ?                     //   if q is numeric (2nd to 4th iteration):
v / 2 / q                 //     yield v / 2q
:                           //   else (1st iteration, with v = a):
q = (                     //     compute x (as defined above) and store it in q
(v + Math.hypot(...q))  //     we use Math.hypot(...q) to compute:
/ 2                     //       (q[0]**2 + q[1]**2 + q[2]**2 + q[3]**2) ** 0.5
) ** .5                   //     yield x
)                             // end of map()


f(a:l)|r<-a+sqrt(sum$(^2)<$>a:l)=(/sqrt(r*2))<$>r:l  Try it online! A direct formula. The main trick to express the real part of the output as r/sqrt(r*2) to parallel the imaginary part expression, which saves a few bytes over: 54 bytes f(a:l)|s<-sqrt$2*(a+sqrt(sum$(^2)<$>a:l))=s/2:map(/s)l


Try it online!

# Charcoal, 32 bytes

≔Ｘ⊗⁺§θ⁰ＸΣＥθ×ιι·⁵¦·⁵η≧∕ηθ§≔θ⁰⊘ηＩθ


Try it online! Link is to verbose version of code. Port of @xnor's Python answer. Explanation:

≔Ｘ⊗⁺§θ⁰ＸΣＥθ×ιι·⁵¦·⁵η


Square all of the elements of the input and take the sum, then take the square root. This calculates $$\ | x + y\vec{u} | = \sqrt{ x^2 + y^2 } = \sqrt{ (a^2 - b^2)^2 + (2ab)^2 } = a^2 + b^2 \$$. Adding $$\ x \$$ gives $$\ 2a^2 \$$ which is then doubled and square rooted to give $$\ 2a \$$.

≧∕ηθ


Because $$\ y = 2ab \$$, calculate $$\ b \$$ by dividing by $$\ 2a \$$.

§≔θ⁰⊘η


Set the first element of the array (i.e. the real part) to half of $$\ 2a \$$.

Ｉθ


Cast the values to string and implicitly print.

# Java 8, 84 bytes

(a,b,c,d)->(a=Math.sqrt(2*(a+Math.sqrt(a*a+b*b+c*c+d*d))))/2+" "+b/a+" "+c/a+" "+d/a


Port of @xnor's Python 2 answer.

Try it online.

Explanation:

(a,b,c,d)->           // Method with four double parameters and String return-type
(a=                 //  Change a to:
Math.sqrt(       //   The square root of:
2*             //    Two times:
(a+          //     a plus,
Math.sqrt(  //     the square-root of:
a*a       //      a  squared,
+b*b      //      b squared,
+c*c      //      c squared,
+d*d))))  //      And d squared summed together
/2                  //  Then return this modified a divided by 2
+" "+b/a            //  b divided by the modified a
+" "+c/a            //  c divided by the modified a
+" "+d/a            //  And d divided by the modified a, with space delimiters


# 05AB1E, 14 bytes

nOtsн+·t©/¦®;š


Port of @xnor's Python 2 answer.

Explanation:

n                 # Square each number in the (implicit) input-list
O                # Sum them
t               # Take the square-root of that
sн+            # Add the first item of the input-list
·           # Double it
t          # Take the square-root of it
©         # Store it in the register (without popping)
/        # Divide each value in the (implicit) input with it
¦       # Remove the first item
®;     # Push the value from the register again, and halve it
š    # Prepend it to the list (and output implicitly)


# Wolfram Language (Mathematica), 28 bytes

{s=#+Norm@{##},##2}/(2s)^.5&


Port of @xnor's Python 2 answer.

Try it online!

# C# .NET, 88 bytes

(a,b,c,d)=>((a=System.Math.Sqrt(2*(a+System.Math.Sqrt(a*a+b*b+c*c+d*d))))/2,b/a,c/a,d/a)


Port of my Java 8 answer, but returns a Tuple instead of a String. I thought that would have been shorter, but unfortunately the Math.Sqrt require a System-import in C# .NET, ending up at 4 bytes longer instead of 10 bytes shorter.. >.>

The lambda declaration looks pretty funny, though:

System.Func<double, double, double, double, (double, double, double, double)> f =


Try it online.

# Perl 6, 49 bytes

{;(*+@^b>>².sum**.5*i).sqrt.&{.re,(@b X/2*.re)}}


Try it online!

Curried function taking input as f(b,c,d)(a). Returns quaternion as a,(b,c,d).

### Explanation

{;                                             }  # Block returning WhateverCode
@^b>>².sum**.5     # Compute B of quaternion written as q = a + B*u
# (length of vector (b,c,d))
(*+              *i)  # Complex number a + B*i
.sqrt  # Square root of complex number
.&{                }  # Return
.re,  # Real part of square root
(@b X/2*.re)  # b,c,d divided by 2* real part
`