# How to calculate Gamma functions values for Quaternions using Mathematica?

How can we calculate the values for the Gamma function using Quaternions on Mathematica?

For example:

<<Quaternions; (* load package *)

Gamma[5] = 24

N[Gamma[I]] = -0.15495 - 0.498016 I

N[Gamma[1 + 2 I]] = 0.151904 + 0.0198049 I

N[Gamma[Quaternion[5, 0, 0, 0]]] = Gamma[Quaternion[5., 0., 0., 0.]] (* which is wrong *)


Specially, how can we calculate

Gamma[Quaternion[1, 2, 3, 4]]?


Attempts:

N[Gamma[FromQuaternion[Quaternion[1, 2, 0, 0]]]] (* works *)
N[Gamma[FromQuaternion[Quaternion[1, 2, 3, 4]]]] (* doesn't work *)

QuaternionPower[x_, y_] := E^(y ** Log[x])
QuaternionGamma[z_] :=
Integrate[QuaternionPower[E, -x + Log[x]*(z - 1)], {x, 0, Infinity}]
N[QuaternionGamma[FromQuaternion[Quaternion[1, 2, 0, 0]]]] =
0.151904 + 0.0198049 I (* works *)
N[QuaternionGamma[FromQuaternion[Quaternion[1, 2, 3, 0]]]] =
ConditionalExpression[Gamma[(1.  + 2. I) + 3. J], Re[J] > -0.333333] (* doesn't work *)

• It works for me in 12.1. What version are you using? May 4, 2020 at 23:59
• o.O, I'm using 12.1 as well, what is the result that you've found? May 5, 2020 at 0:01
• Ugh, it returns Gamma[(1. + 2. I) + 3. J + 4. K] for me and I thought it's some special form. May 5, 2020 at 0:30
• Hum, I see. Unforunately it is not really the case, (1. + 2. I) + 3. J + 4. K is just the open form of N[Quaternion[1, 2, 3, 4]]. May 5, 2020 at 0:32

If the function to be applied to a quaternion is holomorphic/meromorphic, one way to proceed is to switch to the complex matrix representation of a quaternion (see e.g. this), and then apply MatrixFunction[] to it.

For instance,

<< Quaternions

BlockRandom[SeedRandom[42]; (* for reproducibility *)
qq = Quaternion @@ RandomVariate[NormalDistribution[], 4]]
Quaternion[-1.01205, 0.826338, -1.39379, 0.41692]

Tan[qq] (* elementary example *)
Quaternion[-0.0652529, 0.50833, -0.857402, 0.256472]

Quaternion @@ Flatten[ReIm[First[MatrixFunction[Tan,
{{qq[[1]] + I qq[[2]], qq[[3]] + I qq[[4]]},
{-qq[[3]] + I qq[[4]], qq[[1]] - I qq[[2]]}}]]]]
Quaternion[-0.0652529, 0.50833, -0.857402, 0.256472]


Thus, to evaluate a gamma function for quaternion argument:

Quaternion @@ Flatten[ReIm[First[MatrixFunction[Gamma,
{{qq[[1]] + I qq[[2]], qq[[3]] + I qq[[4]]},
{-qq[[3]] + I qq[[4]], qq[[1]] - I qq[[2]]}}]]]]
Quaternion[-0.0568728, 0.0210734, -0.0355446, 0.0106324]


Now, for a warning. The reason I emphasized the requirement of analyticity is that if the function has branch cuts, then the complex matrix method might give results that are the conjugate of the expected result, due to possible differences in branch cut choice. For instance,

Log[qq] (* a consistent example *)
Quaternion[0.670588, 1.04449, -1.76174, 0.526985]

Quaternion @@ Flatten[ReIm[First[MatrixFunction[Log,
{{qq[[1]] + I qq[[2]], qq[[3]] + I qq[[4]]},
{-qq[[3]] + I qq[[4]], qq[[1]] - I qq[[2]]}}]]]]
Quaternion[0.670588, 1.04449, -1.76174, 0.526985]

q2 = Reverse[qq] (* inconsistent example *)

Log[q2]
Quaternion[0.670588, 0.989246, -0.586497, 0.718306]

(* note that this returns the conjugate *)
Quaternion @@ Flatten[ReIm[First[MatrixFunction[Log,
{{q2[[1]] + I q2[[2]], q2[[3]] + I q2[[4]]},
{-q2[[3]] + I q2[[4]], q2[[1]] - I q2[[2]]}}]]]]
Quaternion[0.670588, -0.989246, 0.586497, -0.718306]

• Very interesting approach, thanks a lot =) May 5, 2020 at 14:49

I have decided to write another answer, because the following is more specific to the case of the gamma function. In particular, I will demonstrate how to modify the OP's original approach.

The reason the OP's method failed is that the elementary functions in Mathematica are left unevaluated for symbolic quaternion arguments. For instance,

<<Quaternions

Exp[Quaternion[a, b, c, d]]
E^Quaternion[a, b, c, d]


To make the OP's proposal work, one will need to manually split the quaternion into its scalar and vector parts. Using e.g. the formula here, we can write the following:

ig[Quaternion[a_, b_, c_, d_], x_] :=
With[{l = Log[x], vs = Sqrt[{b, c, d}.{b, c, d}]},
(x^(a - 1) Exp[-x]) Prepend[l {b, c, d} Sinc[l vs], Cos[l vs]]]


Then, we can do this:

qq = Quaternion[1, 2, 3, 0];

Quaternion @@ Integrate[ig[qq, x], {x, 0, Infinity}]
Quaternion[1/2 (Gamma[1 - I Sqrt[13]] + Gamma[1 + I Sqrt[13]]),
(I (Gamma[1 - I Sqrt[13]] - Gamma[1 + I Sqrt[13]]))/Sqrt[13],
(3 I (Gamma[1 - I Sqrt[13]] - Gamma[1 + I Sqrt[13]]))/(2 Sqrt[13]), 0]

N[%] // Chop
Quaternion[-0.003441450173902947, 0.008960473460683197, 0.013440710191024793, 0]


Compare with the matrix function approach:

Quaternion @@ Flatten[ReIm[First[MatrixFunction[Gamma,
{{qq[[1]] + I qq[[2]], qq[[3]] + I qq[[4]]},
{-qq[[3]] + I qq[[4]], qq[[1]] - I qq[[2]]}}]]]] // FullSimplify
Quaternion[1/2 (Gamma[1 - I Sqrt[13]] + Gamma[1 + I Sqrt[13]]),
(I (Gamma[1 - I Sqrt[13]] - Gamma[1 + I Sqrt[13]]))/Sqrt[13],
(3 I (Gamma[1 - I Sqrt[13]] - Gamma[1 + I Sqrt[13]]))/(2 Sqrt[13]), 0]


and we see the two methods give consistent results.

• The approach with the matrices I've found exciting, since I can adapt the calculation to work with reals, complexes or quaternions. May 6, 2020 at 14:28

With my GeometricAlgebra paclet you can work with any geometric number and apply any MatrixFunction supported function on them with MultivectorFunction , that effectively just converts numbers to matrices and back. Quaternions can be simply expressed as GeometricAlgebra[0,2]:

PacletInstall["https://wolfr.am/OkONsyY2"]

<< GeometricAlgebra

MultivectorFunction[Gamma,
Multivector[{a, b, c, d}, {0, 2}]] // Simplify


Convert quaternion to matrix:

MultivectorMatrix[Multivector[{a, b, c, d}, {0, 2}]]["Numeric"]