30 lines
766 B
Text
30 lines
766 B
Text
include builtins/mpfr.e
|
|
procedure bernoulli(mpq rop, integer n)
|
|
sequence a = mpq_init(n+1)
|
|
for m=1 to n+1 do
|
|
mpq_set_si(a[m], 1, m)
|
|
for j=m-1 to 1 by -1 do
|
|
mpq_sub(a[j], a[j+1], a[j])
|
|
mpq_set_si(rop, j, 1)
|
|
mpq_mul(a[j], a[j], rop)
|
|
end for
|
|
end for
|
|
mpq_set(rop, a[1])
|
|
a = mpq_free(a)
|
|
end procedure
|
|
|
|
mpq rop = mpq_init()
|
|
mpz n = mpz_init(),
|
|
d = mpz_init()
|
|
for i=0 to 60 do
|
|
bernoulli(rop, i)
|
|
if mpq_cmp_si(rop, 0, 1) then
|
|
mpq_get_num(n, rop)
|
|
mpq_get_den(d, rop)
|
|
string ns = mpfr_sprintf("%44Zd",n),
|
|
ds = mpfr_sprintf("%Zd",d)
|
|
printf(1,"B(%2d) = %s / %s\n", {i,ns,ds})
|
|
end if
|
|
end for
|
|
{n,d} = mpz_free({n,d})
|
|
rop = mpq_free(rop)
|