[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Axiom-developer] Bernoulli puzzle

From: daly
Subject: [Axiom-developer] Bernoulli puzzle
Date: Mon, 20 Oct 2014 02:39:07 -0500

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
    -- 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)

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?


reply via email to

[Prev in Thread] Current Thread [Next in Thread]