Simplify quaternion multiplication, version 1.0
peter.vlasschaert@gmail.com,23/10/2018

kill(all)$

load(lrats)$
load(vect)$

Maxima to declare a non-commutative variable?
=> use 'dot'
========================================
rem : i.j ≠ j.i , non-commutative product
rem : <a,b >= a*b-b*a ≠ 0.

1e) ordinary multiplication
---------------------------------------

p1:h=i*j+j*i;

2e)  combination of the dot and multiplication '1e'
--------------------------------------------------------------------------

p1a:a*j.i+b*i.j;

rule : dot 

p2:subst([i.j=k,j.i=-k,i.i=-1],a*j.i+b*i.j+c*i.i);

declare([a1,b1,c1,d1,a2,b2,c2,d2],scalar);

p:a1+b1*i+c1*j+d1*k;
q:a2+b2*i+c2*j+d2*k;

 multiplication of two quaternions.
=============================

hamilton : i.i =j.j=k.k=i.j.k=-1
rem : . =  dot operator.

rules : non-commutative for quarternions
--------------------------------------------------------------
[i . j=k,j . i=-k,j . k=i,k . j=-i,k . i=j,i . k=-j,i . i=-1,j . j=-1,k . k=-1]

p3:p.q;

p4:expand(p3);

p5:[i.j=k,j.i=-k,j.k=i,k.j=-i,k.i=j,i.k=-j,i.i=-1,j.j=-1,k.k=-1];
p6:lratsubst(p5,p4);

p7:factorsum(p6);

p8:[a1=1,b1=1,c1=1,d1=1,a2=1,b2=1,c2=1,d2=1];

p9:lratsubst(p8,p7);

product : '1+i+j+k' with '1+i+j+k'
==========================

p10:lratsubst(p5,p9);

 symbolic : computation : z ='p*q' :{p,q} are quaternions
===============================================

p11:args(p7);

p11a:p7;

real part : 'p*q'
============
============

p12a:[i=0,j=0,k=0];
p12b:subst(p12a,p11a);

rem : RE(p) = a1
           RE(q) = a2
p12b = RE(p)*RE(q)  - dotproduct(v1,v2)
==================================
vector  'p' = b1*i+c1*j+d1*k    => v1
vector  'q' = b2*i+c2*j+d2*k    => v2

vector part : 'p*q'
==============
==============

vector  'p' = b1*i+c1*j+d1*k    => v1
vector  'q' = b2*i+c2*j+d2*k    => v2

=>  ' result : vector 'r' = b3*i+c3*j+d3*k '

  find :b3, coef from , 'i'
====================

p13:coeff(p11a,i);

 find :c3, coef from , 'j'
====================

p14:coeff(p11a,j);

 find :d3, coef from , 'k'
====================

p15:coeff(p11a,k);

rem : RE(p)*v2+RE(q)*v1+vectorproduct(v1,v2)

total : result   'z=p*q'
=================

result : [ real-part, r = vector-part ] =[a3,b3,c3,d4]
====================================

p16:r=[p13,p14,p15];

p16a:z = [p12b,p13,p14,p15];

z=[RE(p)*RE(q)  - dotproduct(v1,v2),RE(p)*v1+RE(q)*v2+vectorproduct(v1,v2)]
=================================================================
=================================================================

===============================================================================
RE(z)           =      -d1*d2-c1*c2-b1*b2+a1*a2
vector(z)    =  [c1*d2-c2*d1+a1*b2+a2*b1,-b1*d2+b2*d1+a1*c2+a2*c1,a1*d2+a2*d1+b1*c2-b2*c1]
===============================================================================
===============================================================================
rem : RE =real part 
rem : dotproduct = scalarproduct = innerproduct
rem :vectorproduct = crossproduct = outerproduct.
===============================================================================