Multiply Quaternions 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
Python (83)

r=lambda A,B,R=range(4):[sum(A[m]*B[m^p]*(-1)**(14672&gt;&gt;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.