How to define a new algebra using TagSet?

Question

For example, I would like to define the algebra of quaternions. I coded something like

Unprotect[i, j, k]; ClearAll[i, j, k]
i /: i ** i := -1;
i /: i ** j := k;
i /: j ** i = -k;
i /: i ** k := -j;
i /: k ** i = j;
j /: j ** j := -1;
j /: j ** k := i;
j /: k ** j := -i;
k /: k ** k := -1;
Protect[i, j, k];
ExpandNCM[(h : NonCommutativeMultiply)[a___, b_Plus, c___]] := 
  Distribute[h[a, b, c], Plus, h, Plus, ExpandNCM[h[##]] &];
ExpandNCM[(h : NonCommutativeMultiply)[a___, b_Times, c___]] := 
  Most[b] ExpandNCM[h[a, Last[b], c]];
ExpandNCM[a_] := ExpandAll[a];
QQ = Collect[ExpandNCM[#], {i, j, k}] &;

However, I got two different results when calculating

(a i + b j + c k) ** (d i + e j + f k) // QQ
(a i + b j + c k) ** (l i + m j + n k) // QQ

The first one gave the expected result

k (a e-b d)+j (c d-a f)-a d+i (b f-c e)-b e-c f

while the second one gave a weird answer

i (a i**l+b j**l+c k**l)+j (a i**m+b j**m+c k**m)+k (a i**n+b j**n+c k**n)

I hope the NonCommutativeMultiply ** only operates on those i,j,k. I don't know where the problem is.

I also want to know if there is a handy way to define a new algebra?

Updated

Now the algebra has been set up for +, *, ** but not Power.

Here is the code

Unprotect[i, j, k]; ClearAll[i, j, k, ijkQ, ExpandNCM, Qconj, Qreal, \
Qimag]
i /: i ** i := -1;
i /: i*i := -1;
i /: i ** j := k;
i /: j ** i = -k;
i /: i ** k := -j;
i /: k ** i = j;
j /: j ** j := -1;
j /: j*j := -1;
j /: j ** k := i;
j /: k ** j := -i;
k /: k ** k := -1;
k /: k*k := -1;

i /: Qconj[i] := -i;
j /: Qconj[j] := -j;
k /: Qconj[k] := -k;
Protect[i, j, k];

Qconj[a_: 0] := a /; FreeQ[a, i | j | k];
Qconj[(b_: 1)*(q : i | j | k)] := b*Qconj[q] /; FreeQ[b, i | j | k];
Qconj[a_ + b_] := Qconj[a] + Qconj[b];

Qreal := Select[# + i, (FreeQ[#, i | j | k] &)] &;
Subscript[Qimag, x_] := ({0}~Join~
     Cases[{#}, (a : (_?(FreeQ[#, i | j | k] &)) : 1)*x :> a, 
      Infinity])[[-1]] &

ExpandNCM[(h : NonCommutativeMultiply)[a___, b_Plus, c___]] := 
  Distribute[h[a, b, c], Plus, h, Plus, ExpandNCM[h[##]] &];

ExpandNCM[(h : NonCommutativeMultiply)[a___, 
    b_?(FreeQ[#, i | j | k] &)*(q : i | j | k), c___]] := 
  b*ExpandNCM[h[a, q, c]];
ExpandNCM[(h : NonCommutativeMultiply)[a___, 
    b_?(FreeQ[#, i | j | k] &), c___]] := b*ExpandNCM[h[a, c]];

ExpandNCM[(h : (Times | Plus))[a___, b_NonCommutativeMultiply, 
    c___]] := h[a, c, ExpandNCM[b]];
ExpandNCM[(h : NonCommutativeMultiply)[a_: 1]] := a;

ExpandNCM[a_] := ExpandAll[a];

QQ = Collect[ExpandNCM[#], {i, j, k}] &;

Now the problem is how to set up the rules for Power, such as Exp[] , Log[] and Power[]. For example, I want to define Power[q,n] as q**q** ... **q and then use it to define Exp[q]. How should I do?

Note that the algebra of quaternions here is just an example. What I want is to define the rules for the algebra and let the Mathematica do the calculations. These rules will be changed.

Answer 1

This is an answer for the "I don't know where the problem is" question. The problem is this definition

ExpandNCM[(h:NonCommutativeMultiply)[a___,b_Times,c___]]:=Most[b]*ExpandNCM[h[a,Last[b],c]];

and more precisely the pattern b. To understand the problem, enter

d*i
(* d * i *)

m*i
(* i * m *)

In the second case, Times has reordered the variables and put i before m. This is because Times is Orderless and i comes before m in the canonical order.

Therefore, if the b_Times in the definition above is matched with m*i then

Most[b] == i
Last[b] == m

and the i is moved out of the NoncommutativeMultiply.

One possible fix is to replace the definition above by (but please clear old definitions before doing this):

ExpandNCM[(h:NonCommutativeMultiply)[a___,b_*(q:i|j|k),c___]]:=b*ExpandNCM[h[a,q,c]];

Now

(a i+b j+c k)**(l i+m j+n k) // QQ
(* -a l-b m+k (-b l+a m)-c n+j (c l-a n)+i (-c m+b n) *)

Answer 2

UPDATE

Just found out that there's a Quaternions package: http://reference.wolfram.com/language/Quaternions/tutorial/Quaternions.html

ORIGINAL ANSWER

Here's a start:

(* some basic helpers *)
Quaternion[a_] := Quaternion[a, 0, 0, 0];
Quaternion[a_, b_] := Quaternion[a, b, 0, 0];
Quaternion[a_, b_, c_] := Quaternion[a, b, c, 0];
QI = Quaternion[0, 1];
QJ = Quaternion[0, 0, 1];
QK = Quaternion[0, 0, 0, 1];

Now a bit of meat:

Quaternion /: 
  NonCommutativeMultiply[Quaternion[a1_, b1_, c1_, d1_], Quaternion[a2_, b2_, c2_, d2_]] :=
    Quaternion[
      a1 a2 - b1 b2 - c1 c2 - d1 d2, 
      a1 b2 + b1 a2 + c1 d2 - d1 c2, 
      a1 c2 - b1 d2 + c1 a2 + d1 b2, 
      a1 d2 + b1 c2 - c1 b2 + d1 a2]

Some checks:

Quaternion[3, 0, 0, 0] ** Quaternion[7, 0, 0, 0]
(* Quaternion[21, 0, 0, 0] *)

QJ ** QK
(* Quaternion[0, 1, 0, 0] *)

QI ** QK
(* Quaternion[0, 0, -1, 0] *)

QI ** QJ ** QK
(* Quaternion[-1, 0, 0, 0] *)

Here's a simplistic formatting rule. This is just for demonstration--you should make this more robust.

Format[Quaternion[a_, b_, c_, d_]] := a + b i + c j + d k

Which gives,

Quaternion[1, 1, 1, 1] ** Quaternion[1, 2, 3, 4]
(* dispays as -8 + 4 i + 2 j + 6 k *)

Quaternion[1, 1, 1, 1] ** Quaternion[1, 2, 3, 4] // FullForm
(* Quaternion[-8, 4, 2, 6] *)

Note

I'm only vaguely familiar with quaternions and non-commutative multiplication, so there may be semantic errors above. It should all be made robust and then tested.