How to define a new algebra using TagSet?

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.

• Rather than try to define the algebra directly on the special symbols i,j,k, I'd suggest you define a special wrapper Quaternion[...]. this would be analogous to Complex[...] (Note that FullForm[a+b I] is Plus[a, Times[Complex[0, 1], b]] and just FullForm[I] is Complex[0,1]. Then you would define your operations for pairs of Quaternion elements. For your special symbols, you could define i = Quaternion[0,1,0,0] and similarly for j and k. Oct 19, 2022 at 23:16
• I should add that you can then use Format to provide a friendly display form for your Quaternion[...]s. Oct 19, 2022 at 23:28
• Hi, there is an analogous question and answer with the algebra of dual numbers here. The first answer there, at the present time, uses the quaternion package to make an algebra using TagSet Oct 20, 2022 at 0:21
• Sweet! The Quaternions package should supercede my answer. It's nice, though, to see that I went down a similar path. Oct 20, 2022 at 0:24
• Just adding a point on what @lericr mentioned, FullForm[a+b I] is Plus[a, Times[Complex[0, 1], b]] and not Complex[a,b] because a and b could be complex themselves. FullForm[1+2*I] is Complex[1,2] as expected because 1 and 2 are real. Oct 20, 2022 at 0:33

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


UPDATE

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

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.

• There is also a Quaternion resource function by Wolfram staff here. Oct 20, 2022 at 0:39
• Thanks! The algebra of quaternions here is only an example. I want to define the rules for the algebra and let Mathematica do the calculations so that I don't need to define everything for the algebra. Oct 22, 2022 at 17:43
• Isn't that exactly what I illustrated? Oct 22, 2022 at 18:19
• No, the quaternionic product is already done by hand and given as a function 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]. Suppose we have a different rule for the algebra, say i is the cube root of the -1. Then you have to redefine those quantities inside Quaternion[...] Oct 22, 2022 at 18:42
• You're going to have to define the rules somewhere. Oct 22, 2022 at 19:54