Adding quaternion numbers
Just add the real and imaginary components independently as follows:
(a + i b + j c + k d)+(e + i f + j g + k h) = (a+e) + i ( b+f) + j ( c+g) + k ( d+h)
This operation will be coded in the sfrotation class (see this class here).
Alternatively, if we want to use scalar and vector notation for quaternions, as defined on this page then addition is:
(sa,va) + (sb,vb) = (sa+sb,va+vb)
where:
 (sa,va) = quaternion a
 (sb,vb) = quaternion b
Subtracting quaternion numbers
Just subtract the real and imaginary components independently as follows:
(a + i b + j c + k d)(e + i f + j g + k h) = (ae) + i ( bf) + j ( cg) + k ( dh)
This operation will be coded in the sfrotation class (see this class here).
Alternatively, if we want to use scalar and vector notation for quaternions, as defined on this page then subtraction is:
(sa,va)  (sb,vb) = (sasb,vavb)
where:
 (sa,va) = quaternion a
 (sb,vb) = quaternion b
Multiplying quaternion numbers
The multiplication rules for the imaginary operators are the same as for other numbers. We just put each quaternion in brackets and multiply out all the terms: (a + i b + j c + k d)*(e + i f + j g + k h). When we are multiplying the imaginary operators we use the following rules:
 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
Note that the order of multiplication is significant, in other words q1 * q2 is not necessarily equal to q2 * q1, we might expect this because quaternions can be used to represent rotations and the order of rotations is significant, for example, if you rotate 90 degrees about the xaxis and then 90 degrees about the yaxis you get a different result than if you rotate about the yaxis then the xaxis. See here for more information about the rules of rotation.
In mathematical terms quaternion rotation is not necessarily commutative, some subsets of quaternions may commute, for example:
 real part: a*e = e*a
 real times imaginary: a * i f = i f * a
 same imaginary part: i b * i f = i f * i b ( = b*f )
Other subsets of quaternions anticommute, in other words reversing the order multiplies the result by 1. This apples when different imaginary quantities, as can be seen from the list above, so for example:
 i*j =  j*i = k
In the general case quaternions containing all these parts would be neither commutative nor anticommutative. In order to understand more about the effect of reversing the operands it is useful to introduce the concept of the conjugate explained here. If the quaternion is unit length (normalised) as when used to represent rotations:
conj(q1 * q2) = conj(q2) * conj(q1)
In other words, to reverse the order, take the conjugate of each part and then take the conjugate of the whole thing.
The following table summarises these issues, it shows the type and sign for each permutation of parts of the multiplying quaternions.
a*b 
b.1  b.i  b.j  b.k 
a.1  1  i  j  k 
a.i  i  1  k  j 
a.j  j  k  1  i 
a.k  k  j  i  1 
This table summarises the character of quaternions. We can compare it with, say, geometric algebra.
Notice that squaring any linear combination of imaginary terms will give a negative real number because the leading diagonal is negative and the other terms cancel out because they anticommute.
So back to the general case of multiplying any two quaternions, we can just expand out the terms and group as follows:
(a + i b + j c + k d)*(e + i f + j g + k h)  =a(e + i f + j g + k h) + i b(e + i f + j g + k h) + j c(e + i f + j g + k h) + k d(e + i f + j g + k h) 

= 


= 


= 
a*e  b*f  c*g  d*h 
so the result is:
z1 * z2= a*e  b*f  c*g d*h + i (b*e + a*f + c*h  d*g) + j (a*g  b*h + c*e + d*f) + k (a*h + b*g  c*f + d*e)
notes:
 Other quaternion multiplication issues.
 This operation will be coded in the sfrotation class (see this class here).
Alternatively, if we want to use scalar and vector notation for quaternions, as defined on this page then multiplication is:
(sa,va) * (sb,vb) = (sa*sbva•vb,va × vb + sa*vb + sb*va)
where:
 (sa,va) = quaternion a
 (sb,vb) = quaternion b
 • = vector dot product
 × = vector cross product
Division
We don't tend to use the notation for division, since quaternion multiplication is not commutative we need to be able to distinguish between q1*q2^{1} and q2^{1}*q1. So instead of a divide operation we multiply by the inverse.
If the quaternion is unit length (normalised, as it will be if we are using quaternions to represent rotations) then:
q^{1} = conj(q)
or
q * conj(q) = 1
If the quaternion is not unit length then divide the conjugate by a scalar value which is the square of the magnitude of the quaternion:
q^{1} = conj(q) / q^{2}
So expanding out the terms gives:
(a + i b + j c + k d)^{1}= (a  i b  j c  k d) / (a^{2} + b^{2} + c^{2} + d^{2})
We can easily prove this because:
q * conj(q) = q^{2}
Which we can check by expanding out the terms.
Alternatively, if we want to use scalar and vector notation for quaternions, as defined on this page then division is:
(sa,va) / (sb,vb) = (sa*sb+va•vb,va x vb  sa*vb + sb*va)
where:
 (sa,va) = quaternion a
 (sb,vb) = quaternion b
 • = vector dot product
 x = vector cross product