How can I avoid complex components in a quaternion? (Seems like a bug)

GROUPS:
 In the Quaternion package tutorial, it says:A four-dimensional analog of de Moivre's theorem is used for calculating powers of quaternionsIn[15]:= Quaternion[1, 2, 0, 1]^2.5Out[15]= Quaternion[-9.0604, 2.20741, 0., 1.1037]That's fine, but if I take the 2.5 root again, I getIn[16]:= Quaternion[-9.0604, 2.20741, 0., 1.1037]^2.5Out[16]:=Quaternion[ 6.51112*10^-14 + 212.669 I, -4.56519*10^-14 - 149.11 I, -1.39769*10^-30 - 4.56519*10^-15 I, -2.2826*10^-14 - 74.5552 I]The quaternion components are now complex (e.g. the first component above is 6.51112*10^-14 + 212.669 I). This seems to occur when the real part of the quaternion is negative. This seems like a bug. The output ought to be another well formed quaternion.I thought I could simplify this by going:Quaternion[ 6.51112*10^-14 + 212.669 I, -4.56519*10^-14 - 149.11 I, -1.39769*10^-30 - 4.56519*10^-15 I, -2.2826*10^-14 - 74.5552 I]= (6.51112*10^-14 + 212.669 I) + (-4.56519*10^-14 - 149.11 I)*I + (-1.39769*10^-30 - 4.56519*10^-15 I)*J + (-2.2826*10^-14 - 74.5552 I)*K=  6.51112*10^-14 + 212.669 I -4.56519*10^-14 I - 149.11 I^2  - 1.39769*10^-30 J - 4.56519*10^-15 IJ  - 2.2826*10^-14 K - 74.5552 IK= 6.51112*10^-14 + 212.669 I -4.56519*10^-14 I + 149.11  - 1.39769*10^-30 J - 4.56519*10^-15 K - 2.2826*10^-14 K + 74.5552 J= (6.51112*10^-14 + 149.11) + (212.669 - 4.56519*10^-14) I  + (74.5552 - 1.39769*10^-30) J - (4.56519*10^-15 + 2.2826*10^-14) K= Quaternion[149.11, 212.669, 74.5552, -2.73911*10^-14]But this is incorrect as Quaternion[149.11, 212.669, 74.5552, -2.73911*10^-14]^0.4 <> Quaternion[-9.0604, 2.20741, 0., 1.1037]What is going on here? How can I either avoid getting complex components when I take quaternion roots or eliminate the complex components when they appear?
 Udo Krause 1 Vote Hi Marc, this is a bugIn[1]:= Needs["Quaternions"]In[2]:= Quaternion[1, 2, 0, 1]^2.5Out[2]= Quaternion[-9.0604, 2.20741, 0., 1.1037]In[3]:= QuaternionQ[(Quaternion[1, 2, 0, 1]^(2.5))^(2.5)]Out[3]= Falsebut, normally, if one speaks about the root of a number, one  has   to mean the principal root (cmp. http://mathworld.wolfram.com/PrincipalSquareRoot.html).This seems to be implemented in the Quaternions packageIn[4]:= (Quaternion[1, 2, 0, 1]^(1/2))^(5) // NOut[4]= Quaternion[-9.0604, 2.20741, 0., 1.1037]In[5]:= (Quaternion[1, 2, 0, 1]^5)^(1/2) // NOut[5]= Quaternion[9.0604, -2.20741, 0., -1.1037]tryQuaternion[1,2,0,1]^(1/n) // Nfor integers n. Therefore one has simply to modify the rules for Power in the packageQuaternion[a,b,c,d]^(p/q) => (Quaternion[a,b,c,d]^p)^(1/q)adding to the original ones (* *************** Rules for Power *************** *) (* Power is essentially de Moivre's theorem for quaternions. *)  Quaternion /:     Power[a:Quaternion[__?ScalarQ],0]:= 1  Quaternion /:     Power[a:Quaternion[__?ScalarQ],n_]:= Power[1/a, -n] /; n < 0 && n != -1 Quaternion /:    Power[a:Quaternion[__?ScalarQ],-1]:= Conjugate[a] ** 1/Norm[a]Quaternion /:    Power[a:Quaternion[__?ScalarQ], n_]:=    Module[{angle, pqradius},    angle = toAngle[a];    pqradius = angle[[2]]^n Sin[n angle[[3]]];    Quaternion[        angle[[2]]^n Cos[n angle[[3]]],        Cos[angle[[4]]] pqradius,        Sin[angle[[4]]] Cos[angle[[5]]] pqradius,        Sin[angle[[4]]] Sin[angle[[5]]] pqradius    ]  // QSimplify    ] /; ScalarQ[n] && n > 0one more rule and modifying another one (* *************** Rules for Power *************** *) (* Power is essentially de Moivre's theorem for quaternions. *)  Quaternion /:     Power[a:Quaternion[__?ScalarQ],0]:= 1  Quaternion /:     Power[a:Quaternion[__?ScalarQ],n_]:= Power[1/a, -n] /; n < 0 && n != -1 Quaternion /:    Power[a:Quaternion[__?ScalarQ],-1]:= Conjugate[a] ** 1/Norm[a]Quaternion /:    Power[a:Quaternion[__?ScalarQ], n_]:=    Module[{angle, pqradius, r},    angle = toAngle[a];    pqradius = angle[[2]]^n Sin[n angle[[3]]];    Quaternion[        angle[[2]]^n Cos[n angle[[3]]],        Cos[angle[[4]]] pqradius,        Sin[angle[[4]]] Cos[angle[[5]]] pqradius,        Sin[angle[[4]]] Sin[angle[[5]]] pqradius    ]  // QSimplify    ] /; ScalarQ[n] && (Numerator[Rationalize[n]]==1 || Denominator[Rationalize[n]]==1)Quaternion /:    Power[a:Quaternion[__?ScalarQ], n_]:=Block[{r=Rationalize[n], y},        y=If[Element[n,Rationals],1,1.,1.];        Power[Power[a,Numerator[r]],y/Denominator[r]]    ] /; ScalarQ[n] && n > 0finally (after restarting Mma and reloading the package) arriving at In[1]:= Needs["Quaternions"]  In[2]:= Quaternion[1, 2, 0, 1]^2.5 Out[2]= Quaternion[9.0604, -2.20741, 0., -1.1037]  In[3]:= (Quaternion[1, 2, 0, 1]^2.5)^(2.5) Out[3]= Quaternion[212.669, -149.11, 4.56519*10^-15, -74.5552]  In[6]:= QuaternionQ[(Quaternion[1, 2, 0, 1]^2.5)^(2.5)]Out[6]= Truebe aware that the arithmetics of quaternions is not too similar to the arithmetics of complex numbers. For example there is a continuum of square rootsof -1 in the skew field of quaternions showing another weakness of the Quaternions packageIn[4]:= Sqrt[Quaternion[-1, 0, 0, 0]]During evaluation of In[4]:= Power::infy: Infinite expression 1/0 encountered. >>During evaluation of In[4]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>Out[4]= Quaternion[Sqrt[\[Pi]/2], 0, Indeterminate, Indeterminate]Compare this to the well-known In[5]:= Sqrt[-1]Out[5]= Iagain, note that it says I, not -I  ...RegardsUdo.