Tim Daly wrote: > > Axiom's bernoulli algorithm reads: <snip> > > I'm unable to decode the equations supporting the inner loop > computation (t and b) which support computing the B(i) term. > > Since the loop iterates by 2 it looks like someone worked out > the details of the computation taking 2 steps at a time. > Unfortunately I can't figure out how to unwind this into any > of the usual bernoulli formulas. > > Does anyone recognize the formula that is being used? > Can you translate it into TeX format?
Usualy Bernoulli numbers B_n are defined as coefficients of the following expansion: t/(exp(t) - 1) = \sum_{n=0}^\infty B_n t^n/n! Since t/(exp(t) - 1) + (1/2)*t = t(exp(t) + 1)/(2(exp(t) - 1) = t/(2tanh(t/2)) we have B_0 = 1, B_1 = -1/2 and because t/(2tanh(t/2)) is an even function for odd n different than 1 we have B_n = 0. Given this information we can rewrite the defining formula: t/(exp(t) - 1) = - (1/2)t + \sum_{n=0}^\infty B_{2n}t^{2n}/(2n}! Multiply both sides by (exp(t) - 1): t = - t(exp(t) - 1)/2 + (exp(t) - 1)\sum_{n=0}^\infty B_{2n}t^{2n}/(2n}! Compute coefficient before t^{2n+1} for both sides: 0 = -1/(2(2n)!) + \sum_{k=0}^{n}B_{2k}/((2k)!(2n+1 - 2k)!) so B_{2n}/(2n)! = -( -1/(2(2n)!) + \sum_{k=0}^{n-1}B_{2k}/((2k)!(2n+1 - 2k)!)) Multiplying by (2n+1)!: B_{2n}(2n+1) = -(-(2n+1)/2 + \sum_{k=0}^{n-1}B_{2k}binomial(2n+1, 2k)) Since B_0 = 1 and binomial(2n+1, 0) = 1 the first term of sum is 1 and we can write: B_{2n}(2n+1) = -((1 - 2n)/2 + \sum_{k=1}^{n-1}B_{2k}binomial(2n+1, 2k)) Now, the variable b is used to compute expression in the parenthesis. More precisely, table of Bernoulli numbers always has odd length, to variable l is odd and consequently l + 1 is even. So i goes trough even numbers, starting from first unknown position in the tables. In other words, i in the inner loop is 2n in the formula above. Consequently, initialization b := (1-i)/2 computes the (1 - 2n)/2 term before sum. Variable t is used to compote binomial(2n+1, 2k) using formula binomial(2n+1, 2k) = (2n + 1 - (2k - 1))(2n + 1 - (2k-2))/ ((2k)(2k-1))binomial(2n+1, 2k-2) and binomial(2n+1, 0) = 1. Since j in the inner loop is 2k from above assignment to t implements this formula. Once experession in parenthesis from formula for B_{2n} is computed the B(i) := -b/((i+1)::RN) line gets value of Bernoulli number. PS. This is essentially using the recursive definition given in Wikipedia. Wikipedia calls this inefficient, but AFAICS due to memoization it is much better than Akiyama-Tanigawa presented in Wikipedia. -- Waldek Hebisch hebi...@math.uni.wroc.pl _______________________________________________ Axiom-developer mailing list Axiom-developer@nongnu.org https://lists.nongnu.org/mailman/listinfo/axiom-developer