Sum of powers from uniform distribution moments

Journal of Statistical Planning and Inference 101 (2002) 57–79

www.elsevier.com/locate/jspi

Sum of powers from uniform distribution moments Larry Ericksen P.O. Box 172, Millville, NJ 08332, USA Received 18 November 1998; received in revised form 16 August 1999

Abstract The sum of powers formulas for the 0rst (n) consecutive integers may be viewed as the moments of discrete uniform distributions. These sum of powers equations are extended here for arbitrary moment con0gurations of the symmetric coe3cients in a generalized Pascal’s triangle. Moment equations are also derived for alternating signed versions of these coe3cient moments. Moment generating functions for the sum of powers equations are presented in closed and in exponential sum formats. The moment equations and their exponential generating functions are expressed in terms of generalized Bernoulli, Euler and Chebyshev Polynomials, and by the Lerch transcendentals of the Hurwitz zeta and the Dirichlet eta functions. Following Faulhaber’s lead, moment equations are written as power series of the term number (n). Finally, sum of powers tables are listed to facilitate moment calculations and to highlight c 2002 Elsevier Science B.V. All rights coe3cient characteristics in the moment equations.
reserved. MSC: 05A15; 60E05 Keywords: Moments; Sum of powers; Uniform distributions; Pascal–DeMoivre coe3cients; Generating functions; Faulhaber; Bernoulli numbers; Euler polynomials; Lerch transcendentals; Chebyshev polynomials

1. The objective The well-known formulas (Comtet, 1970; Nielsen, 1923; Zwillinger, 1996) for the
n sum of powers equations are shown for the 0rst few R powers of h=0 hR as n


n h = (n + 1); 2 h=0

n


n h2 = (n + 1)(2n + 1); 6 h=0

n
h=0

h3 =

n2 (n + 1)2 : 4

(1)

This paper presents the sum of powers of integers shown by (2) in both positive and alternating signed formats. The 0rst n + 1 consecutive squares constitute a series such as (02 + 12 + 22 + 32 ) whose sum is 14. And a similar series of alternating signs E-mail address: [email protected] (L. Ericksen). c 2002 Elsevier Science B.V. All rights reserved. 0378-3758/02/$ - see front matter
PII: S 0 3 7 8 - 3 7 5 8 ( 0 1 ) 0 0 1 5 4 - 9

58

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

like (02 − 12 + 22 − 32 ) has a sum of −6. This report derives formulas that generate these summation results n n

(−1)h hR : (2) hR and h=0

h=0

The sum of powers equations are further expanded to include sequences with spacing d and initial term k, as described by Eq. (3). The series terms in this category may be considered to be numbers of k(mod d), such as the sum of the cubes of the 0rst three odd integers (13 + 33 + 53 ) whose sum is 153 and whose alternating signed counterpart totals 99. n n

(3) (dh + k)R and (−1)h (dh + k)R : h=0

h=0

The method used to evaluate these sum of powers equations is to treat them as the Rth moments for the Pascal–DeMoivre distributions of Ch (N; J ) coe3cients from Ericksen (1998). With the relationship between distribution variable J and term variable n 0xed at J = n + 1, the most general cases in this paper for these sum of powers equations are Nn
h=0

(dh + k)R Ch (N; J )

and

Nn
h=0

(−1)h (dh + k)R Ch (N; J ):

(4)

2. The outline In general, the sections of this report describe either term de0nitions, generating functions, or moment equations. The de0nitions for basic terminology and notation used in the formulas are covered in sections on Pascal–DeMoivre coe3cients, distribution moment types, and spaced coe3cient moments. The de0nitions for generalized polynomials are presented in sections for Bernoulli and Euler polynomials, and Lerch transcendentals. Formulas for the moment generating functions and the moment equations are derived in sections on exponential generating functions, moments by Bernoulli and Euler sums, Lerch transcendentals, and Faulhaber formulas. Related topics are covered in sections on inverse equations for integer powers and Chebyshev polynomials of the second kind. Where informative, references are given for known formulas for the sum of powers equations. Especially for the simple formulations of (1) and (2), there is much history and documentation on the sum of powers formulas, corresponding to the moment generators and moment equations in this paper. 3. Pascal–DeMoivre coecients The Pascal–DeMoivre coe3cients Ch (N; J ) are de0ned from the expansion of (1 + z + z 2 + · · · + z J −1 )N as the coe3cients of a power series in z h , for all positive values

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

59

Table 1 Pascal–DeMoivre coe3cients Ch (N; 3) N

h=0

1

2

3

4

5

6

1 2 3

1 1 1

1 2 3

1 3 6

2 7

1 6

3

1

of the distribution variables (N; J ). The generating function for this expansion and its Ch (N; J ) coe3cients is 

zJ − 1 z−1

N =

J −1 N N ( J −1)
h
z = Ch (N; J )z h : h=0

h=0

(5)

A formula (Ericksen, 1998) for the multinomial Ch (N; J ) coe3cients was developed by DeMoivre in 1756 as    h=J h − aJ + N − 1 N
 a Ch (N; J ) = ; (−1) N −1 a a=0 where h=J  is a Joor function for the largest integer not exceeding h=J . The coe3cients Ch (N; J ) of a Pascal–DeMoivre triangle with J = 3 are shown in Table 1, for row numbers 1 6 N 6 3. As seen in Table 1, uniform Ch (N; J ) distributions of unitary terms are created in row N = 1, since Ch (1; J ) = 1 for all 0 6 h 6 J − 1. With variable J rewritten for simplicity as n + 1, the unitary case of Eq. (5) reduces to   n+1 n n

−1 z = Ch (1; n + 1)z h = zh : (6) z−1 h=0 h=0 These unitary Ch (1; n + 1) sequences are used to generate the series in Eqs. (2) and (3). Additionally, the more general series of (4) can be obtained by using higher levels of N ¿ 1 for the Ch (N; n + 1) coe3cients. From Table 1 at (N; J ) = (2; 3), one such series (13 ∗ 1 + 33 ∗ 2 + 53 ∗ 3 + 73 ∗ 2 + 93 ∗ 1) can be solved for its sum of 1845, with an alternating signed series total of 365. All of these series arrangements in (2) – (4) may be interpreted as moment equations of the Ch (N; n + 1) distributions. Thus, the solutions to these sum of powers problems can be derived from the variety of moment equations presented in this paper.

4. Distribution moment types The Rth moment equations about a point (M -) for the Ch (N; n + 1) distributions are uniquely de0ned by (7). With either positive or alternating signed terms, these coe3cient moments cm from Ericksen (2000) are distinguished by their respective g

60

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

values of plus or minus one Nn
cm()(R; T ) = gh (h − (M − ))R Ch (N; n + 1); h=0

(7)

where  is an arbitrary separation from the mean M . Three distinct types of distribution moments are speci0ed in this paper by the notation T for types C; D; or E, according to conditions on g and n:   Type C for g = + 1;   Type T = Type D for g = − 1 and n = even; (8)    Type E for g = − 1 and n = odd: 5. Spaced coecient moments The coe3cient moments cm()(R; T ) in Eq. (7) can be calculated for a separation  of (M + k=d) with arbitrary values of k and d, to give a revised moment equation R    Nn
k k = gh h + Ch (N; n + 1); (9) cm()(R; T ) = cm M + d d (R;T ) h=0 where the mean M = Nn=2. After we multiply Eq. (9) by a factor dR , we are ready to de0ne the “spaced coe3cient moment” scm     Nn k Nn k R scm()(R; T ) = scm = d cm + + 2 d (R;T ) 2 d (R;T ) =

Nn
h=0

gh (dh + k)R Ch (N; n + 1):

(10)

Thus at each moment type (C; D; or E), the spaced coe3cient moment formulas (10) are the same as the general equations of (4). For the unitary case with N = 1, Eq. (10) yields the spaced sum of powers equations in (3). And under the special constraints of (N; d; k) = (1; 1; 0), Eq. (10) simpli0es further to the sum of powers equations in (2). 6. Exponential generating functions The moment equations may be derived from the expansions of their exponential generating functions. The Rth moments are the coe3cients of t R =R! in such an expansion. As shown by the equations in (11), each Rth spaced coe3cient moment scm is obtained as a coe3cient in the expansion of scm(egf; T ) . Using the relationship in Eq. (10), the exponential generating function scm(egf; T ) can be related to the Rth coe3cient moments cm(R; T ) with the variable t R replaced by (dt)R ∞ ∞
tR
(dt)R scm()(egf; T ) = scm()(R; T ) = cm()(R; T ) ; (11) R! R=0 R! R=0 where  = Nn=2 + (k=d).

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

61

Table 2 Variables FT and S(r; T ) , by type T Type T

FT

S(r; T )

Type C Type D Type E

(n + 1)N 1 (−dt=2)N (n + 1)N

((n + 1)r − 1)Br =r ((n + 1)r − 1)(2r − 1)Br =r ((n + 1)r − 2r + 1)Br =r

where Br is the rth Bernoulli number

The exponential generating functions for the coe3cient moments cm(R; T ) were derived in closed form in (Ericksen, 2000). The scm(egf; T ) equation in (12) extends these 0ndings by replacing the variable t R by (dt)R as described in (11). Here the exponent kt comes from ( − M ) dt in these closed form equations, with appropriate g values of ±1 for moment types (C; D; E) from (8)  N   (n+1) (n+1) dt Nn k e −1 g scm + = ekt : (12) 2 d (egf;T ) gedt − 1 To show that Eq. (12) is true, we apply Howard’s approach in (Howard, 1996) for the N = 1 case and extend it to the general N case by induction. So in the N = 1 case
n with sign term g = ± 1, the moment generating function for h=0 gh (dh + k)R becomes     (n+1) (n+1) dt n ∞
tR
e − 1 kt n k g + = e : scm gh e(dh+k)t = (13) 2 d (R;T ) R! h=0 gedt − 1 R=0 This formula is the same as the generating function for the Pascal–DeMoivre coef0cients (6) with the variable z being replaced by gedt and with the added location factor ekt . In fact, this proof was generalized in (Ericksen, 2000) for any polynomial P(gz) ≡
Nn h h h=0 g Ch z with sign term g and variable z, so that the moment generating function for the coe3cients of z h is simply   Nn ∞

tR Nn k + = (dh − k)R gh Ch = ekt P(gedt ): (14) scm 2 d (egf;T ) R=0 h=0 R! Ericksen (1998, 2000) also derived equivalent representations for exponential generating functions of the coe3cient moments cm(egf; T ) in an exponential sum format of Eq. (15). Again the requisite adjustments for (dt)R and Exp(kt) made the conversion to the format of the spaced coe3cient moments. Here the notation Exp(y) means ey .  

 ∞
(dt)r Nn k r ; (15) = FT Exp kt + N (−1) S(r; T ) scm + 2 d (egf;T ) r! r=1 where variables FT and S(r; T ) are given in Table 2 for each moment type T . The normalization factors FT represent the transformation of the probability moments m(R; T ) to the coe3cient moments scm(R; T ) , as outlined in (Ericksen, 1998, 2000). For moment equations with (N; d; k) = (N; 1; 0), the Rth moment about the origin for the Pascal–DeMoivre coe3cients is expressed by an expansion of (15) as a polynomial

62

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

Table 3 A(R; i) coe3cients of N i R 0 1 2

DR 1 2 12

C(N; 1; 0)

D(N; 1; 0)

E(N; 1; 0)

i=0

1

2

0

1

2

0

1

2

1 0 0

0 n n2 + 2n

0 0 3n2

1 0 0

0 n 3n2 + 6n

0 0 3n2

1 0 0

0 n n2 + 2n − 2

0 0 3n2

in N with coe3cients A(R; i) shown in Table 3:   Nn R

Nn scm = gh hR Ch (N; n + 1) = FT A(R; i) N i : 2 (R;T ) h=0 i=0

(16)

Tables of coe3cients, like Table 3, will display the numerators of the coe3cients in matrix or triangular form. The common denominator of the coe3cients in the Rth row will be listed separately as DR . A notation of T (N; d; k) will identify the moment type T as (C; D; E), along with the values of the variables (N; d; k) for each coe3cient table. With the aid of Eq. (15) and Table 2, we can readily decipher some of the ith term values for the A(R; i) coe3cients in Table 3. At i = R for example, the coe3cient A(R; R) of N R is always (n=2)R , which is 1 at R = 0 and n=2 at R = 1. Also the 0rst non-zero coe3cients of N or N 2 may be obtained from the list: (i) (i)th term of A(R; i) for T (N; 1; 0) (condition) 1

(−1)R S(R; T )

for R = even ¿ 0;

2

nR 2 S(R−1; T )

for R = odd ¿ 1:

(17)

The appearance of Bernoulli numbers in the exponential sum (15), Table 2, and list (17) provides motivation for the further study of Bernoulli polynomials. 7. Bernoulli and Euler polynomials Generalized Bernoulli polynomials Br(N ) (x) of NNorlund (1924) are created as the co(N ) e3cients in expansion (18) of the exponential generating function B(egf) (x). Simple

Bernoulli numbers Br of Zwillinger (1996) are de0ned by Br(1) (0) from Eq. (18) with x = 0. The 0rst few Bernoulli numbers Br for r ¿ 0 are {1; −(1=2); (1=6); 0; −(1=30); : : :} N  ∞
t tr (N ) (18) B(egf) (x) = t ext = (Br(N ) (x)) : e −1 r! r=0 According to Ericksen (2000), the generalized Bernoulli polynomials Br(N ) (x) can also be generated as the coe3cients in the expansion of an exponential sum composed only of simple Bernoulli numbers Br : 

r ∞ ∞

tr (N ) r Br t : (19) (Br (x)) = Exp xt − N (−1) r! r r! r=0 r=1

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

63

Additionally, an explicit form (20) of the generalized Bernoulli numbers Br(N ) at x = 0 was presented by Srivastava and Todorov (1988) as      q r N +r N +q−1
q
r! (N ) j (j)(r+q) ; (−1) (20) Br = r−q q j q=0 (r + q)! j=0 with B0(N ) = 1 since 0r = 1 at r = 0. Also, many explicit formulas for Bernoulli numbers Br in the N = 1 case are given in Dilcher’s paper (Dilcher, 1998). Next, a generalized Euler polynomial Er(N ) (x) will be de0ned as the coe3cients in (N ) the expansion of the exponential generating function E(egf) (x), as shown in Eq. (21). Following the simple Bernoulli number notation, we will identify “Euleros numbers E0r ” as the values of the simplest Euler polynomial Er(1) (x) at x = 0, whose initial values for r ¿ 0 are {1; −(1=2); 0; (1=4); 0; −(1=2); : : :} N  ∞
tr 2 (N ) (21) ext = (Er(N ) (x)) : E(egf) (x) = t r! e +1 r=0 Also, the generalized Euler polynomials Er(N ) (x) can be generated from the simple Bernoulli numbers Br or Euleros numbers E0r , according to Ericksen (2000): 

∞ ∞

tr Br t r r r (N ) (−1) (2 − 1) (Er (x)) = Exp xt − N r r! r! r=1 r=0

 r ∞ E0
t = Exp xt + N ; (22) (−1)r (r−1) 2 r! r=1 since E0(r−1) = 2(1 − 2r )Br =r, according to Abramowitz and Stegun (1972). In Dilcher’s paper (Dilcher, 1998), an explicit formula for the Euler polynomial at x = 0 is given for the N = 1 case as   q (−1)j r
q
(j)r : E0r = (23) q 2 j q=0 j=0 We state identities (24) from Abramowitz and Stegun (1972) and NNorlund (1924) that transform Bernoulli polynomials into power series with Bernoulli number coe3cients, because they will be useful in extending sum of powers equations throughout this paper     Q Q Q Q

p (B(Q−p) )t = (Bp )t (Q−p) : (24) BQ (t) = p p p=0 p=0 8. Moments by Bernoulli and Euler sums 8.1. Positive coe3cient moments—Type C Formula (12) for the exponential generating function of the spaced coe3cient moments is simpli0ed for positive coe3cient moments with g = 1 as   N  (n+1) dt Nn k −1 e scm + = ekt : (25) 2 d (egf;C) edt − 1

64

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

If we expand the numerator of (25) by the binomial theorem, we get    N   N N 1 Nn k kt
(−1)(N −i) e(n+1) idt : = dt e scm + 2 d (egf;C) e −1 i i=0

(26)

The Rth spaced coe3cient moments in (27) are obtained as the coe3cients in the expansion of the exponential generating function (26), after appropriately adjusting for the summation index when using Eq. (18). Thus formula (4) for
Nn R h=0 (dh + k) Ch (N; n + 1) can be evaluated as       N N
k R! Nn k (N ) R (N −i) + : B(R+N i(n + 1) + =d (−1) scm ) 2 d (R;C) (R + N )! i=0 d i (27) In the (N; d; k) = (1; d; k) case, we show in Eq. (28) the Hurwitz zeta equivalent of formula (27) as another representation of the sum of powers equation (3) for
n R h=0 (dh + k) , since B(R+1) (a) = − (R + 1)%(−R) (a)        n
k k n k R R + %(−R) ; = (dh + k) = d −%(−R) n + 1 + scm + 2 d (R;C) h=0 d d (28)
∞ where the Hurwitz zeta is %s (a) = %(s; a) = h=0 (h + a)(−s) . The validity and rationale for using the Hurwitz zeta version of the sum of powers formula (28) is made
∞ clear by the partition of the in0nite sums in identity (29), since h=0 dR (h + (n + 1)
∞ + (k=d))R = h=n+1 (dh + k)R with s = − R   n ∞ ∞


n k + = (dh + k)R = (dh + k)R − (dh + k)R : (29) scm 2 d (R;C) h=0 h=0 h=n+1 The occurrence of the generalized Bernoulli polynomial in Eq. (27) implies that a generalized Hurwitz zeta %s(N ) (a), in the next section, could also extend the sum of powers equations of (3) to the general formulas in (4). For the special case with (N; d; k) = (1; 1; 0), Eq. (27) reduces to the known formula
n (Zwillinger, 1996) of Eq. (30) for the sum of powers h=0 hR in (2) n n
1 = hR = (30) scm (B(R+1) (n + 1) − B(R+1) (0)): 2 (R;C) h=0 R+1 8.2. Alternating signed coe3cient moments—types (D; E) The exponential generating function for the spaced coe3cient moments with alternating signs from (12) with g = − 1 can be written as  N   (−1)n e(n+1) dt + 1 Nn k + = ekt : (31) scm 2 d (egf;T ) edt + 1

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

Expanding the numerator of Eq. (31) by the binomial theorem gives   N    N N Nn k 1 kt
+ (−1)in e(n+1)idt : scm = dt e 2 d (egf;T ) e +1 i i=0

65

(32)

The Rth spaced coe3cient moment in (33) is obtained as the coe3cient of t R =R! in the expansion of the exponential generating function (32), thus invoking Eq. (21). So
Nn Eq. (33) evaluates the general formula from (4) for h=0 (−1)h (dh + k)R Ch (N; n + 1) as       N N k dR
Nn k in : (33) ER(N ) i(n + 1) + = N (−1) scm + d 2 d (R;T ) 2 i=0 i In the (N; d; k) = (1; d; k) case, Eq. (34) gives a Dirichlet eta equivalent of formula
n (33) in another format of the sum of powers equation (3) for h=0 (−1)h (dh + k)R , since ER (a) = 2'(−R) (a)        k k n k + (−1)n '(−R) n + 1 + ; (34) + = dR '(−R) scm d d 2 d (R;T )
∞ where the Dirichlet eta is 's (a) = '(s; a) = h=0 (−1)h (h + a)(−s) . Using the approach in (29), Eq. (34) is veri0ed by the in0nite sums identity: n ∞ ∞


(−1)h (dh + k)R = (−1)h (dh + k)R + (−1)n (−1)h (dh + k)R : (35) h=0

h=0

h=n+1

Here too, the generalized Euler polynomials in Eq. (33) imply that a generalized ) Dirichlet eta '(N s (a) would also su3ce to extend the sum of powers formulas in (3) to the general case (4). Again in the special case of (N; d; k) = (1; 1; 0), Eq. (33) reduces to Eq. (36) for the alternating signed sum of powers from (2). This formula corresponds to a similar known formula (Zwillinger, 1996), diOering only by index selection n n
1 scm = (−1)h (h)R = (ER (0) + (−1)n ER (n + 1)): (36) 2 (R;T ) h=0 2 9. Lerch transcendentals The Lerch transcendental is an expression for the sums of reciprocal powers in (37), shown in the index format (s (z; x) of this report  ∞ %s (x) for z = + 1;
zl (s (z; x) = ((z; s; x) = (37) = s 's (x) for z = − 1: l=0 (l + x) The Hurwitz zeta %s (x) and the Dirichlet eta 's (x) are special cases of Lerch transcendentals, under the conditions for the z values shown in (37). These functions are generalized for power term N , as de0ned by  (N ) %s (x) for z = + 1; (s(N ) (z; x) = (38) ) '(N s (x) for z = − 1:

66

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79 Table 4 Variables S(r; T ) , by type T T

S(r; T ) by Zeta,Eta

S(r; T ) by Phi

C D E

((n + 1)r − 1)%(−r+1) (0) ((n + 1)r − 1)'(−r+1) (0) ((n + 1)r %(−r+1) (0)) + '(−r+1) (0)

((n + 1)r − 1)((−r+1) (1; 0) ((n + 1)r − 1)((−r+1) (−1; 0) ((n + 1)r ((−r+1) (1; 0)) + ((−r+1) (−1; 0)

In this paper, we will be concerned with the negative integer values for the variable s in these generalized Hurwitz zeta and the Dirichlet eta functions; i.e., with s taken as −R for R ¿ 0. Hurwitz Zeta. We showed in Eq. (28) that the Hurwitz zeta could be used in moment equations in the N = 1 case, when we can substitute the identity B(R+1) (x) = − (R + 1)%(−R) (x). As in the generalized Bernoulli forms (18) and (19) with N ¿ 1, (N ) a generalized Hurwitz zeta polynomial %(−R) (x) will be de0ned by the coe3cients in expansion (39) of the exponential generating functions in its closed form (39) and its exponential in0nite sum form (40) N  ∞
tr 1 (N ) (N ) %(egf) (x) = (%(−r+1) (x)) = t ext (39) r! e −1 r=0  

 ∞
N tr r t − N log t + N : (40) = Exp x− (−1) (%(−r+1) (0)) 2 r! r=2 Dirichlet Eta. We showed in Eq. (34) that the Dirichlet eta could be used in moment equations in the N = 1 case, when we can substitute the identity E(R) (x) = 2'(−R) (x). Following the generalized Euler polynomial format in (21) and (22) with N ¿ 1, a ) generalized Dirichlet eta polynomial '(N (−R) (x) will be de0ned by the coe3cients in the expansion (41) of the exponential generating functions in closed form (41) and exponential in0nite sum form (42)

t r  1 N ∞
(N ) (N ) '(egf) (x) = '(−r) (x) ext (41) = t r! e +1 r=0 

∞
tr r : (42) = Exp xt − N log (2) + N (−1) ('(−r+1) (0)) r! r=1 Moment generating functions. As was done in Eq. (15), we derive the exponential generating functions of the spaced coe3cient moments scm(egf; T ) in an exponential sum format:  

 ∞
(dt)r Nn k r ; (43) = FT Exp kt − N (−1) S(r; T ) scm + 2 d (egf;T ) r! r=1 with normalization factor FT from Table 2. And the variables S(r; T ) are given in two equivalent formats in Table 4 for each moment type T . Moment equations. Eq. (44) displays the Hurwitz zeta and its Lerch transcendental form of the general formula in (27) for the positively termed moment equation in (4)

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

of

67


Nn

+ k)R Ch (N; n + 1).       N N
k Nn k (N ) R (N −i) %(−R+N ) i(n + 1) + scm =d (−1) + 2 d (R;C) d i i=0

h=0 (dh

R

=d

N
i=0



(−1)

(N −i)

N i



  k (N ) : ((−R+N 1; i(n + 1) + ) d

(44)

Eq. (45) displays the Dirichlet eta and its Lerch transcendental version of the general
Nn h R formula in (33) for the alternating signed moments (4) of h=0 (−1) (dh + k) Ch (N; n + 1)       N N k Nn k ) R
in + '(N i(n + 1) + (−1) =d scm (−R) d d (R;T ) 2 i i=0 R

=d

N
i=0



(−1)

in

N i



(N ) ((−R)



k −1; i(n + 1) + d

 :

(45)

10. Faulhaber formulas 10.1. Faulhaber moment equation—type C
n In 1631, Faulhaber developed a formula for h=1 hR in the form of a polynomial in n with Bernoulli number coe3cients. His Eq. (46) needed only a sign adjustment in one Bernoulli term B1 to be valid   n R+1

1 R+1 R (R+1−p) (B(R+1−p) )np ; h = (−1) (46) R + 1 p p=1 h=1 since (1 + B1 ) = − B1 = 12 and since B(2i+1) = 0 for i ¿ 1. Later, in 1713, Jacques Bernoulli presented Eq. (46) in his Arc conjectandi, according to Nielsen (1923).
n In our current terminology, the Faulhaber formula for h=0 hR becomes   n R+1 R + 1
1 (B(R+1−p) )np ; = 0R + (−1)(R+1−p) (47) scm 2 (R;C) R + 1 p=1 p where 0R = 1 for R = 0 and otherwise 0R = 0, as with the Kronecker delta.
n The sum of powers formula for h=0 hR can also be displayed by (48) as a fully expanded power series in n with coe3cients *(R; i) . With denominators DR at each R level, Table 5 summarizes these well-known sum of powers equations (Edwards, 1987; Zwillinger, 1996) by listing the *(R; i) coe3cients for each ni term n n R R+1


= hR = (n + 1) a(R; i) ni = *(R; i) ni ; (48) scm 2 (R;C) h=0 i=0 i=0 where *(R; i) = (n + 1) a(R; i) , since the normalization factor FT = (n + 1) at N = 1.

68

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79 Table 5 *(R; i) for C(1,1,0) R

DR

i=0

0 1 2 3 4

1 2 6 4 30

1 0 0 0 0

1

2

3

4

5

1 1 1 0 −1

0 1 3 1 0

0 0 2 2 10

0 0 0 1 15

0 0 0 0 6

Table 6 a(R; i) for C(1; 1; 0) R

DR

i=0

0 1 2 3 4

1 2 6 4 30

1 0 0 0 0

1

2

3

4

0 1 1 0 −1

0 0 2 1 1

0 0 0 1 9

0 0 0 0 6

For the type C moments with (N; d; k) = (1; 1; 0), the formula for the ith terms of the coe3cients *(R; i) from Table

5 is the Faulhaber equation (47) with p = i. So (R+1−i) R+1 *(R; i) = 1=(R+1)(−1) B(R+1−i) for 1 6 i 6 R+1 with R ¿ 1, and *(0; i) = 1 1 for 0 6 i 6 1 at R = 0. In another power series arrangement for this (N; d; k) = (1; 1; 0) case, the coe3cients a(R; i) of ni of Eq. (48) are shown in Table 6. The formula (49) for the coe3cients a(R; i) in (48) is obtained from Eq. (30), by applying the identity (24) as in (Comtet, 1970), extracting the common factor (n + 1), expanding the remaining binomial (n + 1)R , and collecting the coe3cients of ni    R+1 i+q
1 R−i B(R−i−q) : a(R; i) = (49) R + 1 q=0 R − i − q i
n R In the case of (N; d; k) = (1; d; k), the equation for h=0 (dh + k) can be written as (50) with the aid of identity (24) and the generalized equation (27), evaluated at N = 1. Following Faulhaber’s transformation to obtain a polynomial in n, this formula (50) is somewhat simpli0ed by Eq. (51)    p  p   
R+1 n k dR R+1 k k scm + (B(R+1−p) ) = n+1+ − 2 d (R;C) R+1 p=0 d d p (50)
dR R+1 = R + 1 p=0  p  k : − d



R+1 p





(R+1−p)

(B(R+1−p) ) (−1)



k n+ d

p

(51)

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

69

Table 7 *(R; i) for C(1; d; k) Coe3cients a(R; i)

R

i

DR

0

0 1

1

1 1

1

0 1 2

2

2k d + 2k d

2

0 1 2 3

6

6k 2 d2 + 6dk + 6k 2 3d(d + 2k) 2d2

3

0 1 2 3 4

4

4k 3 2k(d + 2k)(d + k) d(d2 + 6dk + 6k 2 ) 2d2 (d + 2k) d3

4

0 1 2 3 4 5

30

30k 4 −d4 + 30d2 k 2 + 60dk 3 + 30k 4 30dk(d + 2k)(d + k) 10d2 (d2 + 6dk + 6k 2 ) 15d3 (d + 2k) 6d4

We note that Eq. (51) reduces to Eq. (47) with (N; d; k) = (1; 1; 0), since the term at p = 0 becomes 1=(R + 1)(B(R+1) )((−1)(R+1) − 1) = 0R . If we now expand the (n + (k=d))p term in (51) by the binomial theorem, we obtain the fully expanded polynomial in n for the (N; d; k) = (1; d; k) case: 

n k scm + 2 d

 (R;C)

=

n
h=0

(dh + k)R = (n + 1)

R
i=0

a(R; i) ni =

R+1
p=0

*(R; p) np ;

(52)

where
dR R−p *(R; p) = R + 1 s=0



R+1 s+p+1



s+p+1 p+1

 (B(R−p−s) )

 s k (−1)(R−p−s) d

for p ¿ 1; and *(R; 0) = k R . In these cases when (N; d; k) = (1; d; k), the Faulhaber equation is the sum of powers equation (3). The ith term formula (52) for the coe3cients *(R; i) is shown in Table 7.

70

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79 Table 8 a(R; i) for C(1; d; k) Coe3cients a(R; i)

R

i

DR

0

0

1

1

1

0 1

2

2k d

2

0 1 2

6

6k 2 d(d + 6k) 2d2

3

0 1 2 3

4

4k 3 2dk(d + 3k) d2 (d + 4k) d3

4

0 1 2 3 4

30

30k 4 d(−d3 + 30dk 2 + 60k 3 ) d2 (d2 + 30dk + 60k 2 ) 3d3 (3d + 10k) 6d4

And a list of a(R; i) coe3cients in moment equation (52), with the normalization factor (n + 1) extracted, is given in Table 8. In the most general (N; d; k) cases for positive coe3cient moments, the formula (27) can be written in Faulhaber style as 

Nn k scm + 2 d



N
R! = dR (−1)(N −i) (R + N )! i=0 (R;C) R+N


×

p=0



R+N p



N

 (N ) (B(R+N −p) )



i 

k i(n + 1) + d

p :

(53)

Here the generalized Bernoulli numbers Br(N ) in (53) are obtained as the coe3cients from the expansion of their exponential generating function (18), evaluated at x = 0. Also from Srivastava and Todorov (1988), the explicit form (20) for the generalized Bernoulli numbers Br(N ) can be inserted into Eq. (53) to give a comprehensive explicit expression for the generalized sum of powers of (4). For these positive coe3cient moments with arbitrary (N; d; k) values, the generalized sum of powers equations (4) are shown here in an equation array (54) at small R

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

71

Table 9 a(R; i) for D(1; 1; 0) R

DR

0 1 2 3 4

1 2 2 4 2

levels

i=0 1 0 0 0 0


Nn

1

2

3

4

0 1 1 0 −1

0 0 1 3 0

0 0 0 2 2

0 0 0 0 1

(dh + k)R Ch (N; n + 1)

R

Sum of powers =

0

(n + 1)N ;   1 (2k + dNn); (n + 1)N 2   1 (n + 1)N (12k 2 + 2d2 Nn + 12dkNn + d2 Nn2 + 3(dNn)2 ): 12

1 2

h=0

(54)

10.2. Alternating signed moment equations—types (D; E)
n h R For both moment types T = (D; E), Eq. (55) for h=0 (−1) h is the alternating signed moment complement to the Faulhaber formula (47) in terms of E0r , the rth Euler polynomial at x = 0   n R R 1 1
(R+n−p) scm (E0(R−p) ) np ; = E0 R + (−1) (55) 2 (R;T ) 2 2 p=0 p since (1 − 12 E0(0) ) = 12 E0(0) = 12 , and since E0(2i) = 0 for i ¿ 1. By extracting the term at p = 0 from the summation in (55), another Faulhaber version for R ¿ 0 can be displayed as    n R R 0 for n = even; 1
(R+n−p) p scm (E0(R−p) )n + = (−1) (56) 2 (R;T ) 2 p=1 p E0R for n = odd; since 12 E0R (1 + (−1)(R+n) ) is 0R at moment type D and (1 − 0R )E0R at type E. For (N; d; k) = (1; 1; 0) in type D moments, the coe3cients a(R; i) of ni from Eq. (57) are shown in a matrix format in Table 9 n n R

scm = (−1)h hR = a(R; i) ni : (57) 2 (R;D) h=0 i=0 The ith terms of the coe3cients a(R; i) from Table 9 for the type D moments when (N; d; k) = (1; 1; 0) are the coe3cients of ni in Eqs. (55) and (56) with p = i. Thus, for 0 6 i 6 R, we have a(R; i) = 12 (−1)(R−i) ( Ri )E0(R−i) for R ¿ 1, and 0R at R = 0. For (N; d; k) = (1; 1; 0) in type E moments, the coe3cients a(R; i) of ni in the middle equation of (58) would come from the expansion of (36), using identity (24). When the

72

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79 Table 10 *(R; i) for E(1; 1; 0) R

DR

0 1 2 3 4

1 2 2 4 2

i=0 0 −1 0 1 0

1

2

3

4

0 −1 −1 0 1

0 0 −1 −3 0

0 0 0 −2 −2

0 0 0 0 −1

factor +R is taken inside the summation, the right side of (58) generates the coe3cients *(R; i) of ni , shown in Table 10 n n R R−1


scm = (−1)h hR = +R a(R; i) ni = *(R; i) ni ; (58) 2 (R;E) h=0 i=0 i=0 where the factor +R in Eq. (58) has the value of − 12 (n + 1). Since the normalization factor from Table 2 is FE = +RN t N , the moment equations with R ¡ N will have zero valued coe3cients, which occurs at R = 0 when N = 1. The ith terms of the coe3cients *(R; i) from Table 10 for the type E moments when (N; d; k) = (1; 1; 0) are the coe3cients of ni in Eqs. (55) and (56) with p = i. Diagonals of zero elements in Table 10 occur at (R + i) = 0(mod 2) for (R − i) ¿ 1.
n The case of (N; d; k) = (1; d; k) for h=0 (−1)h (dh + k)R is evaluated by Eq. (59). The Faulhaber presentation as a polynomial in n is shown in Eq. (60) for the same alternating signed moments    p  p     R R n k k dR
k scm (E0(R−p) ) (−1)n n + 1 + = + + 2 d (R;T ) 2 p=0 p d d (59)   p  p    R R dR
k k = (E0(R−p) ) (−1)(R+n−p) n + : + 2 p=0 p d d (60) p

Also, for the (N; d; k) = (1; d; k) case in (60), we can expand the (n + (k=d)) term by the binomial theorem to get the fully expanded polynomial in n:   n R

n k scm = (−1)h (dh + k)R = *(R; p) np ; (61) + 2 d (R;T ) h=0 p=0 where
dR R−p *(R; p) = 2 s=0



R s+p



s+p p



 s k (E0(R−p−s) ) (−1)(R+n−p−s) d

for p ¿ 1;

with the (p = 0) term *(R; 0) = k R for n = even. And for n = odd, *(0; 0) = 0 at R = 0,
R−1 and *(R; 0) = − s=0 ( Rs )(−d)(R−s) k s E0(R−s) at R ¿ 1. In this special case with N = 1, the explicit formula (23) for the Euler polynomial E0r at x = 0 can be inserted into the equation for *(R; p) in (61) to give an explicit expression for the alternating signed sum of powers in (3).

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

73

Table 11 a(R; i) for D(1; d; k) R

i

DR

Coe3cients a(R; i)

0

0

1

1

1

0 1

2

2k d

2

0 1 2

2

2k 2 d(d + 2k) d2

3

0 1 2 3

4

4k 3 6dk(d + k) 3d2 (d + 2k) 2d3

4

0 1 2 3 4

2

2k 4 d(−d3 + 6dk 2 + 4k 3 ) 6d2 k(d + k) 2d3 (d + 2k) d4

For the sum of powers equations in (3) with arbitrary (d; k) values, a list of coe3cients a(R; i) is given in Table 11 when (N; d; k) = (1; d; k). Since a(R; i) = *(R; i) for type D moments, the coe3cients in Table 11 are derived from Eq. (61). For the sum of powers equations in (3) with arbitrary (d; k) values, a list of coe3cients *(R; i) from (61) is given in Table 12 for the type E moments when N = 1. And 0nally for the variables (N; d; k), the most general expression of a summation in Faulhaber style is given as   Nn k + scm 2 d (R;T )     p  N R N
R dR
k ) (in) (E0(N(R−p) = N (−1) ) i(n + 1) + ; (62) 2 i=0 d i p=0 p where the generalized Euleros numbers E0(Nr ) in (62) are the coe3cients in the expansion of the exponential generating function (21) at x = 0. For alternating signed moments with (N; d; k) of type D, the generalized sum of powers Eqs. (4) are given in the array of Eqs. (63) at various R levels.
Nn R Sum of powers = h=0 (−1)h (dh + k)R Ch (N; n + 1) 0

1;

1

1 2 (2k

2

2 1 4 (4k

+ dNn); + 2d2 Nn + 4dkNn + d2 Nn2 + (dNn)2 ):

(63)

74

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79 Table 12 *(R; i) for E(1; d; k) R

i

DR

Coe3cients *(R; i)

0

0

1

0

1

0 1

2

−d −d

2

0 1 2

2

−2dk −d(d + 2k) −d2

3

0 1 2 3

4

d(d2 − 6k 2 ) −6dk(d + k) −3d2 (d + 2k) −2d3

4

0 1 2 3 4

2

2dk(d2 − 2k 2 ) d(d3 − 6dk 2 − 4k 3 ) −6d2 k(d + k) −2d3 (d + 2k) −d4

For alternating signed moments with (N; d; k) of type E, the generalized sum of powers Eqs. (4) are given in the array of Eqs. (64) at various R levels. However, for simplicity of display, the sum of powers equations are reduced by the common factor +RN , indicated here by an asterix ∗ R

Sum of powers∗ =

0

1;

1

1 2 (2k

2

2 1 12 (12k


Nn

h=0

(−1)h (dh + k)R Ch (N; n + 1)=+RN

+ dNn); − 2d2 N + 2d2 Nn + 12dkNn + d2 Nn2 + 3(dNn)2 );

(64)

where +RN = (−d=2)N (n + 1)N in Table 2 as FE =t N . Since type E moments are zero when R ¡ N , formulas (62) and (64) apply for moment equations when R ¿ N .

11. Inverse equations for integer powers 11.1. Inverses from positive coe3cient moments The inverse of the matrix in Table 5 is well documented in (Edwards, 1987; Riordan, 1979). Similar to Riordan (1979), the integer power nR is represented by the power * * series (65) with the -(R; i) coe3cients of scmi . These combinatorial -(R; i) coe3cients

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

75

Table 13 * -(R; i) for C(1; 1; 0) R

DR

0 1 2 3 4

1 1 1 1 1

i=0 1 0 0 0 0

1

2

3

4

0 1 −1 1 −1

0 0 2 −3 4

0 0 0 3 −6

0 0 0 0 4

1

2

3

4

0 2 −1 1 −1

0 0 3 −3 4

0 0 0 4 −6

0 0 0 0 5

Table 14 a -(R; i) for C(1; 1; 0) R

DR

0 1 2 3 4

1 1 1 1 1

i=0 1 0 0 0 0

in (65) are listed in Table 13 at small R levels   n n R R−1 R

R * (R−1−i) scm = (−1) : n = -(R; i) scm 2 (i;C) i=0 2 (i+1;C) i i=0

(65)

Taking the inverse of the triangular matrix in Table 6, the coe3cients -i of scmi from the inverse Eq. (66) for nR can be displayed as the matrix in Table 14. Rearranging a Eq. (66), the power series expansion (67) shows the combinatorial nature of the -(R; i) coe3cients n R 1
a -(R; ; (66) nR = i) scm n + 1 i=0 2 (i;C)   n n R−1 R
nR (n + 1) = scm scm + (−1)(R−1−i) : (67) 2 (R;C) i=0 2 (i+1;C) i The two factor product in (67) can be written in terms of Bernoulli numbers and polynomials, by substituting any of the formulas for scm(n=2)(R; C) found in this paper. By calling upon Eq. (30) for instance, we obtain nR (n + 1) = B(R+1) (n + 1) − B(R+1) (0)    R−2 R
1 (R−1−i) (B(i+2) (n + 1) − B(i+2) (0)): + (−1) i+2 i i=0

(68)

11.2. Inverses from alternating signed moment equations By taking the inverse of the triangular matrix in Table 9, the coe3cients -i of scmi in the inverse Eq. (69) are described combinatorially and displayed in the matrix form

76

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79 Table 15 -(R; i) for D(1; 1; 0) R

DR

0 1 2 3 4

1 1 1 1 1

i=0 1 0 0 0 0

1

2

3

4

0 2 −2 3 −4

0 0 2 −3 6

0 0 0 2 −4

0 0 0 0 2

Table 16 -(R; i) for E(1; 1; 0) R

DR

0 1 2 3 4

1 1 1 1 1

i=0 −1 1 −1 1 −1

1

2

3

4

0 −2 2 −3 4

0 0 −2 3 −6

0 0 0 −2 4

0 0 0 0 −2

of Table 15 R

n =

R
i=0

-(R; i) scm

n 2

(i;D)

= scm

n 2

(R;D)

+

R
i=1

 (−1)

(R−i)

R i

 scm

n 2

(i;D)

: (69)

Taking the inverse of the triangular matrix in Table 10 with the *(0; 0) term at −1, the coe3cients -(R; i) of scmi in the inverse Eq. (70) have the combinatorial values shown in (70) and enumerated in Table 16   n n n R R R

(R−i) R scm = − scm − (−1) : n = -(R; i) scm 2 (i;E) 2 (R;E) i=0 2 (i;E) i i=0 (70) 11.3. Generating function for -(R; i) coe3cients For all type T moments with variables (N; d; k) = (1; 1; 0), the -(R; i) coe3cient values were listed in Tables (13–16) for various Rth moment levels. These coe3cients -(R; i) in the inverse Eqs. (65), (66), (69) and (70) can be obtained as coe3cients in the generating function equation: n ∞
: (71) -(R; gf) = -(R; i) xi with xi ≡ scm 2 (i;T ) i=0 These generating functions -(R; gf) are displayed in Table 17 with the variable x representing the spaced coe3cient moments scmi , as described by the transform Eq. (71). Here the notation 0R has the same meaning as a Kronecker delta .(R; 0) , whose value is one at R = 0 and zero elsewhere. The generating functions in Table 17 are useful as another method to generate the sum of powers Eqs. (2). In the cases with (N; d; k) = (1; 1; 0), coe3cients a(R; i) and

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

77

Table 17 Generating functions for coe3cients -(R; i) for T (1; 1; 0) Table

Type

Kind

13 14 15 16

C C D E

* a a; * *


∞

i=0

-(R; i) xi

−x(x − 1)R + x(R+1) −x(x − 1)R + x(R+1) + xR (x − 1)R + xR − (−1)R −(x − 1)R − xR


∞

i=0

-(R; i)

1 2 − 0R 1 − (−1)R + 0R −1

*(R; i) in the corresponding Eqs. (48), (57) and (58) for the spaced coe3cient moments scm(R; T ) can be derived by taking the inverse of the generating functions -(R; eg) of Table 17. 12. Chebyshev polynomials of the second kind 12.1. The function U(n; T ) (x) The sum of powers Eqs. (4) and (10) are created from exponential generating functions scm()(egf; T ) of (12) and (15). In this section, we consider the moment equations of (3) at N = 1. For simplicity, we will de0ne a polynomial n k + : (72) 2 d The interpretation for U(n; T ) (x) in (72) follows the same hyperbolic forms obtained for the moment generating functions from Ericksen (2000), with 0 = dt=2:  sinh(n + 1)0   for Type C;   sinh 0    cosh(n + 1)0 U(n; T ) (x) = (73) for Type D;  cosh 0       −sinh(n + 1)0 for Type E: cosh 0 U(n; T ) (x) = e−dt scm()(egf; T )

with  =

12.2. Recursion equation for U(n; T ) (x) From extensive empirical evidence, all of the formulas in (73) can be described at each consecutive n level by the recursion expression: U(n+1; T ) (x) = (2x)U(n; T ) (x) − U(n−1; T ) (x);

(74)

where x = cosh 0 with 0 = dt=2. The form of recursion Eq. (74) is the same as that for the Chebyshev polynomial of the second kind Un (x) in (Zwillinger, 1996). But instead of the common trigonometric cosine, here the variable x is the hyperbolic cosine. To satisfy the recusion Eq. (74) for each diOerent moment type T , only the two initial values of U(−1; T ) (x) and U(0; T ) (x) need to be speci0ed as given in Table 18.

78

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79 Table 18 Initial values of U(n; T ) (x) Type T

U(−1; T ) (x)

U(0; T ) (x)

Type C Type D Type E

0 sech 0 0

1 1 −tanh 0

where x = cosh 0 and 0 = (dt=2)

With type C moments for example, the initial values of U(n; C) (x) are exactly the same as those for the Chebyshev polynomial of the second kind Un (x) according to Zwillinger (1996). 12.3. The generating function (U(n; T ) (x))(gf; T ) The polynomials U(n; T ) (x) can also be created as the coe3cients of z n in the expansion of their generating function. Applying their initial values from Table 18, the generating functions at each moment type (8) are given by  1   for Type C;   1 − 2xz + z2     ∞
1 n (U(n; T ) (x))z = (U(n; T ) (x))(gf; T ) = ( 1 − 1) for Type D;  1 − 2xz + z 2 xz n=0       1   ( x2 −1=x) for Type E: 2 1 − 2xz + z As its hyperbolic equivalent, the type C function U(n; C) (x) has the same generating function as the Chebyshev polynomial of the second kind Un (x). References Abramowitz, M., Stegun, C.A. (Eds.), 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. 9th printing. Dover, New York, pp. 804 –806. Comtet, L., 1970. Analyse Combinatoire. Presses Universitaires De France, Paris, pp. 61– 62, 164 –165. Dilcher, K., 1998. Multiplikationstheorem fNur die Bernoullischen Polynome und explizite Darstellungen der Bernoullischen Zahlen. Abh. Math. Sem. Univ. Hamburg 59, 143–156. Edwards, A.W.F., 1987. Pascal’s Arithmetical Triangle. Charles Gri3n, London, pp. 82–84. Ericksen, L., 1998. The Pascal–De Moivre Triangles. Fibonacci Quart. 36 (1), 20–33. Ericksen, L., 2000. The Pascal–De Moivre moments and their generating functions. In: Howard, F.T. (Ed.), Applications of Fibonacci Numbers, Vol. 8. Kluwer Academic Publishers, Dordrecht. Howard, F.T., 1996. Sum of powers of integers via generating functions. Fibonacci Quart. 34 (3), 244–256. S ementaire des Nombres de Bernoulli. Gauthier-Villars, Paris, pp. 70, 188–191, Nielsen, N., 1923. TraitSe ElS 296. N NNorlund, N.E., 1924. Vorlesungen Uber DiOerenzenrechnung. Verlag Von Julius Springer, Berlin, pp. 19 –25, 121–129, 138–151. Riordan, J., 1979. Combinatorial Identities. Wiley, New York, pp. 159 –160.

L. Ericksen / Journal of Statistical Planning and Inference 101 (2002) 57–79

79

Srivastava, H.M., Todorov, P.G., 1988. An explicit formula for the generalized bernoulli polynomials. J. Math. Anal. Appl. 130 (2), 509–513. Zwillinger, D., 1996. CRC Standard Mathematical Tables And Formulae, 30th Edition. CRC Press, Boca Raton, pp. 17–22, 40 – 45, 476 – 491, 581.