Compute the first few Bernoulli numbers.
Now time the computation of the one thousandth Bernoulli number.
> st := time(): b1000 := bernoulli(1000): t1 := time()-st;
Now use Euler's formula for Zeta(1000) .
> st := time():
> Digits := 2000:
> zeta_1000 := add(1./i^1000,i=1..100):
> B := zeta_1000*2*1000!/evalf((2*Pi)^1000):
> B := -(-1)^irem(1000/2,2)*B:
> Digits := Digits-12:
> B := evalf(B):
> B := convert(B,rational):
> t2 := time()-st;
Check the answer.
> zero_4 := b1000-B;
It was more than 10 times faster to use Euler's formula than it was to use the Maple procedure bernoulli(1000), but the procedure did compute all the even indexed Bernoulli numbers up to 1000.