Re: [Axiom-developer] Bernoulli puzzle

[Waldek Hebisch]

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
address@hidden