faulhaber_s_formula.sf

#!/usr/bin/ruby

# The formula for calculating the sum of consecutive
# numbers raised to a given power, such as:
#    1^p + 2^p + 3^p + ... + n^p
# where p is a positive integer.

# See also:
#   https://en.wikipedia.org/wiki/Faulhaber%27s_formula

# The Faulhaber's formula
# See: https://en.wikipedia.org/wiki/Faulhaber%27s_formula
func faulhaber_formula (n, p) {
    sum(0..p, {|k|
        binomial(p + 1, k) * bernoulli(k) * n**(p + 1 - k)
    }) / (p+1)
}

# Alternate expression using Bernoulli polynomials
# See: https://en.wikipedia.org/wiki/Faulhaber%27s_formula#Alternate_expressions
func bernoulli_polynomials (n, x) {
    sum(0..n, {|k|
        binomial(n, k) * bernoulli(n - k) * x**k
    })
}

func faulhaber_formula_2 (n, p) {
    (bernoulli_polynomials(p + 1, n) - bernoulli(p + 1)) / (p + 1)
}

for p in (0..20) {

    var n = 100.irand
    var r = faulhaber_formula(n, p)

    assert_eq(r, faulhaber_formula_2(n, p))
    assert_eq(r, n.faulhaber_sum(p))            # built-in

    printf("faulhaber_#{p}(%2d) = %s\n", n, r)
}