An algorithm for evaluating the Gamma function and ramifications
11 An algorithm for the evaluation of the Gamma functionand ramifications
D. KARAYANNAKIS
T.E.I. of CreteDept. of Science / Applied Math. [email protected]
A closed type formula is provided for Γ(x), when x is rational in (0,1), based on a recursive algorithm which, along with its equivalent formulation as an infinite product, furnishes easy to obtain estimates and leads to an additional insight on various properties concerning Γ(x), 1/Γ(x), lnΓ(x), B(x,y), ψ(x) and ψ’(x).The method allo ws the production of estimates of non- classical values of Γ(x) and a variety of pertinent facts without resorting to any kind of specialized programming tools beyond a scientific pocket calculator.Formulae some of them seemingly new, involving infinite products and inequalities, are presented as well.Keywords-Gamma function, Infinite product, Psi function
1. INTRODUCTION
The theory of the Gamma function, Γ(x), was historically developed in connection with the problem of extending the value of n! for positive integers n to arbitrary real numbers.The most fundamental functional equation, that can be derived almost immediately from its definition and which is basic for the development of the theory of the Gamma function is:
Γ(x+1) = x Γ(x). (1.1)It represents a generalization of the identity n!=n(n-1)! for non integral values of n. Note that Eq.(1.1) allows us to extend the original classical definition to include negative real numbers, provided we avoid negative integral values and zero. In the rest of the paper we will always assume that knowledge of the values of the Gamma function on the interval (0, 1] results to knowing the values of Γ(x) on \Z .Legendre’s doubling formula states that ([4]): Γ(2x) π Γ(x+1/2)Γ(x). (1.2)
When (1.2) is employed for x=1/2 it gives immediately Γ(1/2)=π , while for x=n can lead to the evaluation of Γ(n+1/2) in closed form ( cf. [6]); but these are basically the most this formula can do to that direction.Eq. (1.3) is a special case of Gauss’ multiplication formula ([4]):
Γ(nx) ) π
2( n Γ(x) Γ(x+1/n)… Γ(x+(n -1)/n). (1.3)Eq. (1.3) in most of the cases leads to identities that at first glance look rather exotic (cf. [6]) and its use is subjected to the same, more or less, limitations described above for Eq. (1.3); nevertheless it will be instrumental for our investigations since it will provide the constant that multiplies the infinite product representation of Γ(1/p), for any integer p, p
2. A CLOSED TYPE FORMULA FOR Γ(x+b)/Γ(x), x>0, 0 b F(α, β, γ; z) = 1 + γ βα z + γ ( γ )1 β ( β )1 α ( α z + (2.1) A hypergeometric seri es terminates if α or β is a negative integer or zero. If we exclude the case γ = -n, n=0,1,2,… the series (2.1) converges in the unit circle |z|<1. The sum-function, the hypergeometric function F, (where for simplicitly we prefer this symbolism instead of the usual 2F1) has a branch point at z=1. Still, whenever Re(α+β - γ)<0, the series converges absolutely throu ghout the entire unit circle ([10]). From now on we will consider only the case of real parameters α, β, γ and real z= t and we will be concerned mainly with the special case t=1. We have the following formula known as Gauss Hypergeometric (or Summation) Theorem ([10]):
F(α, β, γ; 1) βγΓαγΓ βαγΓγΓ (2.2) under the condition γ>α+β; (2.2) is valid also for complex va riables whose real parts obey the analogous condition.It is also known ([2]) that the (real) hypergeometric series is one of the solutions of the differential equation:t(1- t) y΄΄(t) + [γ – (α+β+1)t] y΄(t) – αβy(t) = 0. (2.3)We recall that for the function psi, i.e. the logarithmic derivative of Γ(x), we have the formula ([6]): ψ(u) – ψ(t) = ku 1kt 1 (2.4)whenever u and t are nonnegative integers or zero.We will first present a “formal” (heuristic) proof of the following:LEMMA 2.1. F (α, b, 1; t) =
1n nn tb, α gexp for α+b<1, 0 t
1, where g (α,b)=αb, g (α,b)=αb(1+α+b - αb)/4, and the g n ’s (for n
2) are defined through the recursive algorithm:(n+1) g n+1 = n(n+α+b -g )g n + kn1kn1k2n 0k gkng1kng1k . (2.5) PROOF. The differentia l equation (2.3) for γ=1 leads, after the substitution: F (α, b, 1; t) = C exp g(t)to the equation: t(1-t)[(g’) + g’’] + [1- (α+b+1)t]g’ – αb = 0. (2.6)Solving (2.6)formally by the power-seriesmethod ([8]) we obtain g(t)=g + g n (α,b) t n , where g n (α,b) can be retrieved via (2.5) for n
2. But, since: F (α, b, 1; 0)=1 we can conclude immediately that kexpg =1, which leads to the announced result.REMARK 2.1. From now on we will consider 1- α -b=x > n will be considered functions of x and b, and, in particular, g =b(1-b-x) and . Observe also that the coefficient of g n in (1.5) takes now the form n[n+1-x-b(1-b-x)]. Combining now Lemma 2.1 and (2.2) we can conclude immediately : PROPOSITION 2.1. Γ(x+b) = f(x,b) )b1( Γ )x( Γ , for x > 0, 0 b < 1, where f(x,b)=
1n n b,xgexp . COROLLARY 2.1. f(x,b) sin(πb) = f(b,x) sin(πx), 0 < x, b < 1.
PROOF. Using the formula of Proposition 2.1 we automatically see that for 0 Γ(x+b)=f(x,b) )b1( Γ )x( Γ = Γ(b+x) = f(b,x) )x1( Γ )b( Γ f(x,b) Γ(x) Γ(1 - x) = f(b,x) Γ(b) Γ(1 -b). But since Γ(z) Γ(1 -z) )z π sin( π (cf. [4]) for 0 1n n = kx 1kbx 1 . Integrating now for x we obtain directly that: g n (x,b)=C + ln kx kbx . Observing now that g (x,1-b)=g (x,1-b)=0, the recursive formula (2.5) leads to g n (x,1-b)=0 for all n. We conclude that C = _ ln kb1 k1 and thus g n (x,b)= ln )1k)(kx( )b1k)(1bkx( , which is equivalent to the announced result.REMARK 2.2.(i) In the above proof we have interchanged heuristically the order of summation with differentiation and integration; still this interchange is in fact legitimate because the polynomial series (in x) involved, converges uniformly on any closed subinterval of (0, + ), due to the very own construction of the g n ’s via the power series method. In addition it is evident that Γ ( x+b )/ Γ (x) is real analytic (in x) on (0 , ). (ii) It is also evident that the infinite product representation of f(1-b,b) should be considered equal to 1, a convention compatible with the statement of Proposition 2.1 and due to the obvious fact that g n =0 for all n. (iii) For notational convenience we will set the product of the first m terms in the infinite product representation of f(x,b) as f m (x,b) which we will call the m-truncate of the joint factor.Using the functional equation of Proposition 2.1 for b=1/2 we obtain:COROLLARY 2.2. (Legendre’s formula revisited) Γ(2x) π f(x,1/2)[Γ(x)] , x > satisfy the functional equation B(x,y)=Γ(x)Γ(y)/Γ(x+y) for x>0, y>0; in the open strip (0 , ) x (0 , π y), we can obtain immediately:COROLLARY 2.3.B(x,y) = (1/y) .)x1k)(yk(/)yx1k(k REMARK 2.3. Corollary 2.2 provides us with a closed-type formula that not only recovers, but also simplifi es many infinite product representations of [Γ(x)] , scattered throughout the pertinent literature ([1]). For example, it provides us with a representation of [Γ(1/4)] “simpler” than the classical one (cf. [5]), which states (after squaring off) that π Γ . The aforementioned new representation reads π π Γ , and it is a direct application of (2.8). The classical formula has a better “starting point” and, as simple algebra can demonstrate, it always “stays ahead” of its new counterpart. Still one gets compensated by the mere fact that the latter involves exclusively rational fractions and that for very large k’ s we achieve numerical approximations of the same order.REMARK 2.4. Based on the functional equation Γ(1 -x)= - xΓ( -x) , which is valid for non integral x>0, we now see that Corollary 2.2 can easily be extended to negative arguments in the following form: Γ (-2x) = -[ 2 -2x cosec(2 π x)/xf(x,1/2) ][ π / Γ (x)] for 0 3. A CLOSED- TYPE FORMULA FOR Γ(x), 1/Γ(x) AND ln Γ(x), x RATIONAL For evident reasons, we will restrict ourselves to (0,1). We first prove a lemma which, by use of the joint factors introduced in Proposition 2.1, leads to the main result of this chapter:LEMMA 3.1. p p/1 Γ p) π 2p 1k )p/k,p/1(f 1 , for p any integer 3. (3.1)PROOF. By use of Gauss’ multiplication formula and setting n=p and x=1/p in (1.4), we obtain in succession, via our functional equation of Proposition 2.1, that p/1 Γ = p/1 Γ p/2 Γ = f(1/p,1/p) p/11 Γ p/1 Γ p/3 Γ = f(1/p,2/p) p/21 Γ p/1 Γ …….. p/)1p( Γ = f(1/p,(p-2)/p) p/)2p(1 Γ p/1 Γ .Multiplying termwise the above equations and since: 1p 1k Γ(k/p) p) π we obtain the announced result. The functional equation (1.2) and the above lemma gives now directly the following:COROLLARY 3.1. Γ (-1/p) p = -1/sin( π /p) p π /(2p) p-1 2p 1k f(1/p,k/p). We are ready now to prove the subsequentTHEOREM 3.1 Γ(q/p) = C p,q 1q 1k f(k/p, 1/p) 2p 1k p/q )p/k,p/1(f 1 (3.2)where C p,q p/q 1q p π π sin2 π for p, q integers and 1 1q 1k p/1,p/kf and thus, by use of (3.1), we arrive at the announced result.One can now automatically obtain the following results:COROLLARY 3.2.(i) p/q Γ q,p C1 2p 1k [f(1/p,k/p)] q/p 1q 1k p/1,p/kf 1 (ii) lnΓ(q/p) = (lnC p,q )+ 1q 1k μ k (p) – (q/p) 2p 1k ν k (p), where μ k (p) = lnf(k/p, 1/p), 1 k q-1 and ν k (p) = lnf(1/p, k/p), 1 k p-2.REMARK 3.1(a) It is now evident that, via the classical continuity argument, (3.2) provides in principle the mechanism for the evaluation of Γ(x) at any number x in (0,1) and thus at any real x not in . But the actual numerical implementation of this remark would demand, aside to a sophisticated programming, the aid of numerical analysis methods and the pertinent approximation theory and thus lies beyond the scope of the present work.(b) Theorem 3.1 provides an alternative way to reproduce a rather big part of the explcit formulae presented in Section 2 of [9] (pp 269-270).4. APPLICATIONSIn this chapter we demonstrate the method that allows us, by use of the algorithm (2.5) as modified in Remark 2.1, to produce interesting representations of non classical values of Γ(x) and other transcendental functions, as well as some presumably new inequalities and features concerning Γ(x) and its applications : APPLICATION 1.It is rather evident that the joint factor f(1/p,k/p) has f m (1/p,k/p) as a strict lower bound for all positive integers p, k and m such that p>2 and k ≥p -2; in particular, for m=1, Lemma 3.1 and Corrolary 3.1 combined lead now to the following pair of inequalities:(i) Γ (1/p) < (2 π /p) (p-1)! (ιι) Γ( -1/p) - π /sin( π /p) (2 π ) Note that, numerically, (i) and (ii) provide rather crude bounds because of the particularly small value of m; the larger the value of m the more refined inequalities will be produced. We illustrate this point extensively in the reasoning of Applications 9 & 10. APPLICATION 2.For x>0, using Proposition 2.1 twice and the resulting product of the two joint factors as a single infinite product,we see that x Γ(x) = Γ(x+1) = Γ(x+1/2) f(x+1/2,1/2) / π = Γ(x)f(x,1/2)f(x+1/2,1/2)/π which results to x π xk/)x1k()1k2/k2( 1k 2 ; the latter when combined with the infinite product representation of (sinπx) /πx produces, another representation of sin(πx), namely sin(πx) = .)x1k/()xk()1k2/2( 1k 2 .APPLICATION 3.Setting x=b, 0 < b < 1, in Proposition 2.1 and using Corrolary 2.2 one gets immediately that sin( π b)f(b,b)/f(b,1/2) = 2 ; then employing the joint products involved, we obtain the subsequent formula that, modulo the formulae that one could randomly encounter when surfing on the web, we were not able to trace in the pertinent literature: (n-1/2)(n-1/2 + b)/(n-b)(n-1+2b)= 2 /( sin π b ). Note that in particular for b=1/4 and b=1/2 we set by convention the corresponding infinite product of 1’s equal to 1. APPLICATION 4.By Proposition 2.1 for b=1/2 we obtain that x Γπ Γ ; we also have that x21 Γπ Γ x1 Γ for 0 1n 22 )x2(1n2 )x2n2)(x22n2( . The corresponding classical infinite product representation -directly from the classical math literature, or indirectly when one combines the formula stated for sin(πx) with the corresponding one for cos(πx), (e.g. π x) = πx x21n2/xnn/1n2 .APPLICATION 5.Let us combine now the basic fact of the strict convexity of Γ (x) over (0,+ ) with the factorization presented in Proposition 2.1 in order to produce a strict upper bound for Γ ( α ), 0< α < 1, provided α stays away from zero. We know that the classical inequality Γ ( α x+(1- α )y)< αΓ (x)+(1- α ) Γ (y), is true for x>0, y>0, x y; if we set x=1+y and y=1- α + ε with 0< ε < α <1, then using the fact that the joint factor f(1- α , ε ) is strictly bounded from above by any of its m-truncation (compare to the similar reasoning of Application 2), we obtain: Γ ( α ) < [sin( πε ) / sin( πα )][(1+ εα α , ε ) / ε] for all m є N. In particular, when we set m=1, we have the following inequality: Γ ( α ) < [sin( πε ) / sin( πα )][(1- α )(1+ εα - α ) / ε (1- ε )(1+ ε – α ) ]; letting ε=α we have the indicative approximation Γ(α) < 1/α on (0,1), which evidently performs quite satisfactory whwnever α is close to the endpoints. It is also clear that for fixed ε and α (ε≠α), the larger the m the better the m - ub for Γ(α). APPLICATION 6.It is immediate now that Cor.2.4 leads to the following inequalities:(i) If 0 1 and any m є N, I α < [2 m (m!)(α/2) m ] / [(2m- m ], where (.) m denotes Pochammer’s symbol ([6]).(ii) For α>1 the inequality in (i) is reversed. Note that even in the case of the crudest possible estimates, obtained for m=1, the above inequalities provide a better l.b. and u.b., respectively, of Wallis’s integral in the corresponding interval than the ones provided by the classical trigonometric inequalities. For example, when 0<θ<π/2, we know that 2θ/π < sinθ < θ and thus we conclude that π /2(α+1)< I α < (π/2) α+1 / (α+1). On the other hand, for 0<α<1, (i) states that: I α < 2 /(α+1) < (π/2) α+1 /(α+1),while for α>1 (ii) states that: I α > 2 /( α +1) > π /2( α +1). place accuracy; on the other hand the “old classical formula“ for the same upper index produces exactly the same accuracy.APPLICATION 9.In ([2]) H. Alzer, among other sharp inequality conditions about the Gamma function, obtained the following improvement of Gamma inequalities discovered in 1997 by Anderson and Qiu: Let α=1 - γ and β=(π - (i) If 0 1 and 0.37 > 3(β+γ -1). Q.E.D. Proof of (f): we check that s(1.562) < ε < s(1.561) and that s(x) is strictly decreasing on (1, ) which is true since s(x)’ = s(x) [ln(1+0.5/x) – γ / (x+0.5) – γ / x] < s(x) [(0.5 – γ) / (x + 0.5) – γ / x] < 0. Q.E.DAPPLICATION 10.We conclude this chapter by examining a numerical impact of the joint factor approach in the revived research for improvement of Stirling’s formula; we will make use of Nörlund’s asymptotic representation of the Gamma function ([7]) and of arelatively recent result due to W. Schuster ([8]). Stirling’s formula states that Γ(x) = (2π) x x-1/2 e -x e μ(x) , with 0 < μ(x) < 1/12x, while Nörlund’s more effective asymptotic relation reads Γ(x+1/2)=(2π) x x e -x e v(x) ,where the bounds of v(x) have been improved (see (17) in [7]) as follows:(-1/24)[(1/x) + (1/120x )] ≤ v(x) ≤ ( -1/24)[(1/x) – (1/8x )].Combining now the joint factor f(x,1/2) and Stirling’s formula in Nörlund’s relation, we easily obtain that v(x) = μ(x) + ln[f(x,1/2) / x π ] and thus for 0 5. THE DIGAMMA AND THE TRIGAMMA FUNCTION.We conclude this work by considering two more general applications of the concept of the joint product, as well as of algorithm (2.5) concerning the psi (ordigamma ) function ψ(z), and the t rigamma function ψ’(z) respectively. The strate gy to be presented can be repeated with the analogous results in case of any polygammafunction restricted to (0,1) and and it can lead to various rather tedious series whose formulation entails theoretical aspects also. It is an immediate consequence of the proof of Theorem 2.1 that: ψ(x+b)=ψ(x) - g’ n (x,b) for x>0, 0 b <1 and g’ n (x,b) b,xf b,xf (where we restrict b to (0,1) to avoid indeterminancy for the case b=0). Ι f we resort, after differentiating (2.5) w.r.t to x, to the resulting new recursive formula that involves the sequence b,xg n =h n (x,b), we will be able to obtain ψ(x+b) and also, whenever necessary, the form of f ’(x, b). This result besides being helpful in producing and/or reproducing identities as well as app roximating ψ(x) and ψ’(x) for x>0, seems to bear an interest of its own. We present this proposed algorithmic strategy by providing in detail an illuminating abstract example concerning the case x=1-b, 0
Let us recall that g (1-b,b)=g (1-b,b)=0. For n g n+1 (x, b) = = n[n+1-x-b(1-b-x)]g n (x, b)+ b,xgknb,xg1knb,xg1k kn1kn1k2n 0k and thus g n (1-b, b)=0 for all n. Preserving the x notation described in Remark 2.1 and after differentiating wrt x, we obtain at the point x=1-b(n+1) h n+1 = n(n+b)h n (5.1)for n 2, where h b and h b1b41 ; note that we opressed the dependence of h n from b for reasons of simplicity. Multiplying (5.1) for n 2 and after telescoping we arrive at h n b Γ bn Γ !nn 1 (5.2)for n 3. Observe also that for n=1,2 we also obtain the correct values for h n and thus (5.2) holds for all n. Since ψ(x+b) = ψ(x) - b,xh 1n n for all x>0, 0
1n 2 n t1,nf π t π sin γ t ψ (5.4)REMARK 5.1. A rather interesting improvement of the precision that can be obtained from (5.4) is feasible by use of the evident fact that f(n,1- t) = o(Γ(t)n ):Suppose that n is very large in the identity n t1,nfn t1,nfn t1,nf . Thus we see that n 1)t( Γ n t1,nfn t1,nf and so we conclude that: n1)t1( Γ 1n t1,nf π t π sin γ t ψ (5.5) REMARK 5.2. It is rather evident that for n very large (5.5) furnishes one more “interconnection” between the psi function and the Riemann’s zeta function when the former is restricted to (0,1), with the aid of the joint factors namely: )t1( Γ )t1( ζ n t1,nf π t π sin)t1( Γ γ t ψ n 1n 2 .In order to obtain a similar expression for the trigamma function, let us now differentiate (5.4) wrt t. We then have that: 1n 21n 2 n t1,nf π t π sinn t1,nft π cost ψ .Since we are dealing once more with the derivative of an infinite product of functions (in t) we eleviate this difficulty through the classical logarithmic differentiation, namely t1k t1,nft1,nf . By substitution in the above formula of t ψ we arrive at: 1n 2n1k n tk 1t1,nf π t π sint ψ (5.6)Finally let us construct another approximation formula that will improve the calculating efficiency of (5.6) working along the same lines as with (5.5); we can immediately conclude that, for n very large: n tk 1)t1( Γ 1n tk 1t1,nf π t π sint ψ (5.7)REFERENCES1. M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972.2. Η. Alzer, Inequalities for the Gamma Function, Proc. A.M.S (2000), v.128, 141 -147.3. G. Andrew’s, R. Askey and R. Roy, Special Functions, Cambridge Univ. Press, 1999.4.