//
// 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

//