sidef/scripts/faulhaber_s_formula.sf

#!/usr/bin/ruby

#
## https://rosettacode.org/wiki/Faulhaber%27s_formula
#

func bernoulli({.is_one}) { 1/2 }
func bernoulli({.is_odd}) { 0/1 }

func bernoulli(n) {

    var a = []
    for m in ^(n+1) {
        a[m] = 1/(m + 1)
        for j in (m^..0 + 1) {
            a[j-1] = j*(a[j-1] - a[j])
        }
    }

    return a[0]
}

func faulhaber_s_formula(p) {

    var formula = gather {
        { |j|
            take "(#{binomial(p+1, j) * j.bernfrac -> as_rat})*n^#{p+1 - j}"
        } << 0..p
    }

    formula.grep! { !.contains('(0)*') }.join!(' + ')

    formula -= /\(1\)\*/g
    formula -= /\^1\b/g
    formula.gsub!(/\(([^+]*?)\)/, { _ })

    "1/#{p + 1} * (#{formula})"
}

{ |p|
    printf("%2d: %s\n", p, faulhaber_s_formula(p))
} << ^10