blandin-diaz_compositional_bernoulli_numbers_B_S_1.sf

#!/usr/bin/ruby

# Author: Daniel "Trizen" Șuteu
# Date: 23 February 2018
# https://github.com/trizen

# A new recurrence for computing Blandin-Diaz compositional Bernoulli numbers (B^S)_1,n.

# Formula:
#  a(0) = 1
#  a(n) = -Sum_{k=0..n-1} a(k) / ((n-k+1)!)^2

# Which gives us the nth Blandin-Diaz compositional Bernoulli number as:
#   (B^S)_1,n = a(n) * n!

# See also:
#   https://arxiv.org/abs/0708.0809

# OEIS entries:
#   https://oeis.org/A133002   (numerators)
#   https://oeis.org/A133003   (denominators)

func a((0)) { 1 }

func a(n) is cached {
    -sum(^n, {|k| a(k) / (n - k + 1)!**2 })
}

for n in (0..30) {
    printf("(B^S)_1(%2d) = %45s / %s\n", n, a(n) * n! -> nude)
}

__END__
(B^S)_1( 0) =                                             1 / 1
(B^S)_1( 1) =                                            -1 / 4
(B^S)_1( 2) =                                             5 / 72
(B^S)_1( 3) =                                            -1 / 48
(B^S)_1( 4) =                                           139 / 21600
(B^S)_1( 5) =                                            -1 / 540
(B^S)_1( 6) =                                           859 / 2540160
(B^S)_1( 7) =                                            71 / 483840
(B^S)_1( 8) =                                         -9769 / 36288000
(B^S)_1( 9) =                                           233 / 896000
(B^S)_1(10) =                                      -6395527 / 31614105600
(B^S)_1(11) =                                        145069 / 1149603840
(B^S)_1(12) =                                 -304991568097 / 7139902049280000
(B^S)_1(13) =                                  -95164619917 / 2196892938240000
(B^S)_1(14) =                                  119780081383 / 941525544960000
(B^S)_1(15) =                                   -3046785293 / 15216574464000
(B^S)_1(16) =                              4002469707564917 / 16326052949606400000
(B^S)_1(17) =                              -102407337854027 / 443241256550400000
(B^S)_1(18) =                           1286572077762833639 / 11991344662654156800000
(B^S)_1(19) =                            219276930957009857 / 1100420292929126400000
(B^S)_1(20) =                      -20109624681057406222913 / 25964416811662737408000000
(B^S)_1(21) =                        1651690537394493957719 / 989120640444294758400000
(B^S)_1(22) =                     -317111791627190377990199 / 114672159963196932096000000
(B^S)_1(23) =                    11537868018533936870610343 / 3350421369359492972544000000
(B^S)_1(24) =                -53268794333233082810667038099 / 27744839359665961305636864000000
(B^S)_1(25) =               -131365008403523370365114156231 / 22195871487732769044509491200000
(B^S)_1(26) =                  5728528307727220295170552267 / 204884967579071714257010688000
(B^S)_1(27) =               -575909751690138690377372607797 / 7588332132558211639148544000000
(B^S)_1(28) =           36309996434261828839678688809299961 / 230769988036606221013588377600000000
(B^S)_1(29) =            -348737364778474586752969456387259 / 1446833780793769410743500800000000
(B^S)_1(30) =    128926955111338066596432445962308410078213 / 870537324556366868721681061505925120000000