Fast computation of Bernoulli, Tangent and Secant numbers

F ast Computation of Berno ulli, T angent and Secant Numbers Richard P . Brent and David Harvey Abstract W e consider the computatio n of Bernoulli, T angent (zag), and Secant (zig or Euler) numbers. In particular, we give asymp totically fast algorithms for compu t- ing the fir st n such nu mbers in O ( n 2 ( log n ) 2 + o ( 1 ) ) bit-oper ations. W e also give very short in-place algorithms for computing the first n T angent or Secant numbers in O ( n 2 ) integer operations. These algorithms are e xtremely simple, and fast for mod- erate v alues of n . They are f aster and use less space than the algorithms of Atkinson (for T angent and Secant number s) and Akiyama an d T an igawa (for Bernoulli nu m- bers). 1 Intr oduction The Bernoulli numbers are rational numbers B n defined by the generating function ∑ n ≥ 0 B n z n n ! = z exp ( z ) − 1 . (1) Bernoulli numbers are of interest in numb er theo ry and are related to special val- ues of the Riemann zeta fu nction (see § 2). They also occu r as coefficients in the Euler-Maclaurin for mula, so are relevant to h igh-precisio n computa tion of special function s [7, § 4.5]. Richard P . Brent Mathematical Sciences Institute, Australian N ational Univ ersity , Canberra, A CT 0200, Australia. e-mail: Tangent@ rpbrent.c om David Harv ey School of Mathematics and Statistics, Univ ersity of Ne w South W ales, Sydne y , NSW 2052, Australia. e-mail: D.Harvey @unsw.edu .au 1 2 Richard P . Brent and David Har vey It is sometimes con venient to consider s caled Berno ulli numbers C n = B 2 n ( 2 n ) ! , (2) with generating functio n ∑ n ≥ 0 C n z 2 n = z / 2 tanh ( z / 2 ) . (3) The gen erating functions (1) and (3) only differ by the single term B 1 z , since the other odd terms vanish. The T angent numbers T n , and Secant numbers S n , are defined by ∑ n > 0 T n z 2 n − 1 ( 2 n − 1 ) ! = tan z , ∑ n ≥ 0 S n z 2 n ( 2 n ) ! = sec z . (4) In this paper, which is based on an a talk given by the first autho r at a worksho p held to mark Jonathan Borwein’ s sixtieth birthday , we consider some algorithm s for computin g Bernoulli, T ange nt and Secant numb ers. For backgroun d, comb inatorial interpretatio ns, and ref erences, see Abram owitz and Stegun [1, Ch. 23] (where the notation differs from ours, e.g. ( − 1 ) n E 2 n is used for our S n ), and Sloane’ s [29] sequences A00036 7, A000 182, A0003 64. Let M ( n ) be the numb er o f bit-operation s requ ired for n -bit integer multipli- cation. Th e Sch ¨ onhag e-Strassen alg orithm [27] gi ves M ( n ) = O ( n log n log log n ) , and F ¨ urer [1 7] has recently given an improved bound M ( n ) = O ( n ( log n ) 2 log ∗ n ) . For simplicity we merely assume that M ( n ) = O ( n ( log n ) 1 + o ( 1 ) ) , w here th e o ( 1 ) term depend s on the pr ecise algorithm used for multiplication. For example, if the Sch ¨ onh age-Strassen algorithm is used, th en the o ( 1 ) ter m can be r eplaced b y log log log n / log log n . In §§ 2 – 3 we mention some relevant and gener ally well-k nown facts con cerning Bernoulli, T ang ent and Secant numbers. Recently , Harve y [20] showed that the sing le n umber B n can b e c omputed in O ( n 2 ( log n ) 2 + o ( 1 ) ) bit-oper ations using a modu lar algorithm. In this paper we sho w that all th e Bernoulli numbers B 0 , . . . , B n can be computed with the same complexity bound (and similarly for Secant and T an gent numbers). In § 4 we gi ve a relati vely simple algorithm that a chiev es the slightly weaker bound O ( n 2 ( log n ) 3 + o ( 1 ) ) . In § 5 we describe the improvement to O ( n 2 ( log n ) 2 + o ( 1 ) ) . The idea is similar to that espoused by Steel [31], althoug h we reduce the p roblem to division rather than multiplication. It is an open question whether the single number B 2 n can be computed in o ( n 2 ) bit-oper ations. In § 6 we give very short in-place algor ithms fo r co mputing the first n Secant o r T ang ent number s using O ( n 2 ) integer opera tions. These algor ithms are extremely simple, an d fast for moder ate values of n (say n ≤ 1000), althou gh asympto tically not as fast as the algor ithms giv en in §§ 4 – 5. Berno ulli numbers can easily be de- duced from the correspon ding T an gent numbers using the relation (14) below . Fast Compu tation of Bernoulli, T angent and Secan t Numbers 3 2 Bern oulli Numbers From the generating function (1) it is easy to see that the B n are ratio nal numb ers, with B 2 n + 1 = 0 if n > 0. The first fe w nonzero B n are: B 0 = 1, B 1 = − 1 / 2, B 2 = 1 / 6, B 4 = − 1 / 3 0, B 6 = 1 / 42 , B 8 = − 1 / 3 0, B 10 = 5 / 66 , B 12 = − 691 / 273 0, B 14 = 7 / 6. The denomin ators of the Bernoulli num bers ar e gi ven b y the V on S taudt – Clausen Theor em [12, 30], which states that B ′ 2 n : = B 2 n + ∑ ( p − 1 ) | 2 n 1 p ∈ Z . Here the sum is over all primes p for which p − 1 di vides 2 n . Since th e “cor rection” B ′ 2 n − B 2 n is ea sy to com pute, it m ight b e convenient in a progr am to store the integers B ′ 2 n instead of the rational numb ers B 2 n or C n . Euler fou nd that the Riemann zeta-fu nction fo r even non -negativ e integer argu- ments can be expressed in terms of Bernoulli numbers – the relation is ( − 1 ) n − 1 B 2 n ( 2 n ) ! = 2 ζ ( 2 n ) ( 2 π ) 2 n . (5) Since ζ ( 2 n ) = 1 + O ( 4 − n ) as n → + ∞ , we see that | B 2 n | ∼ 2 ( 2 n ) ! ( 2 π ) 2 n . From Stirling’ s approximation to ( 2 n ) !, the number of b its in the in teger part o f B 2 n is 2 n lg n + O ( n ) ( we write lg for log 2 ). Thus, it takes Ω ( n 2 log n ) space to store B 1 , . . . , B n . W e can n ot expect any algorithm to comp ute B 1 , . . . , B n in fewer th an Ω ( n 2 log n ) bit-op erations. Another c onnection betwe en the Bernou lli number s and th e Riemann zeta- function is the identity B n + 1 n + 1 = − ζ ( − n ) (6) for n ∈ Z , n ≥ 1. This follows fro m (5) and the fu nctional equatio n for the zeta- function , o r directly from a contour integral re presentation of the zeta-functio n [ 33]. From the ge nerating function ( 1), multiplying bo th sides by exp ( z ) − 1 and equ at- ing coefficients of z , we ob tain the recurrence k ∑ j = 0  k + 1 j  B j = 0 for k > 0 . (7) This recurr ence h as traditionally been used to compute B 0 , . . . , B 2 n with O ( n 2 ) arith- metic opera tions, for example in [22]. Howe ver, this is un satisfactory if floating- point numbers are used, becau se the recurrence is numerically unstab le : the relati ve 4 Richard P . Brent and David Har vey error in the com puted B 2 n is of ord er 4 n ε if the floatin g-poin t arithmetic ha s p reci- sion ε , i.e. lg ( 1 / ε ) b its. Let C n be d efined by (2). Then, multiplying each side of (3) by sinh ( z / 2 ) / ( z / 2 ) and equating coefficients gi ves the recu rrence k ∑ j = 0 C j ( 2 k + 1 − 2 j ) ! 4 k − j = 1 ( 2 k ) ! 4 k . (8) Using this recurrence to e valuate C 0 , C 1 , . . . , C n , the relati ve error in the compu ted C n is only O ( n 2 ε ) , which is satisfactory from a numerical point of vie w . Equation ( 5) can be used in several ways to compu te Bernou lli numb ers. If we want just one Bernoulli number B 2 n then ζ ( 2 n ) on the right-han d-side of (5) can be ev aluated to suf ficient accuracy using th e Eu ler produ ct: this is the “zeta-f unction” algorithm f or co mputing Bernou lli numb ers mentione d (with several referen ces to earlier work) by Harve y [20]. On the other hand, if we w ant several Bernoulli num- bers, then we can use the generatin g functio n π z tanh ( π z ) = − 2 ∞ ∑ k = 0 ( − 1 ) k ζ ( 2 k ) z 2 k , (9) computin g the coefficients of z 2 k , k ≤ n , to sufficient accu racy , as mentio ned in [3, 9, 10]. This is similar to the fast algorithm that we describe in § 4. The similarity can be seen more clearly if we replace π z by z in (9), gi ving z tanh ( z ) = − 2 ∞ ∑ k = 0 ( − 1 ) k ζ ( 2 k ) π 2 k z 2 k , (10) since it is the r ational n umber ζ ( 2 n ) / π 2 n that we need in order to compute B 2 n from (5). In fact, it is easy to see that (10) is equi valent to (3). There is a vast literature on Bernoulli, T angent an d Secant numbers. For exam- ple, the bibliog raphy of Dilcher and Sla vutskii [15] con tains more than 2000 items. Thus, we do no t attempt to giv e a comple te l ist o f references to related work. How- ev er, we briefly m ention the pro blem of c omputing irr e gular primes [8, 11], which are o dd prim es p such that p divides the class numbe r of the p -th cyclotomic field. The algorithms that we present in §§ 4 – 5 belo w are not suitable for this task because they take too much memory . It is much more space-efficient to use a mo dular algo- rithm where the com putations ar e perform ed modulo a single pr ime (o r m aybe the produ ct of a small num ber of primes), as in [8, 11, 14, 20]. Space can also be saved by the technique of “multisectioning”, which is described by Cran dall [13, § 3.2] and Hare [19]. Fast Compu tation of Bernoulli, T angent and Secan t Numbers 5 3 T angent and Secant Numbers The T angent numbers T n ( n > 0 ) (also called Zag numb ers) are defined by ∑ n > 0 T n z 2 n − 1 ( 2 n − 1 ) ! = tan z = sin z cos z . Similarly , the S ecant nu mbers S n ( n ≥ 0 ) (also called Euler o r Zig number s) are defined by ∑ n ≥ 0 S n z 2 n ( 2 n ) ! = sec z = 1 cos z . Unlike th e Bernou lli numb ers, the T angent and Secant numbers are positive integers. Because tan z and sec z have p oles at z = π / 2, we expect T n to grow roughly like ( 2 n − 1 ) ! ( 2 / π ) n and S n like ( 2 n ) ! ( 2 / π ) n . T o obtain mo re precise estimates, let ζ 0 ( s ) = ( 1 − 2 − s ) ζ ( s ) = 1 + 3 − s + 5 − s + · · · be the odd zeta-fun ction. Then T n ( 2 n − 1 ) ! = 2 2 n + 1 ζ 0 ( 2 n ) π 2 n ∼ 2 2 n + 1 π 2 n (11) (this can be proved in the same way as Euler’ s relation (5) for the Berno ulli nu m- bers). W e also have [1, (23.2.22 )] S n ( 2 n ) ! = 2 2 n + 2 β ( 2 n + 1 ) π 2 n + 1 ∼ 2 2 n + 2 π 2 n + 1 , (12) where β ( s ) = ∞ ∑ j = 0 ( − 1 ) j ( 2 j + 1 ) − s . (13) From (5) and (11), we see that T n = ( − 1 ) n − 1 2 2 n ( 2 2 n − 1 ) B 2 n 2 n . (14) This can also be proved d irectly , with out inv olving the zeta-func tion, by using the identity tan z = 1 tan z − 2 tan ( 2 z ) . Since T n ∈ Z , it follows from (14) that the odd pr imes in the denomina tor of B 2 n must divide 2 2 n − 1 . This is compatible with the V on Staudt–Clausen theorem, since ( p − 1 ) | 2 n implies p | ( 2 2 n − 1 ) by Fermat’ s little theorem. T n has about 4 n more bits than ⌈ B 2 n ⌉ , but bo th h av e 2 n lg n + O ( n ) bits, so asymp- totically there is not much dif feren ce between the sizes of T n and ⌈ B 2 n ⌉ . Thus, if our 6 Richard P . Brent and David Har vey aim is to compute B 2 n , we d o not lo se much by fir st comp uting T n , an d this may be more con venient since T n ∈ Z , B 2 n ∈ Q . 4 A Fast Al gorithm f or Bernoulli Numbers Harvey [20] showed how B n could be computed exactly , using a mod ular algorith m and the Chinese remainder theorem, in O ( n 2 ( log n ) 2 + o ( 1 ) ) bit-oper ations. The same complexity can be obtained using (5) an d the Euler prod uct for th e zeta-f unction (see the discussion in Harvey [20, § 1]). In th is section we show how to comp ute all of B 0 , . . . , B n with almost the same complexity bound ( only larger b y a factor O ( lo g n ) ). In § 5 we giv e an ev en faster algorithm , which av oids the O ( log n ) factor . Let A ( z ) = a 0 + a 1 z + a 2 z 2 + · · · be a power series with coe ffi cients in R , with a 0 6 = 0. Let B ( z ) = b 0 + b 1 z + · · · be the r ecip r oc al power series, so A ( z ) B ( z ) = 1. Using the FFT , we can multiply p olynomia ls of degree n − 1 with O ( n log n ) real operation s. Usi ng Newton’ s metho d [24, 28], we can compute b 0 , . . . , b n − 1 with the same complexity O ( n log n ) , up to a constant factor . T akin g A ( z ) = ( e xp ( z ) − 1 ) / z and working with N -b it flo ating-po int numbe rs, where N = n lg ( n ) + O ( n ) , we ge t B 0 , . . . , B n to sufficient accuracy to d educe the exact (ration al) resu lt. (Alternatively , use (3) to av oid comp uting the terms with odd subscripts, since the se vanish excep t for B 1 .) The work in volved is O ( n log n ) floating-poin t oper ations, each of which can be don e with N -b it ac- curacy in O ( n ( log n ) 2 + o ( 1 ) ) bit-oper ations. Thus, overall we get B 0 , . . . , B n with O ( n 2 ( log n ) 3 + o ( 1 ) ) bit-op erations. Similar ly for Secan t an d T a ngent nu mbers. W e omit a p recise specification of N and a detailed error analysis of the algor ithm, since it is improved i n the following section. 5 A Faster Algorithm for T angent and Bernoulli Numbers T o improve the algorithm of § 4 for B erno ulli nu mbers, we use the “Kronecker– Sch ¨ onhage trick ” [7, § 1.9]. In stead of working with power serie s A ( z ) (or po ly- nomials, which can be regarded as trunca ted power series), we work with b inary number s A ( z ) where z is a suitable (negative) po wer of 2. The id ea is to comp ute a single re al n umber A which is defined in such a w ay that th e num bers that we want to co mpute ar e enco ded in the bin ary rep resentation of A . For example, consider the series ∑ k > 0 k 2 z k = z ( 1 + z ) ( 1 − z ) 3 , | z | < 1 . Fast Compu tation of Bernoulli, T angent and Secan t Numbers 7 The r ight-hand side is an easily- computed ratio nal function of z , say A ( z ) . W e use decimal rather than binary for expository purposes. W ith z = 10 − 3 we easily find A ( 10 − 3 ) = 10010 00 99700 2999 = 0 . 001 004 009 01 6 025 036 049 06 4 081 100 · · · Thus, if we a re interested in the fin ite sequen ce of square s ( 1 2 , 2 2 , 3 2 , . . . , 10 2 ) , it is sufficient to compute A = A ( 10 − 3 ) co rrectly rounded to 30 decimal places, and we can then “read off ” the squares from the decimal representatio n of A . Of course, this example is purely fo r illustrative p urposes, because it is easy to compute the sequence of square s d irectly . Howe ver , we use the same idea to com- pute T angen t numb ers. Suppose we want the fir st n T angent numbers ( T 1 , T 2 , . . . , T n ) . The generatin g functio n tan z = ∑ k ≥ 1 T k z 2 k − 1 ( 2 k − 1 ) ! giv es us alm ost what we need , but no t q uite, b ecause the coefficients are rationa ls, not integers. Instead, consider ( 2 n − 1 ) ! tan z = n ∑ k = 1 T ′ k , n z 2 k − 1 + R n ( z ) , (15) where T ′ k , n = ( 2 n − 1 ) ! ( 2 k − 1 ) ! T k (16) is an integer for 1 ≤ k ≤ n , and R n ( z ) = ∞ ∑ k = n + 1 T ′ k , n z 2 k − 1 = ( 2 n − 1 ) ! ∞ ∑ k = n + 1 T k z 2 k − 1 ( 2 k − 1 ) ! (17) is a rem ainder term which is small if z is sufficiently small. Th us, choosing z = 2 − p with p sufficiently large, the first 2 n p binar y places o f ( 2 n − 1 ) ! tan z de- fine T ′ 1 , n , T ′ 2 , n , . . . , T ′ n , n . Once we hav e compu ted T ′ 1 , n , T ′ 2 , n , . . . , T ′ n , n it is easy to deduce T 1 , T 2 , . . . , T n from T k = T ′ k , n ( 2 n − 1 ) ! / ( 2 k − 1 ) ! . For t his idea to work , tw o conditio ns must be satisfied. First, we need 0 ≤ T ′ k , n < 1 / z 2 = 2 2 p , 1 ≤ k ≤ n , (18) so w e can read of f the T ′ k , n from the binary representation o f ( 2 n − 1 ) ! tan z . Since we have a g ood asympto tic estimate for T k , it is no t hard to choose p sufficiently large for this condition to hold. Second, we need the rem ainder term R n ( z ) to be sufficiently small that it does not influence the estimation of T ′ n , n . A sufficient condition is 8 Richard P . Brent and David Har vey 0 ≤ R n ( z ) < z 2 n − 1 . (19) Choosing z suf ficiently small (i.e. p sufficiently large) guarantee s that condition (19) holds, since R n ( z ) is O ( z 2 n + 1 ) as z → 0 with n fixed. Lemmas 3 and 4 below give suffi cient cond itions for (18) and (19) to hold. Lemma 1. T k ( 2 k − 1 ) ! ≤  2 π  2 ( k − 1 ) for k ≥ 1 . Pr o of. From (11), T k ( 2 k − 1 ) ! = 2  2 π  2 k ζ 0 ( 2 k ) ≤ 2  2 π  2 k ζ 0 ( 2 ) ≤  2 π  2 k π 2 4 =  2 π  2 ( k − 1 ) ⊓ ⊔ Lemma 2. ( 2 n − 1 ) ! ≤ n 2 n − 1 for n ≥ 1 . Pr o of. ( 2 n − 1 ) ! = n n − 1 ∏ j = 1 ( n − j )( n + j ) = n n − 1 ∏ j = 1 ( n 2 − j 2 ) ≤ n 2 n − 1 with equality iff n = 1. ⊓ ⊔ Lemma 3. I f k ≥ 1 , n ≥ 2 , p = ⌈ n lg ( n ) ⌉ , z = 2 − p , and T ′ k , n is as in (16), then z ≤ n − n and T ′ k , n < 1 / z 2 . Pr o of. W e have z = 2 − p = 2 −⌈ n lg ( n ) ⌉ ≤ 2 − n lg ( n ) = n − n , which p roves the first p art of the Lemma. Assume k ≥ 1 and n ≥ 2. From Lemma 1, we have T ′ k , n ≤ ( 2 n − 1 ) !  2 π  2 ( k − 1 ) ≤ ( 2 n − 1 ) ! , and from Lemma 2 it follows that T ′ k , n ≤ n 2 n − 1 < n 2 n . From the first part of the Lemma, n 2 n ≤ 1 / z 2 , so the second part follows. ⊓ ⊔ Lemma 4. I f n ≥ 2 , p = ⌈ n lg ( n ) ⌉ , z = 2 − p , an d R n ( z ) is as define d in (17), then 0 < R n ( z ) < 0 . 1 z 2 n − 1 . Pr o of. Since all the terms in the sum defining R n ( z ) are positi ve, it is im mediate that R n ( z ) > 0. Fast Compu tation of Bernoulli, T angent and Secan t Numbers 9 Since n ≥ 2, we have p ≥ 2 and z ≤ 1 / 4. No w , using Lemma 1, R n ( z ) = ∞ ∑ k = n + 1 T ′ k , n z 2 k − 1 ≤ ( 2 n − 1 ) ! ∞ ∑ k = n + 1  2 π  2 ( k − 1 ) z 2 k − 1 ≤ ( 2 n − 1 ) !  2 π  2 n z 2 n + 1 1 +  2 z π  2 +  2 z π  4 + · · · ! ≤ ( 2 n − 1 ) !  2 π  2 n z 2 n + 1 , 1 −  2 z π  2 ! . Since z ≤ 1 / 4, we ha ve 1 / ( 1 − ( 2 z / π ) 2 ) < 1 . 0 26. Also, from Lemma 2, ( 2 n − 1 ) ! ≤ n 2 n − 1 . Thus, we have R n ( z ) z 2 n − 1 < 1 . 026 n 2 n − 1  2 π  2 n z 2 . Now z 2 ≤ n − 2 n from the first part of Lemma 3, so R n ( z ) z 2 n − 1 < 1 . 026 n  2 π  2 n . (20) The right-hand side is a mo notonic de creasing f unction o f n , so is bounded above by its value when n = 2, gi ving R n ( z ) / z 2 n − 1 < 0 . 1 . ⊓ ⊔ A h igh-level description of the resulting Algorithm FastT ang entNumb ers is gi ven in Figure 1. The algor ithm computes the T ange nt numbers T 1 , T 2 , . . . , T n using the Kronecker-Sch ¨ onhage trick as described above, and deduces the Bernoulli numbers B 2 , B 4 , . . . , B 2 n from the relation (14). Input: integer n ≥ 2 Output: T angent numbers T 1 , . . . , T n and (optional) Bernou lli numbers B 2 , B 4 , . . . , B 2 n p ← ⌈ n lg ( n ) ⌉ z ← 2 − p S ← ∑ 0 ≤ k < n ( − 1 ) k z 2 k + 1 × ( 2 n ) ! / ( 2 k + 1 ) ! C ← ∑ 0 ≤ k < n ( − 1 ) k z 2 k × ( 2 n ) ! / ( 2 k ) ! V ← ⌊ z 1 − 2 n × ( 2 n − 1 ) ! × S / C ⌉ (here ⌊ x ⌉ means roun d x t o neares t integer ) Extract T ′ k , n = T k ( 2 n − 1 ) ! / ( 2 k − 1 ) !, 1 ≤ k ≤ n , from the binary representation of V T k ← T ′ k , n × ( 2 k − 1 ) ! / ( 2 n − 1 ) !, k = n , n − 1 , . . . , 1 B 2 k ← ( − 1 ) k − 1 ( k × T k / 2 2 k − 1 ) / ( 2 2 k − 1 ) , k = 1 , 2 , . . . , n (optional) ret urn T 1 , T 2 , . . . , T n and (optional) B 2 , B 4 , . . . , B 2 n Fig. 1 Algorithm FastT angentNumber s (also optionally computes Bernoulli numbers) 10 Richard P . Brent and David Har vey In o rder to achieve the best comp lexity , the alg orithm m ust be implemented carefully using binary arithmetic. The computa tions of S (an approximatio n to ( 2 n ) ! sin z ) a nd C (an ap proximatio n to ( 2 n ) ! cos z ) inv olve compu ting ratios of fac- torials such as ( 2 n ) ! / ( 2 k ) !, where 0 ≤ k ≤ n . This can b e done in time O ( n 2 ( log n ) 2 ) by a straightfor ward algorithm. The N -bit division to compu te S / C (an app roxima- tion to tan z ) can be done in time O ( N log ( N ) log log ( N )) by the Sch ¨ onh age–Strassen algorithm comb ined with Ne wton’ s m ethod [7, § 4.2 .2]. Her e it is sufficient to ta ke N = 2 n p + 2 = 2 n 2 lg ( n ) + O ( n ) . Note that V = n ∑ k = 1 2 2 ( n − k ) p T ′ k , n (21) is just the finite sum in (15) scaled b y z 1 − 2 n (a power of two), and the integers T ′ k , n can simply be “re ad off ” from the b inary re presentation of V in n blo cks of 2 p con secutiv e bits. The T ′ k , n can th en be scaled by ratios of factorials in time O ( n 2 ( log n ) 2 + o ( 1 ) ) to giv e the T angent nu mbers T 1 , T 2 , . . . , T n . The correctn ess of the computed T ang ent numbers follows f rom Lemmas 3 – 4, apart from possible er rors introduce d by S / C being only an approximatio n to tan ( z ) . Lemma 5 shows that this erro r is suf ficiently small. Lemma 5. S uppose th at n ≥ 2 , z, S and C as in Algorithm F astT angentNum bers. Then z 1 − 2 n ( 2 n − 1 ) !     S C − tan z     < 0 . 0 2 . (22) Pr o of. W e use the inequality     A B − A ′ B ′     ≤ | A | · | B − B ′ | + | B | · | A − A ′ | | B | · | B ′ | . (23) T ake A = sin z , B = co s z , A ′ = S / ( 2 n ) !, B ′ = C / ( 2 n ) ! in (2 3). Since n ≥ 2 we h av e 0 < z ≤ 1 / 4. Th en | A | = | sin z | < z . Also, | B | = | cos z | > 31 / 3 2 from the T aylo r series c os z = 1 − z 2 / 2 + · · · , which has terms of alternating sign and decr easing magnitud e. By similar arguments, | B ′ | ≥ 31 / 32, | B − B ′ | < z 2 n / ( 2 n ) !, and | A − A ′ | < z 2 n + 1 / ( 2 n + 1 ) !. Combinin g these ineq ualities and using (23), we obtain     S C − tan z     < 6 · 32 · 32 5 · 31 · 31 z 2 n + 1 ( 2 n ) ! < 1 . 28 z 2 n + 1 ( 2 n ) ! . Multiplying bo th sides b y z 1 − 2 n ( 2 n − 1 ) ! and using 1 . 2 8 z 2 / ( 2 n ) ≤ 0 . 02, we obtain the inequality (22). This completes the proof of Lemma 5. ⊓ ⊔ In vie w of the constant 0 . 0 2 in ( 22) and the co nstant 0 . 1 in Lemma 4, the effect of all sources of error in computing z 1 − 2 n ( 2 n − 1 ) ! tan z is at most 0 . 12 < 1 / 2 , which is Fast Compu tation of Bernoulli, T angent and Secan t Numbers 11 too s mall to change the computed integer V , that is to say , t he computed V is indeed giv en by (21). The c omputation of the Bernoulli n umbers B 2 , B 4 , . . . , B 2 n from T 1 , . . . , T n , is straightfor ward (details depen ding on exactly how rational numbers are to be repre- sented). The entire computa tion takes time O ( N ( lo g N ) 1 + o ( 1 ) ) = O ( n 2 ( log n ) 2 + o ( 1 ) ) . Thus, we have proved: Theorem 1. The T angent numbers T 1 , . . . , T n and Bernoulli numbers B 2 , B 4 , . . . , B 2 n can be computed in O ( n 2 ( log n ) 2 + o ( 1 ) ) bit-operations using O ( n 2 log n ) space. A small modificatio n of the above ca n be used to comp ute the Secant number s S 0 , S 1 , . . . , S n in O ( n 2 ( log n ) 2 + o ( 1 ) ) bit-o perations and O ( n 2 log n ) spac e. The boun d on T an gent numbers gi ven by Lemma 1 can be replaced by the bound S n ( 2 n ) ! ≤ 2  2 π  2 n + 1 which follows from (12) since β ( 2 n + 1 ) < 1. W e remark th at an efficient imple mentation of Algorithm F astT an gentNumb ers in a high-level lang uage such as Sage [32] or Magma [5] is nontr i vial, because it r equires access to th e intern al binary r epresentation of high-precision in tegers. Everything can b e done u sing (implicitly scaled ) integer arithm etic – ther e is n o need for floating-point – b ut for the sake of clarity we did not include the scaling in Figure 1. If floating-po int arithm etic is u sed, a precision of N bits is sufficient, where N = 2 n p + 2 . Comparing our Algor ithm FastT angentNum bers with Harvey’ s modular algo- rithm [2 0], we see that ther e is a spac e-time trade-off: Harve y’ s algorithm uses less space (by a factor of order n ) to compute a single B n , b ut more time (again by a fac- tor of o rder n ) to comp ute all of B 1 , . . . , B n . Harve y’ s algorithm has b etter locality and is readily parallelisable. In the following section we gi ve much simpler algorithms which are fast enoug h for most practical purposes, and are based on three-term recurren ce relations. 6 Algorithms Based On Thr ee-T erm Recurren ces Akiyama and T anig awa [21] g av e an algor ithm for compu ting Bernoulli numbers based on a three-term recurren ce. Ho wever , it is only useful for e xact computations, since it is numerically unstable if applied using floating-poin t arithmetic. It is f aster to use a stable recur rence for compu ting T angent number s, a nd then d educe the Bernoulli numbers from (14). 12 Richard P . Brent and David Har vey 6.1 Bernoulli and T angent numbers W e now g i ve a stab le three-term r ecurrenc e and correspo nding in-place algorithm for c omputing T angent nu mbers. The algorithm is perfectly stable since all o pera- tions are on positi ve integers and there is no can cellation. Also, it in volv es less arith- metic than the Akiyam a-T an igawa algorith m. This is pa rtly b ecause th e operation s are on integers rather th an ratio nals, and p artly b ecause there are fewer operation s since we take advantage of zeros. Bernoulli nu mbers can be computed using Alg orithm T angen tNumbers and the relation (14). The time required for the application of (14) is negligible. The r ecurrenc e (24) th at we u se was given by Buck holtz an d Kn uth [ 23], but they did no t g iv e our in-p lace Alg orithm T angen tNumbers explicitly . Related recurren ces with applications to parallel computatio n were considered by Hare [ 19]. Write t = tan x , D = d / d x , so Dt = 1 + t 2 and D ( t n ) = nt n − 1 ( 1 + t 2 ) for n ≥ 1. It is clear that D n t is a polynomial in t , say P n ( t ) . Write P n ( t ) = ∑ j ≥ 0 p n , j t j . Then deg ( P n ) = n + 1 and, from the formula for D ( t n ) , p n , j = ( j − 1 ) p n − 1 , j − 1 + ( j + 1 ) p n − 1 , j + 1 . (24) W e are in terested in T k = ( d / d x ) 2 k − 1 tan x | x = 0 = P 2 k − 1 ( 0 ) = p 2 k − 1 , 0 , which c an be computed from the recurr ence (24) in O ( k 2 ) o perations using the obvious b oundar y condition s. W e save work by noticing that p n , j = 0 if n + j is even. The resulting algorithm is given in Figure 2. Input: positiv e integer n Output: T angent numbers T 1 , . . . , T n T 1 ← 1 f or k fr om 2 to n T k ← ( k − 1 ) T k − 1 f or k fr om 2 to n f or j fr om k to n T j ← ( j − k ) T j − 1 + ( j − k + 2 ) T j ret urn T 1 , T 2 , . . . , T n . Fig. 2 Algorithm T angentNumbers The first for loop initializes T k = p k − 1 , k = ( k − 1 ) !. The variable T k is th en u sed to store p k , k − 1 , p k + 1 , k − 2 , . . . , p 2 k − 2 , 1 , p 2 k − 1 , 0 at successive iterations of the second for loop. Thus, when the algorith m terminates, T k = p 2 k − 1 , 0 , as expected. The p rocess in the case n = 3 is illustrated in Figure 3, where T ( m ) k denotes the value of the variable T k at successive iterations m = 1 , 2 , . . . , n . It is instructive to compare a similar Figure for the Akiyama- T anigawa algorithm in [21]. Algorithm T angentNum bers takes Θ ( n 2 ) o perations on positi ve integers. The in- tegers T n have O ( n log n ) bits, othe r integers have O ( log n ) bits. Thu s, the overall Fast Compu tation of Bernoulli, T angent and Secan t Numbers 13 T ( 1 ) 1 = p 0 , 1 ւ ց T ( 1 ) 1 = p 1 , 0 T ( 1 ) 2 = p 1 , 2 ց ւ ց T ( 2 ) 2 = p 2 , 1 T ( 1 ) 3 = p 2 , 3 ւ ց ւ T ( 2 ) 2 = p 3 , 0 T ( 2 ) 3 = p 3 , 2 ց ւ T ( 3 ) 3 = p 4 , 1 ւ T ( 3 ) 3 = p 5 , 0 Fig. 3 Dataflow in Algorithm T angentNumbers for n = 3 complexity is O ( n 3 ( log n ) 1 + o ( 1 ) ) bit-o perations, or O ( n 3 log n ) word-op erations if n fits in a single word. The algo rithm is not optimal, but it is good in practice f or modera te v alues of n , and much simpler than a symptotically faster algorithms such as those de scribed in §§ 4 – 5. For examp le, usin g a straightforward Magma implemen tation of Alg orithm T ang entNumbe rs, we co mputed the first 1000 T angen t n umbers in 1.50 sec o n a 2.26 GHz I ntel Cor e 2 Duo . For com parison, it takes 1.92 sec for a sing le N -bit division computing T in Algorith m F astT an gentNumb ers (wh ere N = 1993156 8 correspo nds to n = 1 000). Thu s, w e exp ect the cr ossover point where Algor ithm FastT angentNu mbers actually becomes faster to be sligh tly larger than n = 1000 (but dependent on implementation details). 6.2 Secant numbers A similar algo rithm may be used to compute Secant numb ers. Let s = sec x , t = tan x , and D = d / d x . Then Ds = s t , D 2 s = s ( 1 + 2 t 2 ) , and in general D n s = sQ n ( t ) , where Q n ( t ) is a polyn omial of degree n in t . The Secant n umbers are giv en by S k = Q 2 k ( 0 ) . Let Q n ( t ) = ∑ k ≥ 0 q n , k t k . From D ( s t k ) = st k + 1 + k s t k − 1 ( 1 + t 2 ) we obtain the three-term recurr ence q n + 1 , k = kq n , k − 1 + ( k + 1 ) q n , k + 1 for 1 ≤ k ≤ n . (25) By av oiding the comp utation o f terms q n , k that are known to b e zero ( n + k odd), and orderin g the computation in a manner analogous to that used f or Algorithm T ang ent- Numbers, we obtain Algorithm SecantNumbers (see Figure 4), which computes the Secant number s in place using non -negativ e integer arithmetic. 14 Richard P . Brent and David Har vey Input: positiv e integer n Output: Secant numbers S 0 , S 1 . . . , S n S 0 ← 1 f or k fr om 1 to n S k ← kS k − 1 f or k fr om 1 to n f or j fr om k + 1 to n S j ← ( j − k ) S j − 1 + ( j − k + 1 ) S j ret urn S 0 , S 1 , . . . , S n . Fig. 4 Algorithm SecantNumbers 6.3 Comparison with Atkinson’ s algorit hm Atkinson [ 2] g av e an elegant algor ithm for com puting both the T angent n umbers T 1 , T 2 , . . . , T n and the Secant nu mbers S 0 , S 1 , . . . , S n using a “Pascal’ s triang le” style of alg orithm tha t only in volves a dditions of non- negati ve integers. Since a trian- gle with 2 n + 1 rows in inv olved, Atkinson’ s algorithm requ ires 2 n 2 + O ( n ) integer additions. Th is can b e com pared with n 2 / 2 + O ( n ) add itions and n 2 + O ( n ) multi- plications (by small integers) for our Algorithm T an gentNumb ers, and similarly for Algorithm SecantNumb ers. Thus, we might expect Atkinson’ s algor ithm to be slower than Algorith m T an- gentNumb ers. Computation al experiments confirm this. W ith n = 1000, Algorithm T ang entNumbe rs prog rammed in Magma takes 1 . 50 seconds on a 2.2 6 GHz Intel Core 2 Duo , algorithm Secan tNumbers also takes 1 . 5 0 second s, and Atkin son’ s al- gorithm takes 4 . 51 second s. Thus, e ven if both T an gent and Secant numb ers are re- quired, Atkinson’ s algorithm is slightly slo wer . It also requires about twice as much memory . Acknowled gements W e thank Jon Borwein for encouraging the belief that high-precision com- putations are useful in “experimental” mathematics [4 ], e.g. in the PSLQ algorithm [16]. Ben F . “T ex” Loga n, Jr . [25] sugges ted t he use of T angent numbers to compute Bernou lli numbers. Chris- tian Reinsch [26] pointed out the numerical instability of the recurre nce (7) and sugges ted the use of the numerically stable recu rrence (8). Christopher Heckman kindly dre w our attention to Atkin- son’ s algorithm [2]. W e tha nk Paul Z immermann f or his comments . Some of the material presented here i s drawn from t he recent book Modern Computer Arithmetic [7] (and as-yet-unpublished so- lutions to exerc ises in the book). In particular , see [7, § 4.7.2 and ex ercises 4.35–4.41]. Finall y , w e thank David Bailey , Richard Crandall, and two anonymous referees for suggestions and pointers to additional refe rences. Fast Compu tation of Bernoulli, T angent and Secan t Numbers 15 Refer ences 1. Milton Abramo witz an d Irene A. S tegun, Handbook of Mathematical Functions , Dove r , 19 73. 2. M. D. Atkinson, How to compute the series expansions of sec x and tan x , Amer . Math. Monthly 93 (1986) , 387–389. 3. David H. Bailey , Jonathan M. Borwe in and Richard E. Crand all, On the Khintchine constant, Math. Comput. 66 (1997), 417–431. 4. Jonathan M. Borwein and R. M. Corless , Emer ging tools for experimental mathematics, Am. Math. Mon. 106 (1999), 899–909. 5. W ieb Bosma, John Cannon and Catherine Playoust, The M agma algebra system. I. The user language, J . Symbolic Comput. 24 (1997), 235-265 6. Richard P . Brent, Unrestricted algorithms for elemen tary and special functions, in Information Pr ocess ing 80 , North-Holland, Amsterdam, 1980, 613–619. arXi v:1004.3621v1 7. Richard P . Brent and Paul Zimmermann, Modern Computer Arithmetic , Cambridge Univ er- sity Press, 2010, 237 pp. arXiv :1004.4710v1 8. Joe Buhler , Richard Crandall, Reijo Ern v all, T auno Mets ¨ anky l ¨ a and M. Amin Shokrollahi, Irre gular primes and c yclotomic in v ariants to twelve milli on, J. Symbolic Computation 31 (2001) , 89–96. 9. Joe Buhler , Richard Crandall, Reijo Ern v all and T auno Mets ¨ ankyl ¨ a, Irre gular primes to four million, Math. Comput. 61 (1993), 151–153. 10. Joe Buhler , Richard Crandall and R. Sompolsk i, Irregular primes to one million, Math. Co m- put. 59 (1992), 717–722. 11. Joe Buhler and Dav id Harv ey , Irregu lar primes to 163 m illion, Math. Comput. 80 (201 1), 2435–244 4. 12. Thomas Cl ausen, Theorem, Astr on. Nac hr . 17 (1840), cols. 351–352. 13. Ri chard E. Crandall, T opics i n Advan ced Scientific Computation , Springer-V erlag, 1996. 14. Ri chard E. Crandall and Carl Po merance, Prime numbers : A Computational P erspective , Springer -V erlag, 2001. 15. Karl Dilcher and Ilja Sh . Sla vutskii, A Bibli ogr aphy of Berno ulli Numbe rs ( last updated March 3, 2007), http://ww w.mscs.da l.ca/%7Edilcher/bernoulli.html . 16. Helaman R. P . Fer guson, David H. Bail ey and Steve Arno , Analysis of PSLQ, an integer relation finding algorithm, Math. Comput. 68 (1999), 351–369. 17. M artin F ¨ urer , Fas ter integer multiplication, P r oc. 39th Annual A CM Sympo sium on Theory of Computing (STOC), A CM, San Di ego, Californ ia, 2007, 57–66. 18. Ronald L. Graham, Donald E. Knuth, and Oren Patashnik, Concr ete Mathematics , t hird edi- tion, Addison-W esley , 1994. 19. Ke vin Hare, Multisectioning, Rational P oly-Expone ntial Functions and P arallel Computa- tion , M.Sc. thesis, Dept. of Mathematics an d Statistics, Simon Fraser Univ ersity , Canada, 2002. 20. David Harv ey , A multimodular algorithm for computing Bernoulli numbers, Math. Comput. 79 (2010), 2361–2370 . 21. M asanob u Kanek o, The Ak iyama-T anigawa algorithm for Bernoulli numbers, J . of Inte ger Sequences 3 (2000) . Article 00.2.9, 6 pa ges. http:// www.cs.uw aterloo.ca/journals/JIS/ . 22. Donald E. Knuth, Euler’ s constant to 1271 places, Math. Comput. 16 (1962), 275–281. 23. Donald E. Knuth and Thomas J. Buckholtz, Compu tation of T angent, Euler , and Bernoulli numbers, Math. Comput. 21 (1967), 663–688. 24. H. T . Kung, On computing re ciprocals of power series, Numer . Math. 22 (1974) , 341–348. 25. B. F . Logan, unpublished observ ation, m entioned in [18, § 6.5]. 26. Christian Reinsch, per sonal communication t o R. P . Brent, about 1979, ack nowledg ed in [6]. 27. Arnold Sch ¨ onhage and V olker Stras sen, Schnelle Multipli kation großer Zahlen, Computing 7 (1971) , 281–292. 28. M alte Siev eking, An algorithm for div ision of power series , Computing 10 (1972), 153–156 . 29. Neil J. A. Sloane, The On-L ine Encyclopedia of Inte ger Sequences, http://oeis. org . 16 Richard P . Brent and David Har vey 30. Karl G. C. von Staudt, Be weis eines Lehrsatzes, die Bernoullischen Zahlen betref fend, J. Reine Ange w . Math. 21 (1840) , 372–37 4. http://gdz.su b.uni- goettingen.de . 31. Al lan Steel, Reduce e very thing t o multiplication, prese nted at Com- puting by the Number s: Algorithms, Pr ecision and Complex- ity , W orkshop for Richard Brent’ s 60th birthday , Berlin, 2006, http:// www.mathe matik.hu- berlin.de/%7Egaggle/EVENTS/2006/BRENT60/ . 32. William Stein et al , Sa ge, http:// www.sagem ath.org/ . 33. Edward C. Titchmars h, The Theory of the Riemann Zeta-Function , se cond edition (r evis ed by D. R. Heath-Bro wn), Clarendon Press, Oxford, 1986.