Sidef, 10 bytes

{|a,b|a*b}

Sidef has a built-in Quaternion type, as does Odin. Neither language versions used by TIO support this, however.

var f = {|a,b|a*b}
say f(Quaternion(12, 54, -2, 23), Quaternion(1, 4, 6, -2))
#=> Quaternion(-146, -32, 270, 331)

say f(Quaternion(1, 4, 6, -2), Quaternion(12, 54, -2, 23))
#=> Quaternion(-146, 236, -130, -333)

say f(Quaternion(3.5, 4.6, -0.24, 0), Quaternion(2.1, -3, -4.3, -12))
#=> Quaternion(10059/500, 51/25, 19823/500, -125/2)

Mathematica, 2 bytes

**

Definitely Mathematica has NonCommutativeMultiply for quaternions.

Reference: Official Documentation

Needs["Quaternions`"];

q1 = Quaternion[12, 54, -2, 23];
q2 = Quaternion[1, 4, 6, -2];
q3 = Quaternion[3.5, 4.6, -0.24, 0];
q4 = Quaternion[2.1, -3, -4.3, -12];

q1 ** q2
q2 ** q1
q3 ** q4

Quaternion[-146, -32, 270, 331]
Quaternion[-146, 236, -130, -333]
Quaternion[20.118, 2.04, 39.646, -62.5]

Go, 119 bytes

type F=float64
func f(a,b,c,d,e,f,g,h F)[]F{return[]F{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}}

Attempt This Online!

The very straightforward solution

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.

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.

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.

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}   *)

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

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.              ";

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

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]

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

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]