[Axiom-developer] Bernoulli puzzle

[daly]

Axiom's bernoulli algorithm reads:

  bernoulli n ==
    n < 0 => error "bernoulli not defined for negative integers"
    -- except for 1, all odd terms are zero
    odd? n =>
      n = 1 => -1/2
      0
    -- B is a cache. If we already know the answer, just return it
    l := (#B) :: I
    n < l => B(n)
    -- extend B to hold the new terms we will compute
    -- there are n+1 terms but we may have cached l of them
    concat_!(B, new((n+1-l)::NNI, 0)$IndexedFlexibleArray(RN,0))
    -- compute B(i) i = l+2,l+4,...,n given B(j) j = 0,2,...,i-2
    -- 'by 2' to only compute even entries
    for i in l+1 .. n by 2 repeat
      t:I := 1
      b := (1-i)/2
      for j in 2 .. i-2 by 2 repeat
        t := (t*(i-j+2)*(i-j+3)) quo (j*(j-1))
        b := b + (t::RN) * B(j)
      B(i) := -b/((i+1)::RN)
    B(n)

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?

Tim