93 lines
1.7 KiB
Text
93 lines
1.7 KiB
Text
//
|
|
// Bernoulli Number
|
|
// Using FutureBasic 7.0.35
|
|
//
|
|
// September 2025, R.W
|
|
//
|
|
|
|
include "gmp.incl"
|
|
|
|
_kMaxN = 200 // max Bernoulli terms
|
|
|
|
mpq_t tmp
|
|
dim a( _kMaxN ) as mpq_t // workspace a(0.._kMaxN)
|
|
|
|
// Bernoulli:
|
|
// B_m = -(1 / (m+1)) * sum_{k=0..m-1} binomial(m+1, k) * B_k
|
|
|
|
local fn Bernoulli( rop as mpq_t, n as long )
|
|
long m, j
|
|
if n < 0 then n = 0
|
|
if n > _kMaxN then n = _kMaxN
|
|
for m = 1 to n + 1
|
|
mpq_set_si( a(m-1), 1, m )
|
|
// j = m-1..1
|
|
for j = m - 1 to 1 step -1
|
|
// a[j] = a[j+1] - a[j] with 0-based mapping
|
|
mpq_sub( a(j-1), a(j), a(j-1) )
|
|
// multiply by j
|
|
mpq_set_si( tmp, j, 1 )
|
|
mpq_mul( a(j-1), a(j-1), tmp )
|
|
next
|
|
next
|
|
mpq_set( rop, a(0) )
|
|
end fn
|
|
|
|
|
|
//-------------- Main ---------
|
|
|
|
Window 1, @"Bernoulli Numbers",(0,0,500,500)
|
|
|
|
CFStringRef cfnum, cfden, pLine, filler,sp,sp2
|
|
Int terms
|
|
long i
|
|
|
|
terms = 60
|
|
|
|
mpq_init( tmp )
|
|
for i = 0 to terms
|
|
mpq_init( a(i) )
|
|
next
|
|
filler = @" "
|
|
mpq_t rop
|
|
mpz_t num, den
|
|
|
|
mpq_init( rop )
|
|
mpz_init( num )
|
|
mpz_init( den )
|
|
pLine = @""
|
|
for i = 0 to 60 step 2
|
|
if i = 2 then i--
|
|
if i = 3 then i--
|
|
fn Bernoulli( rop, i )
|
|
mpq_get_num( num, rop )
|
|
if i == 1
|
|
mpz_neg( num, num )
|
|
end if
|
|
if fn mpz_cmp_si( num, 0 ) // nonzero?
|
|
mpq_get_den( den, rop )
|
|
cfnum = fn mpz_cf(10, num)
|
|
cfden = fn mpz_cf(10, den)
|
|
if i < 10
|
|
sp = @" "
|
|
else
|
|
sp = @""
|
|
end if
|
|
sp2 = left(filler,44 - len(cfnum))
|
|
pLine = concat(sp, @"B(",mid(str(i),1),@") = ",sp2 ,cfnum, @" / ", cfden)
|
|
print@ pLine
|
|
end if
|
|
next
|
|
|
|
// clean up
|
|
for i = 0 to terms
|
|
mpq_clear( a(i) )
|
|
next
|
|
mpq_clear( tmp )
|
|
mpz_clear( den )
|
|
mpz_clear( num )
|
|
mpq_clear( rop )
|
|
|
|
HandleEvents
|
|
|
|
//
|