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

Questions
1

In the Quaternion package tutorial, it says:

A four-dimensional analog of de Moivre’s theorem is used for calculating powers of quaternions
In[15]:=
Quaternion[1, 2, 0, 1]^2.5
Out[15]= Quaternion[-9.0604, 2.20741, 0., 1.1037]

That’s fine, but if I take the 2.5 root again, I get

In[16]:= Quaternion[-9.0604, 2.20741, 0., 1.1037]^2.5
Out[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?

2

Hi Marc,

this is a bug

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]:= QuaternionQ[(Quaternion[1, 2, 0, 1]^(2.5))^(2.5)]
Out[3]= False

but, normally, if one speaks about the root of a number, one has to mean the principal root (cmp. Principal Square Root -- from Wolfram MathWorld).
This seems to be implemented in the Quaternions package

In[4]:= (Quaternion[1, 2, 0, 1]^(1/2))^(5) // N
Out[4]= Quaternion[-9.0604, 2.20741, 0., 1.1037]

In[5]:= (Quaternion[1, 2, 0, 1]^5)^(1/2) // N
Out[5]= Quaternion[9.0604, -2.20741, 0., -1.1037]

try

Quaternion[1,2,0,1]^(1/n) // N

for integers n. Therefore one has simply to modify the rules for Power in the package

Quaternion[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 > 0

one 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 > 0

finally (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]= True

be aware that the arithmetics of quaternions is not too similar to the arithmetics of complex numbers. For example there is a continuum of square roots
of -1 in the skew field of quaternions showing another weakness of the Quaternions package

In[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]= I

again, note that it says I, not -I …

Regards
Udo.

3

Hi again,

my proposal does not work in some cases. The strange thing with the original Quaternions package is that one can do things like that:

In[1]:= Needs["Quaternions`"]

In[2]:= (* Orig *) Quaternion[-1, 2, -7, 1]^(1/9.)
Out[2]= Quaternion[1.15911 + 0.421881 I, -0.0507501 - 0.0184715 I, 
 0.177625 + 0.0646503 I, -0.025375 - 0.00923576 I]

In[3]:= % ** % ** % ** % ** % ** % ** % ** % ** % 
Out[3]= Quaternion[-1. + 2.24647*10^-16 I, 
 2. - 8.1532*10^-16 I, -7. + 2.81719*10^-15 I, 1. - 4.0766*10^-16 I]

despite the fact, that Out[2] gives a Bi-Quaternion, non-commutative multiplication brings the original quaternion back. It is necessary to check the Quaternions package in more detail. Here is a reference

http://www.mathematica-journal.com/issue/v8i3/features/turner/contents/html/Links/index\_lnk\_3.html

“Real Quaternions and Rotations” by Robert Piziak and Danny W. Turner, The Mathematica Journal, vol. 8 issue 3.

Using the article’s formulae one gets the ninth zeroes of Quaternion[-1, 2, -7, 1] as follows

In[64]:= m = 9;
q = Quaternion[-1, 2, -7, 1];
w = PolarQtrig[q]
w /. deg -> Degree // Simplify
phi = angle[q]

Out[66]= Sqrt[55] (Cos[
    97.7494 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
      3 Sqrt[6])) Sin[97.7494 deg])
Out[67]= (-1. + 2. I) - 7. J + 1. K
Out[68]= 97.7494 deg

In[69]:= r = Abs[q]^(1/m)
Out[69]= 55^(1/18)

In[72]:= theta = (phi/m + d 360. deg/m);
n = AdjustedSignIJK[q];
ninthrootsofq = 
 Table[Flatten[
   r (Cos[theta] (Quaternion[1, 0, 0, 0] // FromQuaternion) + 
      Sin[theta] (n // FromQuaternion))], {d, 0, m - 1}]
numroots = ninthrootsofq /. deg -> Degree // N // Simplify

Out[74]= {55^(1/18) (Cos[0. + 10.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[0. + 10.861 deg]), 
 55^(1/18) (Cos[50.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[50.861 deg]), 
 55^(1/18) (Cos[90.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[90.861 deg]), 
 55^(1/18) (Cos[130.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[130.861 deg]), 
 55^(1/18) (Cos[170.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[170.861 deg]), 
 55^(1/18) (Cos[210.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[210.861 deg]), 
 55^(1/18) (Cos[250.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[250.861 deg]), 
 55^(1/18) (Cos[290.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[290.861 deg]), 
 55^(1/18) (Cos[330.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[330.861 deg])}

Out[75]= {(1.22698 + 0.0640715 I) - 0.22425 J + 0.0320357 K, 
          (0.788599 + 0.263735 I) - 0.923072 J + 0.131867 K, 
          (-0.0187746 + 0.339994 I) - 1.18998 J + 0.169997 K, 
          (-0.817363 + 0.257166 I) - 0.90008 J + 0.128583 K, 
          (-1.2335 + 0.0540071 I) - 0.189025 J + 0.0270036 K, 
          (-1.07247 - 0.174422 I) + 0.610477 J - 0.087211 K, 
          (-0.409615 - 0.321237 I) + 1.12433 J - 0.160619 K, 
          (0.4449 - 0.317742 I) + 1.1121 J - 0.158871 K, 
          (1.09124 - 0.165572 I) + 0.579501 J - 0.0827858 K}

(* only real coefficients, let's check the reproduction of q *)
Remove[q1, q2, q3, q4, q5, q6, q7, q8, q9]
q1 = ToQuaternion[%75[[1]]];
q2 = ToQuaternion[%75[[2]]];
q3 = ToQuaternion[%75[[3]]];
q4 = ToQuaternion[%75[[4]]];
q5 = ToQuaternion[%75[[5]]];
q6 = ToQuaternion[%75[[6]]];
q7 = ToQuaternion[%75[[7]]];
q8 = ToQuaternion[%75[[8]]];
q9 = ToQuaternion[%75[[9]]];

In[95]:= q1 ** q1 ** q1 ** q1 ** q1 ** q1 ** q1 ** q1 ** q1
Out[95]= Quaternion[-1., 2., -7., 1.]

In[96]:= q2 ** q2 ** q2 ** q2 ** q2 ** q2 ** q2 ** q2 ** q2
Out[96]= Quaternion[-1., 2., -7., 1.]

In[97]:= q3 ** q3 ** q3 ** q3 ** q3 ** q3 ** q3 ** q3 ** q3
Out[97]= Quaternion[-1., 2., -7., 1.]

In[98]:= q4 ** q4 ** q4 ** q4 ** q4 ** q4 ** q4 ** q4 ** q4
Out[98]= Quaternion[-1., 2., -7., 1.]

In[99]:= q5 ** q5 ** q5 ** q5 ** q5 ** q5 ** q5 ** q5 ** q5
Out[99]= Quaternion[-1., 2., -7., 1.]

In[100]:= q6 ** q6 ** q6 ** q6 ** q6 ** q6 ** q6 ** q6 ** q6
Out[100]= Quaternion[-1., 2., -7., 1.]

In[101]:= q7 ** q7 ** q7 ** q7 ** q7 ** q7 ** q7 ** q7 ** q7
Out[101]= Quaternion[-1., 2., -7., 1.]

In[102]:= q8 ** q8 ** q8 ** q8 ** q8 ** q8 ** q8 ** q8 ** q8
Out[102]= Quaternion[-1., 2., -7., 1.]

In[103]:= q9 ** q9 ** q9 ** q9 ** q9 ** q9 ** q9 ** q9 ** q9
Out[103]= Quaternion[-1., 2., -7., 1.]

Finally two solutions are possilbe
(a) you stay with the original Quaternions package but find a clever way to transform its Bi-Quaternions to Quaternions
(b) you use my proposal supplemented with the formulae of Piziak and Turner for the (1/q) part

Sincerely
Udo.