(* Title:
* The Barnes G-function
*
* Author:
* Victor S. Adamchik
* adamchik@cs.cmu.edu
*
* Summary:
* This package provides the numeric and symbolic computation of
*
* 1) the Barnes G-function;
* 2) Glaisher’s constant;
* 3) Clausen function;
* 4) derivatives of the Hurwitz and Riemann Zeta functions.
*
* Source:
* E. W. Barnes, “Genesis of the double gamma function”,
* Proc. London Math. Soc., 31(1900), 358-381.
*
* E. W. Barnes, “The theory of the double gamma function”,
* Philos. Trans. Roy. Soc. London, Ser. A, 196(1901), 265-388.
*
* V. Adamchik, “On the Barnes function”, Proceedings of the 2001
* International Symposium on Symbolic and Algebraic Computation
* (July 22-25, 2001, London, Canada), Academic Press, 2001, pp.15-20.
*
* Keywords:
* Barnes function, Gamma function, Riemann zeta function, Hurwitz
* zeta function, Stirling numbers, harmonic numbers, Catalan’s constant,
* Glaisher’s constant.
*
* History:
* Written by Victor Adamchik, June-October 2002.
*
* Discussion:
*
*
*
*)
BeginPackage[“BarnesG`”]
Unprotect[BarnesG]
BarnesG::usage =
“BarnesG[z] is the Barnes function. It can be viewed as a generalization of
the classical Euler gamma function. The Barnes function satisfies the
following functional relation BarnesG[z+1] = Gamma[z]*BarnesG[z].”
Begin[“BarnesG`Private`”]
(* ====================================================================== *)
(* SYMBOLIC COMPUTATIONS *)
BarnesG[z_/;Precision[z] == Infinity] :=
Module[{tmp},
tmp = BarnesSymb[z];
tmp /; NumericQ[tmp]
]/; NumberQ[z]
BarnesSymb[n_Integer?NonPositive] := 0
BarnesSymb[1] = 1
BarnesSymb[2] = 1
BarnesSymb[3] = 1
BarnesSymb[n_Integer?Positive] := Product[k!,{k,1,n-2}] /; n <= 10^4
(* closed form representations for z = 1/2, 1/4, 3/4, 5/4, 1/6, 5/6, 1/3, 2/3 *)
BarnesSymb[1/2] = 2^(1/24)*E^(1/8)/(Glaisher^(3/2)*Pi^(1/4))
BarnesSymb[1/4] = E^(3/32 - Catalan/(4*Pi))/(Glaisher^(9/8)*Gamma[1/4]^(3/4))
BarnesSymb[3/4] = E^(3/32 + Catalan/(4*Pi))*Gamma[1/4]^(1/4)/(2^(1/8)*Glaisher^(9/8)
*Pi^(1/4))
BarnesSymb[5/4] = (2^(1/8)*Pi^(1/4)*Gamma[3/4])/E^(Catalan/(2*Pi))
BarnesSymb[1/3] = (3^(1/72)*E^((12 + 2*Sqrt[3]*Pi - (3*Sqrt[3]*
PolyGamma[1,1/3])/Pi)/108))/(Glaisher^(4/3)*Gamma[1/3]^(2/3))
BarnesSymb[2/3] = (3^(1/72)*E^((12 + 2*Sqrt[3]*Pi - (3*Sqrt[3]*
PolyGamma[1,2/3])/Pi)/108))/(Glaisher^(4/3)*Gamma[2/3]^(1/3))
BarnesSymb[1/6] = E^((25 + 6*Sqrt[3]*Pi - (3*Sqrt[3]*
PolyGamma[1, 1/6])/Pi)/360)/(2^(1/72)*3^(1/144)*Glaisher^(5/6)*Gamma[1/6]^(5/6))
BarnesSymb[5/6] = E^(5/72 + Pi/(20*Sqrt[3]) - PolyGamma[1, 5/6]/(
40*Sqrt[3]*Pi))/(2^(1/72)*3^(1/144)*Glaisher^(5/6)*Gamma[5/6]^(1/6))
(* shift *)
BarnesSymb[z_?Positive] :=
With[{n = Floor[z]},
Product[Gamma[z+k-n],{k,0,n-1}]*BarnesSymb[z-n]]/;MemberQ[{2,3,4,6},
Denominator[z]]
BarnesSymb[z_?Negative] :=
With[{n = Floor[-z]},
BarnesSymb[z+n+1]/Product[Gamma[z-k+n],{k,0,n}]]/;MemberQ[{2,3,4,6},
Denominator[-
z]]
(* ====================================================================== *)
(* NUMERIC COMPUTATIONS *)
BarnesG[z_/;Precision[z] < Infinity] :=
If[MachineNumberQ[z], BarnesMachineNum[z], BarnesNumeric[z, Precision[z]]]
(* MACHINE PRECISION *)
BarnesMachineNum[z_] := If[Arg[z] == Pi, G1[z], G2[z]]
G1[z_] := (-1)^(Floor[-z/2] - 1)*G2[2 - z](Abs[Sin[Pi z]]/Pi)^(1-z) *
Exp[ClausenCl[2 Pi (-z - Floor[-z])]/(2 Pi)]
G2[z_] := Exp[NIntegrate[x*PolyGamma[x], {x, 0, 1/10, z - 1}] +
(z - 1)((2 - z) + Log[2*Pi])/2]
(* ARBITRARY PRECISION *)
Memoization = {} (* cacheing the results of previous computations *)
BarnesNumeric[z_, prec_] := BarnesNumeric[SetPrecision[z, $MachinePrecision],
prec] =
Module[{zz = SetPrecision[z, $MachinePrecision], max = 0},
If[MemberQ[Memoization, {zz, _}]
,
max = Max[Cases[Memoization, {zz,_}]/.{_,p_}->p];
,
AppendTo[Memoization, {zz, prec}]
];
If[prec < max
,
N[BarnesNumeric[zz, max], prec]
,
If[Arg[z] == Pi, G1[z, prec], G2[z, prec]]
]
]
G1[z_, prec_] := (-1)^(Floor[-z/2] - 1)*G2[2 - z, prec]*
(Abs[Sin[Pi z]]/Pi)^(1-z) * Exp[ClausenCl[2 Pi (-z - Floor[-z])]/(2 Pi)]
G2[z_, prec_] :=
With[{n = Floor[Re[z]] - 1},
Product[Gamma[z+k-n],{k,0,n-1}]*G2[z-n, prec]
]/; Re[z] > 2
G2[z_, prec_] :=
With[{n = Floor[-Re[z]] + 1},
G2[z+n+1, prec]/Product[Gamma[z-k+n],{k,0,n}]
]/; Re[z] < 1
G2[z_, prec_] :=
Exp[(z - 1)((2 - z) + Log[2*Pi])/2 +
NIntegrate[x*PolyGamma[x], {x, 0, z - 1},
PrecisionGoal -> prec,
WorkingPrecision -> prec+10]]
(* ======================== Clausen function =========================== *)
ClausenCl[x_] := – Im[PolyLog[2, Exp[-I x]]];
(* ======================== Glaisher =================================== *)
(* Glaisher’s constant is defined in MMA V4.2 *)
(* This is to support earlier versions *)
If[ !NumericQ[Glaisher]
,
Glaisher::usage =
“Glaisher is Glaisher’s constant defined by Log[Glaisher] = 1/12 –
Zeta'[-1].”;
Attributes[Glaisher] = {Constant, Protected, ReadProtected};
N[Glaisher] := N[Glaisher500];
Glaisher500 =
1.2824271291006226368753425688697917277676889273250011920637400217\
40406308858826461129736491958202374394206461203990007489331577913627752804041
5\
90725738617275221433432714343978733506791525736685690787656114668644999778496
2\
75451817431239465276128213808180219264516851546143919901083573730703504903888
12\
34188136749781330509377083368222249411587483734806439997883007012556700128699
41\
57705432053927585405817315881554817629703847432504677751473746000316160230466
13\
296342991558095879293363438872887019889535;
]
(* Numeric computation.
*
* This algorithm is slightly faster than computing Zeta'[-1]
*
*)
Unprotect[ Glaisher ];
N[Glaisher, p_?NumberQ] :=
If[p <= 500, N[System`MathConstantsDump`Glaisher500, p],
With[{top = Ceiling[p Log[2, 10] /2]},
2^(1/12)*Exp[Sum[(Zeta[2k+N[1, p*1.07]] - 1) *
(28+3/(k+1)-6/(k+2)), {k, 1, top}]/36]
]
]
Protect[ Glaisher ];
(* ======================== Zeta'[-1] ================================== *)
N[Derivative[1][Zeta][-1]] := N[1/12 - Log[Glaisher]];
N[Derivative[1][Zeta][-1], prec_] := 1/12 - Log[N[Glaisher, prec]];
(* ======================== Zeta'[-1, z] =============================== *)
(* Derivatives of the Hurwitz function by means of
* the Barnes function.
*
*)
N[Derivative[1, 0][Zeta][-1, z_/; Re[z] > 0]] :=
N[1/12-Log[Glaisher]-Log[BarnesG[z+1]]+z Log[Gamma[z]]]
N[Derivative[1, 0][Zeta][-1, z_/; Re[z] > 0],prec_] :=
N[1/12-Log[N[Glaisher,prec]]-Log[N[BarnesG[z+1],prec]]+z Log[Gamma[z]], prec]
(* ====================================================================== *)
End[]
Attributes[BarnesG] = {Listable, NumericFunction, Protected, ReadProtected};
Protect[BarnesG]
EndPackage[]