Wolfram Function Repository

Instant-use add-on functions for the Wolfram Language

Function Repository Resource:

 BarnesMultipleGamma

Source Notebook

Compute the Barnes multiple gamma function numerically

Contributed by: Pieter-Jan De Smet

ResourceFunction["BarnesMultipleGamma"][w,{a1,...,an}]

gives Γn(w❘a1,…,an).

Details

The Barnes multiple gamma function is defined as
where the Barnes zeta function ζN is defined as
ζN(s,w❘a1,…,aN)=∑m1,…, mN=0+∞1/(w+m1 a1+⋯+mN aN)s
It is everywhere assumed that a1>0, …, aN>0.
The sum (2) converges for ℜ(w)>0 and ℜ(s) > N. It admits a meromorphic continuation in s that is holomorphic at s=0.
Barnes multiple gamma function ΓN admits a meromorphic continuation in w to the complex plane ℂ. All its poles lie on the interval (-∞,0].
The integer order of ResourceFunction["BarnesMultipleGamma"] is the length of the parameter list.
BarnesMultipleGamma computes a numerical approximation.
Barnes multiple gamma function satisfies a recurrence relation

Examples

Basic Examples (1) 

Evaluate numerically:

In[1]:=
w = -1.2 + 1.3 I; a1 = 1.2; a2 = 4.5; ResourceFunction["BarnesMultipleGamma"][w, {a1, a2}]
Out[4]=

Scope (1) 

BarnesMultipleGamma supports arbitrary integer orders subject to limits in computing accuracy:

In[5]:=
Table[ResourceFunction["BarnesMultipleGamma"][2.5 + 1.4 I, ConstantArray[1, n]], {n, 0, 5}]
Out[5]=

Applications (1) 

The Barnes multiple gamma function can be used to define the double sine function:

In[6]:=
doubleSine[w_, a_List] := ResourceFunction["BarnesMultipleGamma"][Total[a] - w, a]/ ResourceFunction["BarnesMultipleGamma"][w, a] doubleSine[.7 + .3 I, {1, 2}]
Out[7]=

Properties and Relations (4) 

If the order is zero then Γ0(w) = 1/w:

In[8]:=
w = 1.2; {ResourceFunction["BarnesMultipleGamma"][w, {}], 1/w}
Out[9]=

If the order is one then Γ1 is related to the usual gamma function Γ1(w❘a)=ⅇ(w/a - 1/2) log aΓ(w/a)(2 π)-1/2:

In[10]:=
w = -4.07 - 1.41 I; a = 1.2; {ResourceFunction["BarnesMultipleGamma"][w, {a}], Exp[(w/a - 1/2) Log[a]] Gamma[w/a] (2 Pi)^(-1/2)}
Out[12]=

If the order is two and a1=a2=1 then Γ2 is related to the Barnes G function, Γ2(w❘1,1)=ⅇ(w/2 - 1/2) log (2 π)eζ′(-1)G(w)-1:

In[13]:=
w = 4.07 - 1.41 I; {ResourceFunction["BarnesMultipleGamma"][w, {1, 1}], Exp[(w/2 - 1/2) Log[2 Pi]] Exp[Zeta'[-1]] BarnesG[w]^(-1)}
Out[14]=

BarnesMultipleGamma satisfies a recurrence relation:

In[15]:=
w = 2.5 + 3 I; a1 = 1.2; a2 = 1.5; {ResourceFunction["BarnesMultipleGamma"][w + a2, {a1, a2}], ResourceFunction["BarnesMultipleGamma"][w, {a1, a2}]/ ResourceFunction["BarnesMultipleGamma"][w, {a1}]}
Out[18]=

Possible Issues (3) 

The parameters should be strictly positive real numbers:

In[19]:=
ResourceFunction["BarnesMultipleGamma"][.4, {1, -2}]
Out[19]=

The recurrence method in use can lead to exceeding recursion depth limitations:

In[20]:=
ResourceFunction["BarnesMultipleGamma"][1 + 10^5 I, {1, 12}]
Out[20]=

BarnesMultipleGamma does not compute exact values even when input is exact:

In[21]:=
ResourceFunction["BarnesMultipleGamma"][2, {1, 2, 3}]
Out[21]=

Neat Examples (1) 

For special values of the parameters one can often calculate closed form expressions involving the Hurwitz zeta function:

In[22]:=
BarnesGamma123[w_] := Module[{b = {3, 4, 5, 6, 7, 8}, c = {1, 1, 2, 3, 4, 5}, q, B, C, zetaSum}, zetaSum[s_] := Sum[q = (w + r)/6; B = b[[r + 1]] - 6 q; C = c[[r + 1]] - b[[r + 1]] q + 3 q^2; 3 HurwitzZeta[s - 2, q] + B HurwitzZeta[s - 1, q] + C HurwitzZeta[s, q], {r, 0, 5}]; Exp[(D[zetaSum[s], s] /. s -> 0) - Log[6] (zetaSum[0])]]; w = 2.5 + 1.4 I; {ResourceFunction["BarnesMultipleGamma"][w, {1, 2, 3}], BarnesGamma123[w]}
Out[24]=

Publisher

Pieter DeSmet

Requirements

Wolfram Language 13.0 (December 2021) or above

Version History

  • 1.0.0 – 23 October 2026

Source Metadata

  • Citation:
    • [R] Ruijsenaars, "On Barnes’ Multiple Zeta and Gamma Functions'', Advances in Mathematics 156 (2000), pp. 107-132 (see equation (1.5) for double sine definition and equation (3.27) for the order one case)

Related Resources

Author Notes

Mathematical background

The numerical evaluation combines a contour-integral representation with the recurrence relation (3). The recurrence is used to shift the argument w to the sector ℜe(w) >ℑm(z).

The numerical integration is based on the following formula:

The contour 𝒞 in (4) is a "Hankel contour" that winds from -∞ - 0 ⅈ in the lower half-plane, around 0, and back to -∞ + 0 ⅈ in the upper half-plane. Formula (4) holds for s ≠ 0, 1, 2,… and ℜ(a)>0.

Applying (4) to the defining series (2) gives (5):

Performing the geometric sums yields (6):

Using definition (1) in equation (6) gives

Here γ is the Euler–Mascheroni constant. The integral converges for ℜ(w)>0. This is the main formula used for numerical calculation.

For t with large negative real part, the integrand in (7) behaves roughly as ⅇwt/t, ignoring the logarithm. Analyzing the steepest descent contours leads to a rotated Talbot contour around a half-line that points in the direction. So, the contour z(θ) has the form

We take R sufficiently small so the Hankel contour can be deformed into the rotated Talbot contour without crossing poles. We also choose the direction of the branch cut of the logarithm in the direction .

Illustration

Here is an illustration for w = 2 + ⅈ, N=3, {a1,a2,a3} = {1,2,3}

In[1]:=
Clear[\[Theta], t, f, a, ai, asorted, an, arest] w = 2 + I; a = {1., 2., 3.}; (* Log[t/Conjugate[w]] - Log[1/ Conjugate[w]] is the logarithm of t, with branch cut in the direction -Overscript[w, _] *) f[t_?NumericQ] := 1/t ( EulerGamma + Log[t/Conjugate[w]] - Log[1/ Conjugate[w]]) Exp[ w t] 1 / Product[1 - Exp[ai t], {ai, a}]; mySincPrime[\[Theta]_?NumericQ] := -SphericalBesselJ[1, \[Theta]]; \[Phi] = Arg[Conjugate[w]]; (* R should be sufficiently small, it is constrained by the nearest nonzero pole *) asorted = Sort[a, Less]; an = asorted[[-1]]; arest = asorted[[1 ;; -2]]; list = {Pi/an Abs[Sin[Pi/2 - \[Phi]]/(Pi/2 - \[Phi])], Pi/an Abs[Sin[-Pi/2 - \[Phi]]/(-Pi/2 - \[Phi])], 1 / Abs[w]}; R = Min[list]; z[\[Theta]_] = R / Sinc[\[Theta]] Exp[I (\[Phi] + \[Theta])]; dz[\[Theta]_] = R Exp[I (\[Phi] + \[Theta])] ( I /Sinc[\[Theta]] - mySincPrime[\[Theta]]/Sinc[\[Theta]]^2);
In[2]:=
\[Theta]min = -2.7; \[Theta]max = 2.7; figTalbot = ParametricPlot[{Re[z[\[Theta]]], Im[z[\[Theta]]]}, {\[Theta], \[Theta]min, \[Theta]max}, PlotStyle -> Green];
In[3]:=
poles = Table[{0, 2 Pi n /ai}, {n, 0, 1}, {ai, a}]; poles = Flatten[poles, 1];
In[4]:=
xmin = -4; xmax = 1; ymin = -2; ymax = 4; poles = Select[poles, #[[2]] < ymax &]; figpoles = ListPlot[poles, PlotStyle -> {PointSize[Medium], Red}];
In[5]:=
figf = ComplexPlot[f[t], {t, xmin + I ymin, xmax + I ymax}, MeshStyle -> {Directive[Thin, White], Black}, Mesh -> 20, MeshFunctions -> {Log[Abs[#2]] &, Arg[#2] &}, RasterSize -> 30];

In the following figure, one sees the rotated Talbot contour (green) and selected poles of the integrand (red) in the complex t-plane. The black lines are the lines of constant phase. The colours indicate the phase of the function f. One sees that the Talbot contour more or less follows the lines of constant phase.

In[6]:=
fig = Show[figf, figTalbot, figpoles]; (* Rasterize to reduce the size of the file *) Rasterize[fig, "Image", ImageSize -> 200]

The following figure shows the integrand. Although there are oscillations, they are sufficiently mild that NIntegrate can deal with them.

In[7]:=
fig2 = ReImPlot[dz[\[Theta]] f[z[\[Theta]]], {\[Theta], -3, 3}, PlotRange -> All]; (* Rasterize to reduce the size of the file *) Rasterize[fig2, "Image", ImageSize -> 200]

License Information