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
_______________________________________________
Axiom-developer mailing list
Axiom-developer@nongnu.org
https://lists.nongnu.org/mailman/listinfo/axiom-developer
