Programming Puzzles & Code Golf Stack Exchange is a question and answer site for programming puzzle enthusiasts and code golfers. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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


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.


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
share|improve this question
up vote 4 down vote accepted

CJam, 49 45 39 bytes


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.              ";
share|improve this answer

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.

share|improve this answer
The XOR part is pure genius. – Dennis Aug 29 '14 at 18:57

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):
    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.)


[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])]
share|improve this answer

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

share|improve this answer

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.


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


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]]


f[a,{b,c,d},e,{f,g,h}]        (* => {x,w,y,z}   *)
share|improve this answer
Could you explain how this works? What are 1, 8, 17, 24? – xnor Aug 28 '14 at 21:33

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]
share|improve this answer

JavaScript ES6 - 86

share|improve this answer

Lua - 99

Might as well.


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

share|improve this answer

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.

share|improve this answer
Very clever! Do you know why the quaternions seem to multiply like complex numbers with complex coefficients, as it looks from your expression? – xnor Aug 30 '14 at 4:13
Never mind, I now understand from your explanation how i and j act like as inner and outer complex coefficients. How fascinating! – xnor Aug 30 '14 at 4:21
It's funny how the conj calls take up more than 2/5 of your chars. I think you can shave a char each using (2*w.real-w). abs(w)**2/w would work but for 0. Maybe even exec with string substitution would be worth it? ` – xnor Aug 30 '14 at 4:37

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.