Quaternion square root

Question

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)  =>  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 code-golf rules apply. The shortest program or function in bytes in each language wins.

Answer 1

APL (NARS), 2 bytes

√

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

Answer 2

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.

Answer 3

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

Answer 4

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()

Answer 5

Haskell, 51 bytes

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!

Answer 6

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.

Answer 7

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

Answer 8

05AB1E, 14 bytes

nOtsн+·t©/¦®;š

Port of @xnor's Python 2 answer.

Try it online or verify all test cases.

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)

Answer 9

Wolfram Language (Mathematica), 28 bytes

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

Port of @xnor's Python 2 answer.

Try it online!

Answer 10

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.

Answer 11

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