Multiply Quaternions

Question

Write a named function or program that computes the quaternion product of two quaternions. Use as few bytes as possible.

Quaternions

Quaternions are an extension of the real numbers that further extends the complex numbers. Rather than a single imaginary unit i, quaternions use three imaginary units i,j,k that satisfy the relationships.

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

(There's also tables of these on the Wikipedia page.)

In words, each imaginary unit squares to -1, and the product of two different imaginary units is the remaining third one with a +/- depending on whether the cyclic order (i,j,k) is respected (i.e., the right-hand-rule). So, the order of multiplication matters.

A general quaternion is a linear combination of a real part and the three imaginary units. So, it is described by four real numbers (a,b,c,d).

x = a + b*i + c*j + d*k

So, we can multiply two quaternions using the distributive property, being careful to multiply the units in the right order, and grouping like terms in the result.

(a + b*i + c*j + d*k) * (e + f*i + g*j + h*k)
= (a*e - b*f - c*g - d*h)    +
  (a*f + b*e + c*h - d*g)*i  +
  (a*g - b*h + c*e + d*f)*j  +
  (a*h + b*g - c*f + d*e)*k

Seen this way, quaternion multiplication can be seen as a map from a pair of 4-tuples to a single 4-tuples, which is what you're asked to implement.

Format

You should write either a program or named function. A program should take inputs from STDIN and print out the result. A function should take in function inputs and return (not print) an output.

Input and output formats are flexible. The input is eight real numbers (the coefficients for two quaternions), and the output consists of four real numbers. The input might be eight numbers, two lists of four numbers, a 2x4 matrix, etc. The input/output format don't have to be the same.The ordering of the (1,i,j,k) coefficients is up to you.

The coefficients can be negative or non-whole. Don't worry about real precision or overflows.

Banned: Function or types specifically for quaternions or equivalents.

Test cases

These are in (1,i,j,k) coefficient format.

[[12, 54, -2, 23], [1, 4, 6, -2]] 
 [-146, -32, 270, 331]

[[1, 4, 6, -2], [12, 54, -2, 23]] 
 [-146, 236, -130, -333]

[[3.5, 4.6, -0.24, 0], [2.1, -3, -4.3, -12]] 
 [20.118, 2.04, 39.646, -62.5]

Reference Implementation

In Python, as function:

#Input quaternions: [a,b,c,d], [e,f,g,h]
#Coeff order: [1,i,j,k]

def mult(a,b,c,d,e,f,g,h):
    coeff_1 = a*e-b*f-c*g-d*h
    coeff_i = a*f+b*e+c*h-d*g
    coeff_j = a*g-b*h+c*e+d*f
    coeff_k = a*h+b*g-c*f+d*e

    result = [coeff_1, coeff_i, coeff_j, coeff_k]
    return result

Answer 1

CJam, 49 45 39 bytes

"cM-^\M-^G-^^KM-zP"256bGbq~m*f{=:*}4/{:-W*}/W*]`

The above uses caret and M notation, since the code contains unprintable characters.

At the cost of two additional bytes, those characters can be avoided:

6Z9C8 7YDXE4BFA5U]q~m*f{=:*}4/{:-W*}/W*]`

You can try this version online: CJam interpreter

Test cases

To calculate (a + bi + cj + dk) * (e + fi + gj + hk), use the following input:

[ d c b a ] [ h g f e ]

The output will be

[ z y x w ]

which corresponds to the quaternion w + xi + yj + zk.

$ base64 -d > product.cjam <<< ImOchy0eS/pQIjI1NmJHYnF+bSpmez06Kn00L3s6LVcqfS9XKl1g
$ wc -c product.cjam
39 product.cjam
$ LANG=en_US cjam product.cjam <<< "[23 -2 54 12] [-2 6 4 1]"; echo
[331 270 -32 -146]
$ LANG=en_US cjam product.cjam <<< "[-2 6 4 1] [23 -2 54 12]"; echo
[-333 -130 236 -146]
$ LANG=en_US cjam product.cjam <<< "[0 -0.24 4.6 3.5] [-12 -4.3 -3 2.1]"; echo
[-62.5 39.646 2.04 20.118]

How it works

6Z9C8 7YDXE4BFA5U]  " Push the array [ 6 3 9 12 8 7 2 13 1 14 4 11 15 10 5 0].         ";
q~                  " Read from STDIN and interpret the input.                         ";
m*                  " Compute the cartesian product of the input arrays.               ";
f                   " Execute the following for each element of the first array:       ";
{                   " Push the cartesian product (implicit).                           ";
    =               " Retrieve the corresponding pair of coefficients.                 ";
    :*              " Calculate their product.                                         ";
}                   "                                                                  ";
4/                  " Split into chunks of 4 elements.                                 ";
{:-W*}/             " For each, subtract the first element from the sum of the others. ";
W*                  " Multiply the last integers (coefficient of 1) by -1.             ";
]`                  " Collect the results into an array and stringify it.              ";

Answer 2

Python (83)

r=lambda A,B,R=range(4):[sum(A[m]*B[m^p]*(-1)**(14672>>p+4*m)for m in R)for p in R]

Takes two lists A,B in [1,i,j,k] order and returns a result in the same format.

The key idea is that with [1,i,j,k] corresponding to indices [0,1,2,3], you get the product's index (up to sign) by XOR'ing the indices. So, the terms that get placed in index p are those who indices XOR to p, and are thus the products A[m]*B[m^p].

It only remains to make the signs work out. The shortest way I found was to simply code them into a magic string. The 16 possibilities for (m,p) are turned into numbers 0 to 15 as p+4*m. The number 14672 in binary has 1's at the places where -1 signs are needed. By shifting it the appropriate number of places , a 1 or 0 winds up at the last digit, making the number odd or even, and so (-1)** is either 1 or -1 as needed.

Answer 3

Python - 90 75 72 69

Pure Python, no libraries - 90:

m=lambda a,b,c,d,e,f,g,h:[a*e-b*f-c*g-d*h,a*f+b*e+c*h-d*g,a*g-b*h+c*e+d*f,a*h+b*g-c*f+d*e]

It's probably pretty hard to shorten this "default" solution in Python. But I'm very curious as to what other might come up with. :)

---

Using NumPy - 75 72 69:

Well, since the input and output are rather flexible, we can use some NumPy functions and exploit the scalar-vector representation:

import numpy
m=lambda s,p,t,q:[s*t-sum(p*q),s*q+t*p+numpy.cross(p,q)]

Input arguments s and t are the scalar parts of the two quaternions (the real parts) and p and q are the corresponding vector parts (the imaginary units). Output is a list containing scalar part and vector part of the resulting quaternion, the latter being represented as NumPy array.

Simple test script:

for i in range(5):
    a,b,c,d,e,f,g,h=np.random.randn(8)
    s,p,t,q=a, np.array([b, c, d]), e, np.array([f, g, h])
    print mult(a, b, c, d, e, f, g, h), "\n", m(s,p,t,q)

(mult(...) being the OP's reference implementation.)

Output:

[1.1564241702553644, 0.51859264077125156, 2.5839001110572792, 1.2010364098925583] 
[1.1564241702553644, array([ 0.51859264,  2.58390011,  1.20103641])]
[-1.8892934508324888, 1.5690229769129256, 3.5520713781125863, 1.455726589916204] 
[-1.889293450832489, array([ 1.56902298,  3.55207138,  1.45572659])]
[-0.72875976923685226, -0.69631848934167684, 0.77897519489219036, 1.4024428845608419] 
[-0.72875976923685226, array([-0.69631849,  0.77897519,  1.40244288])]
[-0.83690812141836401, -6.5476014589535243, 0.29693969165495304, 1.7810682337361325] 
[-0.8369081214183639, array([-6.54760146,  0.29693969,  1.78106823])]
[-1.1284033842268242, 1.4038096725834259, -0.12599103441714574, -0.5233468317643214] 
[-1.1284033842268244, array([ 1.40380967, -0.12599103, -0.52334683])]

Answer 4

Haskell, 85

m a b c d e f g h=[a*e-b*f-c*g-d*h,a*f+b*e+c*h-d*g,a*g-b*h+c*e+d*f,a*h+b*g-c*f+d*e]

Porting it to Haskell saves us a few chars ;)

Answer 5

Mathematica 83 50

Probably can be golfed more..

p = Permutations;
f = #1.(Join[{{1, 1, 1, 1}}, p[{-1, 1, -1, 1}][[1 ;; 3]]] p[#2][[{1, 8, 17, 24}]]) &

Spaces and newlines not counted & not needed.

Usage:

f[{a,b,c,d},{e,f,g,h}]        (* => {x,w,y,z}   *)

EDIT How this works.

The Mathematica function Permutations makes all possible permutations of #2 (the second argument). There are 24 permutations, but we only need {e,f,g,h}, {f,e,h,g}, {g,h,e,f}, and {h,g,f,e}. These are the first, 8th, 17th and 24th permutations. So the code

p[#2][[{1,8,17,24}]]

exactly selects these from the permutations of the second argument and returns them as a matrix. But then they don't have the correct sign yet. The code p[{-1,1,-1,1}][[1;;3]] returns a 3x4 matrix with the correct sign. We prepend it with {1,1,1,1} by using Join, and making a normal multiplication (Times, or as is the case here by just writing them after each other) between two matrices makes an element-by-element multiplication in Mathematica.

So finally, the result of

(Join[{{1, 1, 1, 1}}, p[{-1, 1, -1, 1}][[1 ;; 3]]] p[#2][[{1, 8, 17, 24}]])

is the matrix

 e  f  g  h
-f  e -h  g
-g  h  e -f
-h -g  f  e

Making a matrix multiplication between {a,b,c,d} (the first argument #1) and the former matrix gives the desired result.

EDIT 2 Shorter code

Inspired by the Python code of Falko, I split up the quaternion in a scalar and a vector part, and use Mathematica's built in command Cross to calculate the cross product of the vector parts:

f[a_, A_, b_, B_] := Join[{a*b - A.B}, a*B + b*A + Cross[A, B]]

Usage:

f[a,{b,c,d},e,{f,g,h}]        (* => {x,w,y,z}   *)

Answer 6

Whispers v2, 396 bytes

> 1
> 2
> 0
> 4
> Input
> Input
>> 6ᶠ2
>> 6ᵗ2
>> 7ⁿ3
>> 7ⁿ1
>> 10‖9
>> 8ⁿ3
>> 8ⁿ1
>> 13‖12
>> 7‖8
>> 11‖14
>> 8‖7
>> 14‖11
>> 15‖16
>> 19‖17
>> 20‖18
>> 4⋅5
>> L⋅R
>> Each 23 22 21
> [1,-1,-1,-1,1,1,1,-1,1,-1,1,1,1,1,-1,1]
>> Each 23 24 25
>> 26ᶠ4
>> 26ᵗ4
>> 28ᶠ4
> 8
>> 26ᵗ30
>> 31ᶠ4
>> 31ᵗ4
>> ∑27
>> ∑29
>> ∑32
>> ∑33
>> Output 34 35 36 37

Try it online!

Takes input in the form

[a, b, c, d]
[e, f, g, h]

and outputs as

w
x
y
z

to represent \$q = w + xi + yj + zk\$

The structure tree of this answer is:

A good chunk of this answer comes from two main faults in Whispers:

• No function to reverse an array
• The usage of sets in the computation of the Cartesian product

Therefore, we can split the code into 3 sections.

How it works

We'll use the following definitions for clarity and conciseness:

$$q = a + bi + cj + dk$$ $$p = e + fi + gj + hk$$ $$r = w + xi + yj + zk, \: (q \cdot p = r)$$ $$A = [a, b, c, d]$$ $$B = [e, f, g, h]$$ $$C = [w, x, y, z]$$

Section 1: Permuting \$A\$ and \$B\$

The first section is by far the longest, stretching from line 1 to line 22:

> 1
> 2
> 0
> 4
> Input
> Input
>> 6ᶠ2
>> 6ᵗ2
>> 7ⁿ3
>> 7ⁿ1
>> 10‖9
>> 8ⁿ3
>> 8ⁿ1
>> 13‖12
>> 7‖8
>> 11‖14
>> 8‖7
>> 14‖11
>> 15‖16
>> 19‖17
>> 20‖18
>> 4⋅5

The main purpose of this section is to permute \$B\$ so that simple element-wise multiplication between \$A\$ and \$B\$ is possible. There are four different arrangements of \$B\$ to multiply the elements of \$A\$ with:

$$B_1 = [e, f, g, h]$$ $$B_2 = [f, e, h, g]$$ $$B_3 = [g, h, e, f]$$ $$B_4 = [h, g, f, e]$$

The second input, \$B\$, is stored on line 6. We then split \$B\$ down the middle, as each possible arrangement of \$B\$ is grouped in pairs. In order to reverse these pairs (to get the correct orders in \$B_2\$ and \$B_4\$), we take the first and last element, then concatenate them in reverse order:

>> 7ⁿ3
>> 7ⁿ1
>> 10‖9

(forming \$[f, e]\$) and

>> 8ⁿ3
>> 8ⁿ1
>> 13‖12

(forming \$[h, g]\$). We now have all the halves needed to form the arrangements, so we concatenate them together to form \$B_1, B_2, B_3\$ and \$B_4\$. Finally, we concatenate these four arrangements together to form \$B_T\$. We then quickly make \$A_T\$, defined as \$A\$ repeated \$4\$ times:

$$A_T = [a, b, c, d, a, b, c, d, a, b, c, d, a, b, c, d]$$ $$B_T = [e, f, g, h, f, e, h, g, g, h, e, f, h, g, f, e]$$

When each element of \$B_T\$ is multiplied by the corresponding element in \$A_T\$, we get the (signless) values in \$q \cdot p\$

Section 2: Signs and products

As said in Section 1, the values in \$A_T\$ and \$B_T\$ correspond to the signless (i.e. positive) values of each of the coefficients in \$q \cdot p\$. No obvious pattern is found in the signs that would be shorter than simply hardcoding the array, so we hardcode the array:

> [1,-1,-1,-1,1,1,1,-1,1,-1,1,1,1,1,-1,1]

We'll call this array \$S\$ (signs). Next, we zip together each element in \$A_T, B_T\$ and \$S\$ and take the product of each sub-array e.g. \$[[a, e, 1], [b, f, -1], \dots, [e, f, -1], [d, e, 1]] \to D = [ae, -bf, \dots, -ef, de]\$.

Section 3: Partitions and final sums.

Once we have the array of coefficients of \$q \cdot p\$, with signs, we need to split it into 4 parts (i.e. the four factorised coefficients of \$q \cdot p\$), and then take the sums. This leads us to the only golfing opportunity found: moving the

> 4

to line 4 rather than 26, as it is used 6 times, each time saving a byte by moving it. Unfortunately, this costs a byte changing the 9 to a 10, so only \$5\$ bytes are saved. The next section takes slices of size \$4\$ from the front of \$D\$, saving each slice to the corresponding row and passing on the shortened list of \$D\$. Once 4 slices are taken, we the take the sum of each, before outputting them all.

Answer 7

Python, 94

The most straightforward way isn't too long.

def m(a,b,c,d,e,f,g,h):return[a*e-b*f-c*g-d*h,a*f+b*e+c*h-d*g,a*g-b*h+c*e+d*f,a*h+b*g-c*f+d*e]

Answer 8

JavaScript ES6 - 86

f=(a,b,c,d,e,f,g,h)=>[a*e-b*f-c*g-d*h,a*f+b*e+c*h-d*g,a*g-b*h+c*e+d*f,a*h+b*g-c*f+d*e]

Answer 9

Lua - 99

Might as well.

_,a,b,c,d,e,f,g,h=unpack(arg)print(a*e-b*f-c*g-d*h,a*f+b*e+c*h-d*g,a*g-b*h+c*e+d*f,a*h+b*g-c*f+d*e)

Lua's "unpack()" frees the elements of a table. So the table 'arg' is where all the command line input is stored (including arg[0] which is the program's file name, it gets discarded).

Answer 10

Python, 58 56 chars

m=lambda x,y,z,w:(x*z-y*(2*w.real-w),x*w+y*(2*z.real-z))

I take very liberal use of the input/output format wiggle room. The inputs are 4 complex numbers, encoded thusly:

x = a+b*i
y = c+d*i
z = e+f*i
w = g+h*i

It outputs a pair of complex numbers in a similar format, the first of the pair encodes the real and i part, the second encodes the j and k parts.

To see this works, note that the first quaternion is x+y*j and the second is z+w*j. Just evaluate (x+y*j)*(z+w*j) and realize that j*t = conj(t)*j for any imaginary number t.