Re: computing Bernoulli numbers

"Bart Vandewoestyne" <MyFirstName.MyLastName@xxxxxxxxxx> wrote in message
news:vJkoj.683$6n5.478@xxxxxxxxxxxxxxxx

I would like to implement a routine to compute floating point
Bernoulli numbers.

I have seen the formula that gets them as a linear combination of
preceding Bernoulli numbers quoted, and there is also one that
uses half the terms and a quadratic dependence on lower Bernoulli
numbers, maybe that was quoted but I just didn't read carefully
enough.

The Riemann zeta function approach has the advantage that big
Bernoulli numbers get big like factorials and in this method
you can encapsulate the part that wants to overflow in
ln(Gamma(x)). The problem I see with your method is that the
first omitted term is a poor approximation to the error of a
slowly converging series. Use integral from n-0.5 to oo of
f(x) dx for the approximation of sum([(f(k),k=n,oo)]) or the
Euler-MacLaurin summation formula.

If you want to use the first omitted term as an approximation
still, you can start with the alternating series for
(1-2.0/2**n)*zeta(n); there is also a simple refinement of this
method that improves the approximation rapidly. ISTR it's in
the chapter on sums and series in Scheid's Schaum's outline on
numerical analysis, at least the first edition.

Oh, also i**n overflows quickly. You need real(i,qp)**n or
similar.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end


.